2body

examples/hatp11_scattering/problem.c

@@ -27,7 +27,7 @@ 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 orb_title[100] = "output_orb_2equals";
 //char spin_title[100] = "output_spin_3equals";
 
 int main(int argc, char* argv[]){
@@ -57,28 +57,30 @@ int main(int argc, char* argv[]){
 
     // c
     double ci = 1.0 * (M_PI / 180.);
-    double ca = 8.68624;
+    double ca = 6.91607;
     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 da = 10.643;
     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);
+    //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
@@ -134,8 +136,10 @@ int main(int argc, char* argv[]){
 
     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(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(forb, "#Iconds: %e %e %e %e %e %e\n", com, cOm, cf, dom, dOm, df);
+    fprintf(forb, "t,a1,i1,e1,p_ob,a2,i2,e2,pert_ob,E\n");
     //fprintf(fspin, "t,star_sx,star_sy,star_sz,magstar\n");
     fclose(forb);
     //fclose(fspin);
@@ -144,7 +148,7 @@ int main(int argc, char* argv[]){
         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_particle* scat = &sim->particles[3];
 
         struct reb_vec3d* Omega_sun = rebx_get_param(rebx, sun->ap, "Omega");
         double star_sx = Omega_sun->x;
@@ -191,23 +195,25 @@ int main(int argc, char* argv[]){
         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);
+        //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(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(forb, "%f,%f,%f,%f,%f,%f,%f,%f,%f,%f\n", sim->t,a1,i1,e1,p_ob,a2,i2,e2,pert_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);