|
| 1 | +/** |
| 2 | + * Secular Frequencies |
| 3 | + * |
| 4 | + * This example integrates the outer Solar System and then performs a |
| 5 | + * frequency analysis using the Frequency Modified Fourier Transform |
| 6 | + * to determine the secular frequencies (g-modes). |
| 7 | + */ |
| 8 | +#include <stdio.h> |
| 9 | +#include <stdlib.h> |
| 10 | +#include <math.h> |
| 11 | +#include "rebound.h" |
| 12 | + |
| 13 | +double ss_pos[6][3] = |
| 14 | + { |
| 15 | + {-4.06428567034226e-3, -6.08813756435987e-3, -1.66162304225834e-6}, // Sun |
| 16 | + {+3.40546614227466e+0, +3.62978190075864e+0, +3.42386261766577e-2}, // Jupiter |
| 17 | + {+6.60801554403466e+0, +6.38084674585064e+0, -1.36145963724542e-1}, // Saturn |
| 18 | + {+1.11636331405597e+1, +1.60373479057256e+1, +3.61783279369958e-1}, // Uranus |
| 19 | + {-3.01777243405203e+1, +1.91155314998064e+0, -1.53887595621042e-1}, // Neptune |
| 20 | + {-2.13858977531573e+1, +3.20719104739886e+1, +2.49245689556096e+0} // Pluto |
| 21 | +}; |
| 22 | +double ss_vel[6][3] = |
| 23 | + { |
| 24 | + {+6.69048890636161e-6, -6.33922479583593e-6, -3.13202145590767e-9}, // Sun |
| 25 | + {-5.59797969310664e-3, +5.51815399480116e-3, -2.66711392865591e-6}, // Jupiter |
| 26 | + {-4.17354020307064e-3, +3.99723751748116e-3, +1.67206320571441e-5}, // Saturn |
| 27 | + {-3.25884806151064e-3, +2.06438412905916e-3, -2.17699042180559e-5}, // Uranus |
| 28 | + {-2.17471785045538e-4, -3.11361111025884e-3, +3.58344705491441e-5}, // Neptune |
| 29 | + {-1.76936577252484e-3, -2.06720938381724e-3, +6.58091931493844e-4} // Pluto |
| 30 | +}; |
| 31 | + |
| 32 | +double ss_mass[6] = |
| 33 | + { |
| 34 | + 1.00000597682, // Sun + inner planets |
| 35 | + 1. / 1047.355, // Jupiter |
| 36 | + 1. / 3501.6, // Saturn |
| 37 | + 1. / 22869., // Uranus |
| 38 | + 1. / 19314., // Neptune |
| 39 | + 0.0 // Pluto |
| 40 | +}; |
| 41 | + |
| 42 | +int main(int argc, char* argv[]) { |
| 43 | + struct reb_simulation* r = reb_simulation_create(); |
| 44 | + const double k = 0.01720209895; // Gaussian constant |
| 45 | + r->dt = 120; // Timestep is 120 days. |
| 46 | + r->G = k * k; // These are the same units as used by the mercury6 code. |
| 47 | + r->integrator = REB_INTEGRATOR_WHFAST; |
| 48 | + |
| 49 | + // Initial conditions |
| 50 | + for (int i = 0; i < 6; i++) { |
| 51 | + struct reb_particle p = {0}; |
| 52 | + p.x = ss_pos[i][0]; |
| 53 | + p.y = ss_pos[i][1]; |
| 54 | + p.z = ss_pos[i][2]; |
| 55 | + p.vx = ss_vel[i][0]; |
| 56 | + p.vy = ss_vel[i][1]; |
| 57 | + p.vz = ss_vel[i][2]; |
| 58 | + p.m = ss_mass[i]; |
| 59 | + reb_simulation_add(r, p); |
| 60 | + } |
| 61 | + |
| 62 | + reb_simulation_move_to_com(r); |
| 63 | + |
| 64 | + int Nsamples = 2048; // Number of samples. Must be a power of two. |
| 65 | + // Choose a larger number for better accuracy, e.g. 32768. |
| 66 | + double* inp = malloc(sizeof(double)*2*Nsamples); |
| 67 | + // Start integration |
| 68 | + for (int i=0; i<Nsamples; i++){ |
| 69 | + // Integrate for 1000 steps (120000 days) |
| 70 | + reb_simulation_steps(r, 1000); |
| 71 | + // Calculate orbital elements of Jupiter |
| 72 | + struct reb_orbit o = reb_orbit_from_particle(r->G, r->particles[1], r->particles[0]); |
| 73 | + // Store complex eccentricity in array |
| 74 | + inp[i*2+0] = o.e*cos(o.pomega); |
| 75 | + inp[i*2+1] = o.e*sin(o.pomega); |
| 76 | + } |
| 77 | + |
| 78 | + // Perform frequency analysis |
| 79 | + int nfreq = 5; |
| 80 | + double datasep = 120000.0/365.25*2.0*M_PI; // sampling interval in units of year/2pi |
| 81 | + double minfreq = 60.0/1296000.0*datasep; // min/max frequenxy 60"/year |
| 82 | + double* out = malloc(sizeof(double)*3*nfreq); |
| 83 | + // The next command performs the actual Frequency Modified Fourier Transform (FMFT). |
| 84 | + // Other options are MFT (faster) and FMFT2 (more accurate). |
| 85 | + // See Sidlichovsky and Nesvorny (1996) for more details: |
| 86 | + // https://ui.adsabs.harvard.edu/abs/1996CeMDA..65..137S/abstract |
| 87 | + int error = reb_frequency_analysis(out, nfreq, -minfreq, minfreq, REB_FREQUENCY_ANALYSIS_FMFT, inp, Nsamples); |
| 88 | + if (error){ |
| 89 | + printf("An error occured during the frequency analysis.\n"); |
| 90 | + } |
| 91 | + |
| 92 | + // Output the nfreq most dominate modes |
| 93 | + for (int i=0; i<nfreq; i++){ |
| 94 | + double nu = out[0*nfreq+i]*1296000.0/datasep; // frequency in "/year |
| 95 | + double A = out[1*nfreq+i]; // amplitude error |
| 96 | + double phi = out[2*nfreq+i]/M_PI*180.0; // phase in deg |
| 97 | + printf("nu = %5.2f\"/yr A = %0.6f phi = %5.1f°\n", nu, A, phi); |
| 98 | + } |
| 99 | + |
| 100 | + // Cleanup |
| 101 | + free(out); |
| 102 | + free(inp); |
| 103 | + reb_simulation_free(r); |
| 104 | +} |
| 105 | + |
0 commit comments