problem file

examples/hatp11_scattering/problem.c

@@ -0,0 +1,229 @@
+/**
+ * Kozai cycles
+ *
+ * This example uses the IAS15 integrator to simulate
+ * a Lidov Kozai cycle of a planet perturbed by a distant star.
+ * The integrator automatically adjusts the timestep so that
+ * even very high eccentricity encounters are resolved with high
+ * accuracy.
+ *
+ * This example includes self-consistent spin, tidal & dynamical effects
+ * as well as general relativity
+ */
+ #include <stdio.h>
+ #include <stdlib.h>
+ #include <string.h>
+ #include <unistd.h>
+ #include <math.h>
+ #include "rebound.h"
+ #include "reboundx.h"
+ #include "tides_spin.c"
+
+void heartbeat(struct reb_simulation* r);
+double tmax = 5e4 * 2 * M_PI;
+double DELTA = 3.;
+
+double obl(struct reb_vec3d v1, struct reb_vec3d v2){
+  return acos(reb_vec3d_dot(v1,v2) / (sqrt(reb_vec3d_length_squared(v1)) * sqrt(reb_vec3d_length_squared(v2))));
+}
+
+char orb_title[100] = "output_orb_3equals";
+//char spin_title[100] = "output_spin_3equals";
+
+int main(int argc, char* argv[]){
+
+    int index = 0;
+    if (argc >= 2){
+      strcat(orb_title, argv[1]);
+      //strcat(spin_title, argv[1]);
+    }
+
+    struct reb_simulation* sim = reb_create_simulation();
+    // Setup constants
+    sim->integrator        = REB_INTEGRATOR_IAS15;
+    //sim->heartbeat        = heartbeat;
+
+    // Initial conditions
+
+    struct reb_particle star = {0};
+    star.m  = 0.811;
+    star.r = 0.76 * 0.00465;
+    reb_add(sim, star);
+
+    // All planet masses and eccs the same
+    double pm = 3.60 * 9.55e-4;
+    double pr = 0.3 * 4.676e-4;
+    double pe = 0.05;
+
+    // c
+    double ci = 1.0 * (M_PI / 180.);
+    double ca = 8.68624;
+    double com = reb_random_uniform(sim, 0., 2.*M_PI);
+    double cOm = reb_random_uniform(sim, 0., 2.*M_PI);
+    double cf = reb_random_uniform(sim, 0., 2.*M_PI);
+
+    // d
+    double di = 2.0 * (M_PI / 180.);
+    double da = 13.367;
+    double dom = reb_random_uniform(sim, 0., 2.*M_PI);
+    double dOm = reb_random_uniform(sim, 0., 2.*M_PI);
+    double df = reb_random_uniform(sim, 0., 2.*M_PI);
+
+    // e
+    double ei = 3.0 * (M_PI / 180.);
+    double ea = 20.5702;
+    double eom = reb_random_uniform(sim, 0., 2.*M_PI);
+    double eOm = reb_random_uniform(sim, 0., 2.*M_PI);
+    double ef = reb_random_uniform(sim, 0., 2.*M_PI);
+
+    reb_add_fmt(sim, "m a e inc omega Omega f", pm, ca, pe, ci, com, cOm, cf);
+    reb_add_fmt(sim, "m a e inc omega Omega f", pm, da, pe, di, dom, dOm, df);
+    reb_add_fmt(sim, "m a e inc omega Omega f", pm, ea, pe, ei, eom, eOm, ef);
+
+    //struct reb_orbit orb = reb_tools_particle_to_orbit(sim->G, sim->particles[1], sim->particles[0]);
+    //sim->dt = orb.P/15.12345;     // initial timestep
+
+    // Only care about spin vector of the star?
+    struct rebx_extras* rebx = rebx_attach(sim);
+    struct rebx_force* effect = rebx_load_force(rebx, "tides_spin");
+    rebx_add_force(rebx, effect);
+    // Sun
+    const double solar_spin_period = 20 * 2. * M_PI / 365.;
+    const double solar_spin = (2 * M_PI) / solar_spin_period;
+    const double solar_k2 = 0.028;
+    // const double solar_tau = (0.2 / solar_k2) * (1 / 86400.) * 2 * M_PI / 365.; // seconds to reb years
+    rebx_set_param_double(rebx, &sim->particles[0].ap, "k2", solar_k2);
+    rebx_set_param_double(rebx, &sim->particles[0].ap, "I", 0.08 * star.m * star.r * star.r);
+    rebx_set_param_vec3d(rebx, &sim->particles[0].ap, "Omega", (struct reb_vec3d){.z=solar_spin});
+
+    //double solar_Q = 1e6;
+    //struct reb_orbit orb = reb_tools_particle_to_orbit(sim->G, sim->particles[1], sim->particles[0]);
+    //double solar_tau = 1 / (2 * solar_Q * orb.n);
+    // rebx_set_param_double(rebx, &sim->particles[0].ap, "tau", solar_tau);
+
+    // Spin vector of P1 just for fun...
+    /*
+    const double planet_k2 = 0.3;
+    const double planet_Q = 100.;
+    rebx_set_param_double(rebx, &sim->particles[1].ap, "k2", planet_k2);
+    rebx_set_param_double(rebx, &sim->particles[1].ap, "I", 0.25 * pm * pr * pr);
+
+    // tidally locked & 0 obliquity
+    const double spin_period = orb.P;
+    const double spin_rate = (2. * M_PI / spin_period);
+
+    struct reb_vec3d norm = orb.hvec;
+    struct reb_vec3d nhat = reb_vec3d_normalize(norm);
+    struct reb_vec3d Omega_p = reb_vec3d_mul(nhat, spin_rate);
+
+    rebx_set_param_vec3d(rebx, &sim->particles[1].ap, "Omega", Omega_p);
+    double planet_tau = 1. / (2. * planet_Q * orb.n);
+    rebx_set_param_double(rebx, &sim->particles[1].ap, "tau", planet_tau);
+    */
+
+    reb_move_to_com(sim);
+    struct reb_vec3d newz = reb_vec3d_add(reb_tools_angular_momentum(sim), rebx_tools_spin_angular_momentum(rebx));
+    struct reb_vec3d newx = reb_vec3d_cross((struct reb_vec3d){.z =1}, newz);
+    struct reb_rotation rot = reb_rotation_init_to_new_axes(newz, newx);
+    if isnan(rot.r) {
+      rot = reb_rotation_identity();
+    }
+    rebx_simulation_irotate(rebx, rot); // This rotates our simulation into the invariable plane aligned with the total ang. momentum (including spin)
+
+    rebx_spin_initialize_ode(rebx, effect);
+
+    FILE* forb = fopen(orb_title,"w");
+    //FILE* fspin = fopen(spin_title,"w");
+    fprintf(forb, "#Iconds: %e %e %e %e %e %e %e %e %e\n", com, cOm, cf, dom, dOm, df, ei, eom, eOm, ef);
+    fprintf(forb, "t,a1,i1,e1,p_ob,a2,i2,e2,pert_ob,a3,i3,e3,scat_ob,E\n");
+    //fprintf(fspin, "t,star_sx,star_sy,star_sz,magstar\n");
+    fclose(forb);
+    //fclose(fspin);
+
+    for (int i=0; i<1000; i++){
+        struct reb_particle* sun = &sim->particles[0];
+        struct reb_particle* p1 = &sim->particles[1];
+        struct reb_particle* pert = &sim->particles[2];
+        struct reb_particle* scat = &sim->particles[3];
+
+        struct reb_vec3d* Omega_sun = rebx_get_param(rebx, sun->ap, "Omega");
+        double star_sx = Omega_sun->x;
+        double star_sy = Omega_sun->y;
+        double star_sz = Omega_sun->z;
+        double magstar = sqrt(star_sx * star_sx + star_sy * star_sy + star_sz * star_sz);
+
+        //struct reb_vec3d* Omega_p = rebx_get_param(rebx, p1->ap, "Omega");
+        //double s1x = Omega_p->x;
+        //double s1y = Omega_p->y;
+        //double s1z = Omega_p->z;
+        //double mag1 = sqrt(s1x * s1x +s1y * s1y + s1z * s1z);
+
+        struct reb_orbit o1 = reb_tools_particle_to_orbit(sim->G, *p1, *sun);
+        double a1 = o1.a;
+        double Om1 = o1.Omega;
+        double i1 = o1.inc;
+        double pom1 = o1.pomega;
+        double e1 = o1.e;
+        double f1 = o1.f;
+        struct reb_vec3d norm1 = o1.hvec;
+
+        // Interpret planet spin vector in rotating frame
+        // Transform spin vector into planet frame, w/ z-axis aligned with orbit normal and x-axis aligned with line of nodes
+        //struct reb_vec3d line_of_nodes = reb_vec3d_cross((struct reb_vec3d){.z =1}, norm1);
+        //struct reb_rotation rot = reb_rotation_init_to_new_axes(norm1, line_of_nodes); // Arguments to this function are the new z and x axes
+        //if isnan(rot.r) {
+        //  rot = reb_rotation_identity();
+        //}
+        //struct reb_vec3d srot = reb_vec3d_rotate(*Omega_p, rot); // spin vector in the planet's frame
+
+        // Interpret the spin axis in the more natural spherical coordinates
+        //double mag_1;
+        //double theta_1;
+        //double phi_1;
+        //reb_tools_xyz_to_spherical(srot, &mag_1, &theta_1, &phi_1);
+
+        // Perturber orbital elements
+        struct reb_particle com = reb_get_com_of_pair(sim->particles[0],sim->particles[1]);
+        struct reb_orbit o2 = reb_tools_particle_to_orbit(sim->G, *pert, com);
+        double a2 = o2.a;
+        double Om2 = o2.Omega;
+        double i2 = o2.inc;
+        double pom2 = o2.pomega;
+        double e2 = o2.e;
+        struct reb_vec3d norm2 = o2.hvec;
+
+        struct reb_particle com2 = reb_get_com_of_pair(com,sim->particles[2]);
+        struct reb_orbit o3 = reb_tools_particle_to_orbit(sim->G, *scat, com2);
+        double a3 = o3.a;
+        double i3 = o3.inc;
+        double e3 = o3.e;
+        struct reb_vec3d norm3 = o3.hvec;
+
+        // Stellar obliquities
+        double p_ob = obl(*Omega_sun, norm1);
+        double pert_ob = obl(*Omega_sun, norm2);
+        double scat_ob = obl(*Omega_sun, norm3);
+        //ouble mi = obl(norm1 , norm2);
+
+        FILE* forb = fopen(orb_title,"a");
+        //FILE* fspin = fopen(spin_title,"a");
+        fprintf(forb, "%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f,%f\n", sim->t,a1,i1,e1,p_ob,a2,i2,e2,pert_ob,a3,i3,e3,scat_ob,reb_tools_energy(sim));
+        //fprintf(fspin, "%e,%e,%e,%e,%e\n", sim->t, star_sx, star_sy, star_sz, magstar);
+        fclose(forb);
+        //fclose(fspin);
+        reb_integrate(sim, sim->t+(1000 * 2 * M_PI));
+    }
+   rebx_free(rebx);
+   reb_free_simulation(sim);
+}
+
+void heartbeat(struct reb_simulation* sim){
+    if(reb_output_check(sim, 100.*M_PI)){        // outputs to the screen
+        reb_output_timing(sim, tmax);
+    }
+    /*
+    if(reb_output_check(r, 12.)){            // outputs to a file
+        reb_output_orbits(r, "orbits.txt");
+    }
+    */
+}