Merge branch 'phoebe' of https://github.com/dtamayo/reboundx into dta…

…mayo-phoebe
dtamayo · Aug 1, 2024 · 6e23601 · 6e23601
2 parents 4cc397a + ae2e967
commit 6e23601
Show file tree

Hide file tree

Showing 3 changed files with 65 additions and 33 deletions.
diff --git a/ipython_examples/GasDampingTimescale.ipynb b/ipython_examples/GasDampingTimescale.ipynb
diff --git a/src/core.c b/src/core.c
@@ -63,6 +63,8 @@ void rebx_register_default_params(struct rebx_extras* rebx){
     rebx_register_param(rebx, "coordinates", REBX_TYPE_INT);
     rebx_register_param(rebx, "p", REBX_TYPE_DOUBLE);
     rebx_register_param(rebx, "d_factor", REBX_TYPE_DOUBLE);
+    rebx_register_param(rebx, "cs_coeff", REBX_TYPE_DOUBLE);
+    rebx_register_param(rebx, "d_coeff", REBX_TYPE_DOUBLE);
     rebx_register_param(rebx, "tau_a", REBX_TYPE_DOUBLE);
     rebx_register_param(rebx, "tau_e", REBX_TYPE_DOUBLE);
     rebx_register_param(rebx, "tau_inc", REBX_TYPE_DOUBLE);

diff --git a/src/gas_damping_timescale.c b/src/gas_damping_timescale.c
@@ -67,9 +67,15 @@ static struct reb_vec3d rebx_calculate_gas_damping_timescale(struct reb_simulati
     struct reb_orbit o = reb_orbit_from_particle(sim->G, *planet, *star);
 
     const double* const d_factor = rebx_get_param(rebx, planet->ap, "d_factor");
+    const double* const cs_coeff = rebx_get_param(rebx, force->ap, "cs_coeff");
+    const double* const d_coeff = rebx_get_param(rebx, force->ap, "d_coeff");
+
+    struct reb_vec3d a = {0};
 
-    // initialize damping timescales
-    double invtau_a = 0.0;
+    if (d_factor == NULL || cs_coeff == NULL || d_coeff == NULL){
+        rebx_error(rebx, "Need to set d_factor, cs_coeff, d_coeff parameters.  See examples in documentation.\n");
+        return a;
+    }
 
     // initialize positions and velocities
     const double dvx = planet->vx - star->vx;
@@ -88,33 +94,28 @@ static struct reb_vec3d rebx_calculate_gas_damping_timescale(struct reb_simulati
     const double planetMass = planet->m;
 
     // eccentricity and inclination timescales from Dawson+16 Eqn 16
-    const double cs_coeff = 0.272125; // 1.29 km s^-1 in AU yr^-1
     double coeff;
 
-    double v = sqrt(pow(e0, 2.)+pow(inc0, 2.))*pow(starMass, 1./2.)*pow(a0, -1./2.);
-    double cs = cs_coeff/(2.*M_PI)*pow(a0, -1./4.);
+    double vk = sqrt(sim->G*starMass/a0);
+    double v = sqrt(e0*e0+inc0*inc0)*vk;
+    double cs = *cs_coeff/sqrt(sqrt(a0));
     double v_over_cs = v/cs;
 
     if (v <= cs){
         coeff = 1.;
     }
     else {
-        if (inc0 < cs/v) {
-            coeff = pow(v_over_cs, 3.);
+        if (inc0 < cs/vk) {
+            coeff = v_over_cs*v_over_cs*v_over_cs;
         }
         else {
-            coeff = pow(v_over_cs, 4.);
+            coeff = v_over_cs*v_over_cs*v_over_cs*v_over_cs;
         }
     }
 
-    double tau_e = -0.003*(*d_factor)*pow(a0, 2.)*(starMass/planetMass)*2.*M_PI*coeff;
+    double tau_e = -*d_coeff*(*d_factor)*a0*a0*(starMass/planetMass)*coeff;
     double tau_inc = 2.*tau_e;  // from Kominami & Ida 2002 [Eqs. 2.9 and 2.10]
 
-    struct reb_vec3d a = {0};
-
-    a.x =  dvx*invtau_a/(2.);
-    a.y =  dvy*invtau_a/(2.);
-    a.z =  dvz*invtau_a/(2.);
 
     if (tau_e < INFINITY || tau_inc < INFINITY){
         const double vdotr = dx*dvx + dy*dvy + dz*dvz;
@@ -135,4 +136,4 @@ void rebx_gas_damping_timescale(struct reb_simulation* const sim, struct rebx_fo
     const int back_reactions_inclusive = 1;
     const char* reference_name = "primary";
     rebx_com_force(sim, force, coordinates, back_reactions_inclusive, reference_name, rebx_calculate_gas_damping_timescale, particles, N);
-}
+}