Merge branch 'arya'

dtamayo · Jun 20, 2023 · 8bef4e4 · 8bef4e4
2 parents d5215d4 + cd582cf
commit 8bef4e4
Show file tree

Hide file tree

Showing 9 changed files with 377 additions and 332 deletions.
diff --git a/.github/workflows/python.yml b/.github/workflows/python.yml
@@ -8,7 +8,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: ["2.7", "3.7", "3.8", "3.9", "3.10"]
+        python-version: ["3.7", "3.8", "3.9", "3.10"]
 
     steps:
       - uses: actions/checkout@v3

diff --git a/doc/effects.rst b/doc/effects.rst
@@ -199,6 +199,39 @@ em_afin (double)               Yes          Object's final semimajor axis
 General Relativity
 ^^^^^^^^^^^^^^^^^^
 
+.. _lense_thirring:
+
+lense_thirring
+**************
+
+======================= ===============================================
+Authors                 A. Akmal
+Implementation Paper    None
+Based on                `Park et al. <https://iopscience.iop.org/article/10.3847/1538-3881/abd414/>`_.
+C Example               :ref:`c_example_lense_thirring`
+Python Example          `LenseThirring.ipynb <https://github.com/dtamayo/reboundx/blob/master/ipython_examples/LenseThirring.ipynb>`_.
+======================= ===============================================
+
+Adds Lense-Thirring effect due to rotating central body in the simulation. Assumes the source body is particles[0]
+
+**Effect Parameters**
+
+============================ =========== ==================================================================
+Field (C type)               Required    Description
+============================ =========== ==================================================================
+lt_c (double)                Yes         Speed of light in the units used for the simulation.
+============================ =========== ==================================================================
+
+**Particle Parameters**
+
+============================ =========== ==================================================================
+Field (C type)               Required    Description
+============================ =========== ==================================================================
+I (double)                   Yes         Moment of Inertia of source body 
+Omega (reb_vec3d)            Yes         Angular rotation frequency (Omega_x, Omega_y, Omega_z) 
+============================ =========== ==================================================================
+
+
 .. _gr_potential:
 
 gr_potential
@@ -223,7 +256,7 @@ Nice if you have a single-star system, don't need to get GR exactly right, and w
 ============================ =========== ==================================================================
 Field (C type)               Required    Description
 ============================ =========== ==================================================================
-c (double)                   Yes         Speed of light in the units used for the simulation.
+c (double)                   Yes         Speed of light, needs to be specified in the units used for the simulation.
 ============================ =========== ==================================================================
 
 
@@ -251,7 +284,7 @@ Adding this effect to several bodies is NOT equivalent to using gr_full.
 ============================ =========== ==================================================================
 Field (C type)               Required    Description
 ============================ =========== ==================================================================
-c (double)                   Yes         Speed of light in the units used for the simulation.
+c (double)                   Yes         Speed of light, needs to be specified in the units used for the simulation.
 ============================ =========== ==================================================================
 
 
@@ -276,7 +309,7 @@ This algorithm incorporates the first-order post-newtonian effects from all bodi
 ============================ =========== ==================================================================
 Field (C type)               Required    Description
 ============================ =========== ==================================================================
-c (double)                   Yes         Speed of light in the units used for the simulation.
+c (double)                   Yes         Speed of light, needs to be specified in the units used for the simulation.
 ============================ =========== ==================================================================
 
 **Particle Parameters**
@@ -458,7 +491,7 @@ This effect consistently tracks both the spin and orbital evolution of bodies un
 In all cases, we need to set masses for all the particles that will feel these tidal forces. Particles with only mass are point particles.
 
 Particles are assumed to have structure (i.e - physical extent & distortion from spin) if the following parameters are set: physical radius particles[i].r, potential Love number of degree 2 k2 (Q/(1-Q) in Eggleton 1998), and the spin angular rotation frequency vector Omega.
-If we wish to evolve a body's spin components, the fully dimensional moment of inertia I must be set as well. If this parameter is not set, the spin components will be stationary.
+If we wish to evolve a body's spin components, the fully dimensional moment of inertia I must be set as well. If this parameter is not set, the spin components will be stationary. Note that if the body is a test particle, this is assumed to be the specific moment of inertia.
 Finally, if we wish to consider the effects of tides raised on a specific body, we must set the constant time lag tau as well.
 
 For spins that are synchronized with a circular orbit, the constant time lag can be related to the tidal quality factor Q as tau = 1/(2*n*tau), with n the orbital mean motion.
