diff --git a/examples/hatp11_scattering/problem.c b/examples/hatp11_scattering/problem.c new file mode 100644 index 00000000..77c9741a --- /dev/null +++ b/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 + #include + #include + #include + #include + #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"); + } + */ +}