@@ -476,8 +509,8 @@ Field (C type)               Required    Description
 ============================ =========== ==================================================================
 particles[i].r (float)       Yes         Physical radius (required for contribution from tides raised on the body).
 k2 (float)                   Yes         Potential Love number of degree 2.
-Omega (reb_vec3d)            Yes         Angular rotation frequency
-I (float)                    No          Moment of inertia
+Omega (reb_vec3d)            Yes         Angular rotation frequency (Omega_x, Omega_y, Omega_z)
+I (float)                    No          Moment of inertia (for test particles, assumed to be the specific MoI I/m)
 tau (float)                  No          Constant time lag. If not set, defaults to 0
 ============================ =========== ==================================================================
 

diff --git a/examples/lense_thirring/problem.c b/examples/lense_thirring/problem.c
@@ -1,7 +1,7 @@
 /**
- * Adding gravitational harmonics (J2, J4) to particles
+ * Adding Lense-Thirring effect
  * 
- * This example shows how to add a J2 and J4 harmonic to particles.
+ * This example shows how to add the Lense-Thirring effect to a simulation.
  * If you have GLUT installed for the visualization, press 'w' and/or 'c' for a clearer view of
  * the whole orbit.
  */
@@ -14,38 +14,31 @@
 
 int main(int argc, char* argv[]){
     struct reb_simulation* sim = reb_create_simulation();
+    sim->G = 4*M_PI*M_PI; // units of AU, yr, Msun
     struct reb_particle star = {0};
     star.m     = 1.;
     reb_add(sim, star);
     double omega = 90.361036076; //solar rotation rate in rad/year
     double C_I = 0.06884; //solar moment of inertia prefactor
     double R_eq = 0.00465247264; //solar equatorial radius in AU
-    double p_x = 0.0; //x-comp. of unit spin-pole direction
-    double p_y = 0.0; //y-comp. of unit spin-pole direction
-    double p_z = 1.0; //z-comp. of unit spin-pole direction
-
 
     struct reb_particle planet = {0};  // add a planet on a circular orbit (with default units where G=1)
-    planet.x = 1.;
-    planet.vy = 1.;
-    reb_add(sim, planet);
+    const double mp = 1.7e-7;   // approximate values for mercury in units of Msun and AU
+    const double a = 0.39;
+    const double e = 0.21;
+    reb_add_fmt(sim, "m a e", mp, a, e);
 
     struct rebx_extras* rebx = rebx_attach(sim);  // first initialize rebx
     struct rebx_force* lense  = rebx_load_force(rebx, "lense_thirring"); // add our new force
     rebx_add_force(rebx, lense);
-
-    rebx_set_param_double(rebx, &sim->particles[0].ap, "lt_rot_rate", omega);
-    rebx_set_param_double(rebx, &sim->particles[0].ap, "lt_Mom_I_fac", C_I);
-    rebx_set_param_double(rebx, &sim->particles[0].ap, "lt_R_eq", R_eq);
-    rebx_set_param_double(rebx, &sim->particles[0].ap, "lt_p_hatx", p_x);
-    rebx_set_param_double(rebx, &sim->particles[0].ap, "lt_p_haty", p_y);
-    rebx_set_param_double(rebx, &sim->particles[0].ap, "lt_p_hatz", p_z);
+    rebx_set_param_vec3d(rebx, &sim->particles[0].ap, "Omega", (struct reb_vec3d){.x=0, .y=0, .z=omega});
+    rebx_set_param_double(rebx, &sim->particles[0].ap, "I", C_I*star.m*R_eq*R_eq);
 
     // Have to set speed of light in right units (set by G & initial conditions).  Here we use default units of AU/(yr/2pi)
-    rebx_set_param_double(rebx, &lense->ap, "lt_c", 10065.32);
-
-
+    rebx_set_param_double(rebx, &lense->ap, "lt_c", 63241.077); // speed of light in AU/yr
 
-    double tmax = 100000.;
+    double tmax = 1000.;
     reb_integrate(sim, tmax);
+    struct reb_orbit orb = reb_tools_particle_to_orbit(sim->G, sim->particles[1], sim->particles[0]);
+    printf("Pericenter precession rate = %.3f arcsec/century\n", orb.pomega*206265/10.);
 }
diff --git a/ipython_examples/LT.ipynb b/ipython_examples/LT.ipynb
diff --git a/ipython_examples/LenseThirring.ipynb b/ipython_examples/LenseThirring.ipynb
diff --git a/ipython_examples/TypeIMigration.ipynb b/ipython_examples/TypeIMigration.ipynb
diff --git a/ipython_examples/YarkovskyEffect.ipynb b/ipython_examples/YarkovskyEffect.ipynb
diff --git a/src/lense_thirring.c b/src/lense_thirring.c
@@ -31,30 +31,26 @@
  * Implementation Paper    None
  * Based on                `Park et al. <https://iopscience.iop.org/article/10.3847/1538-3881/abd414/>`_.
  * C Example               :ref:`c_example_lense_thirring`
- * Python Example          `LT.ipynb <https://github.com/dtamayo/reboundx/blob/master/ipython_examples/LT.ipynb>`_.
+ * Python Example          `LenseThirring.ipynb <https://github.com/dtamayo/reboundx/blob/master/ipython_examples/LenseThirring.ipynb>`_.
  * ======================= ===============================================
  * 
- * Adds Lense-Thirring effect due to massive rotating bodies in the simulation.
+ * Adds Lense-Thirring effect due to rotating central body in the simulation. Assumes the source body is particles[0]
  *
  * **Effect Parameters**
  * 
  * ============================ =========== ==================================================================
  * Field (C type)               Required    Description
  * ============================ =========== ==================================================================
- * lt_c (double)                   Yes         Speed of light in the units used for the simulation.
+ * lt_c (double)                Yes         Speed of light in the units used for the simulation.
  * ============================ =========== ==================================================================
  *
  * **Particle Parameters**
  *
  * ============================ =========== ==================================================================
  * Field (C type)               Required    Description
  * ============================ =========== ==================================================================
- * lt_rot_rate (double)         No          rotation rate, omega`
- * lt_R_eq (double)             No          Equatorial radius of source body
- * lt_Mom_I_fac (double)        No          Moment of Inertia of source body over MR^2
- * lt_p_hatx (double)           No          x-component of spin-pole unit vector
- * lt_p_haty (double)           No          y-component of spin-pole unit vector
- * lt_p_hatz (double)           No          z-component of spin-pole unit vector
+ * I (double)                   Yes         Moment of Inertia of source body 
+ * Omega (reb_vec3d)            Yes         Angular rotation frequency (Omega_x, Omega_y, Omega_z) 
  * ============================ =========== ==================================================================
  * 
  */
@@ -66,14 +62,11 @@
 #include "rebound.h"
 #include "reboundx.h"
 
-static void rebx_calculate_LT_force(struct reb_simulation* const sim, struct reb_particle* const particles, const int N, const double omega, const double R_eq, const double C_fac, const double p_hat_x, const double p_hat_y, const double p_hat_z, const int source_index, const double C2){
-    const struct reb_particle source = particles[source_index];
+static void rebx_calculate_LT_force(struct reb_simulation* const sim, struct reb_particle* const particles, const int N, const struct reb_vec3d Omega, const double I, const double C2){
     const double G = sim->G;
     const double gamma = 1.000021;   //hard-coded Eddington-Robertson-Shiff parameter for now
-    for (int i=0; i<N; i++){
-        if(i == source_index){
-            continue;
-        }
+    const struct reb_particle source = particles[0]; // hard-code particles[0] as source particle
+    for (int i=1; i<N; i++){
         const struct reb_particle p = particles[i];
         const double dx = p.x - source.x;
         const double dy = p.y - source.y;
@@ -84,20 +77,25 @@ static void rebx_calculate_LT_force(struct reb_simulation* const sim, struct reb
         const double dvx = p.vx - source.vx;
         const double dvy = p.vy - source.vy;
         const double dvz = p.vz - source.vz;
-        const double Jx = C_fac*source.m * R_eq*R_eq*omega*p_hat_x ;
-        const double Jy = C_fac*source.m * R_eq*R_eq*omega*p_hat_y ;
-        const double Jz = C_fac*source.m * R_eq*R_eq*omega*p_hat_z ;
-        const double Omega_fac = (1.+gamma)*G/2/C2;
-        const double Omega_x = Omega_fac*(-Jx +3.*(Jx*dx+Jy*dy+Jz*dz)*dx/r2)/r3;
-        const double Omega_y = Omega_fac*(-Jy +3.*(Jx*dx+Jy*dy+Jz*dz)*dy/r2)/r3;
-        const double Omega_z = Omega_fac*(-Jz +3.*(Jx*dx+Jy*dy+Jz*dz)*dz/r2)/r3;
+        const double Jx = I*Omega.x; //C_fac*source.m * R_eq*R_eq*omega*p_hat_x ;
+        const double Jy = I*Omega.y; //C_fac*source.m * R_eq*R_eq*omega*p_hat_y ;
+        const double Jz = I*Omega.z; //C_fac*source.m * R_eq*R_eq*omega*p_hat_z ;
+        const double ms = source.m;
+        const double mt = p.m;
+        const double mtot = ms + mt;
+        const double mratio = mt/ms;
+        const double Omega_fac = ms/mtot*(1.+gamma)*G/2/C2/r3;
+        const double Jdotr = Jx*dx+Jy*dy+Jz*dz;
+        const double Omega_x = Omega_fac*(-Jx +3.*Jdotr*dx/r2);
+        const double Omega_y = Omega_fac*(-Jy +3.*Jdotr*dy/r2);
+        const double Omega_z = Omega_fac*(-Jz +3.*Jdotr*dz/r2);
 
-        particles[i].ax += 2.*Omega_y*dvz - Omega_z*dvy;
-        particles[i].ay += 2.*Omega_z*dvx - Omega_x*dvz;
-        particles[i].az += 2.*Omega_x*dvy - Omega_y*dvx;
-        particles[source_index].ax -= 2.*Omega_y*dvz - Omega_z*dvy; 
-        particles[source_index].ay -= 2.*Omega_z*dvx - Omega_x*dvz;
-        particles[source_index].az -= 2.*Omega_x*dvy - Omega_y*dvx;
+        particles[i].ax += 2.*(Omega_y*dvz - Omega_z*dvy);
+        particles[i].ay += 2.*(Omega_z*dvx - Omega_x*dvz);
+        particles[i].az += 2.*(Omega_x*dvy - Omega_y*dvx);
+        particles[0].ax -= mratio * 2. * (Omega_y*dvz - Omega_z*dvy); 
+        particles[0].ay -= mratio * 2. * (Omega_z*dvx - Omega_x*dvz);
+        particles[0].az -= mratio * 2. * (Omega_x*dvy - Omega_y*dvx);
     }
 }
 
@@ -109,25 +107,11 @@ void rebx_lense_thirring(struct reb_simulation* const sim, struct rebx_force* co
     }
     const double C2 = (*c)*(*c);
 
-    for (int i=0; i<N; i++){
-        const double* omega = rebx_get_param(rebx, particles[i].ap, "lt_rot_rate");
-        if(omega != NULL){
-           const double* R_eq = rebx_get_param(rebx, particles[i].ap, "lt_R_eq");
-           if (R_eq != NULL){
-               const double* C_fac = rebx_get_param(rebx, particles[i].ap, "lt_Mom_I_fac");
-               if(C_fac != NULL){
-                  const double* p_hat_x = rebx_get_param(rebx, particles[i].ap, "lt_p_hatx");
-                  if(p_hat_x != NULL){
-                     const double* p_hat_y = rebx_get_param(rebx, particles[i].ap, "lt_p_haty");
-                     if(p_hat_y != NULL){
-                        const double* p_hat_z = rebx_get_param(rebx, particles[i].ap, "lt_p_hatz");
-                        if(p_hat_z != NULL){
-                           rebx_calculate_LT_force(sim, particles, N, *omega, *R_eq, *C_fac, *p_hat_x, *p_hat_y, *p_hat_z, i, C2);
-                        }
-                     }
-                  }
-               }
-           }
+    const double* I = rebx_get_param(rebx, particles[0].ap, "I");
+    if (I != NULL){
+        const struct reb_vec3d* Omega  = rebx_get_param(rebx, particles[0].ap, "Omega");
+        if(Omega != NULL){
+            rebx_calculate_LT_force(sim, particles, N, *Omega, *I, C2);
         }
     }
 }

diff --git a/src/tides_spin.c b/src/tides_spin.c
@@ -57,7 +57,7 @@
  * ============================ =========== ==================================================================
  * particles[i].r (float)       Yes         Physical radius (required for contribution from tides raised on the body).
  * k2 (float)                   Yes         Potential Love number of degree 2.
- * Omega (reb_vec3d)            Yes         Angular rotation frequency
+ * Omega (reb_vec3d)            Yes         Angular rotation frequency (Omega_x, Omega_y, Omega_z)
  * I (float)                    No          Moment of inertia (for test particles, assumed to be the specific MoI I/m)
  * tau (float)                  No          Constant time lag. If not set, defaults to 0
  * ============================ =========== ==================================================================