Implementation of Gas Damping Forces

examples/gas_damping_timescale/Makefile

@@ -0,0 +1,48 @@
+export OPENGL=0
+
+ifndef REB_DIR
+ifneq ($(wildcard ../../../rebound/.*),) # Check for REBOUND in default location
+REB_DIR=../../../rebound
+endif
+ifneq ($(wildcard ../../../../rebound/.*),) # Check for REBOUNDx being inside REBOUND directory
+REB_DIR=../../../
+endif
+endif
+ifndef REB_DIR # REBOUND is not in default location and REB_DIR is not set
+    $(error REBOUNDx not in the same directory as REBOUND.  To use a custom location, you Must set the REB_DIR environment variable for the path to your rebound directory, e.g., export REB_DIR=/Users/dtamayo/rebound.  See reboundx.readthedocs.org)
+endif
+PROBLEMDIR=$(shell basename `dirname \`pwd\``)"/"$(shell basename `pwd`)
+
+include $(REB_DIR)/src/Makefile.defs
+
+REBX_DIR=../../
+
+all: librebound.so libreboundx.so
+	@echo ""
+	@echo "Compiling problem file ..."
+	$(CC) -I$(REBX_DIR)/src/ -I$(REB_DIR)/src/ -Wl,-rpath,./ $(OPT) $(PREDEF) problem.c -L. -lreboundx -lrebound $(LIB) -o rebound
+	@echo ""
+	@echo "Problem file compiled successfully."
+
+librebound.so:
+	@echo "Compiling shared library librebound.so ..."
+	$(MAKE) -C $(REB_DIR)/src/
+	@echo "Creating link for shared library librebound.so ..."
+	@-rm -f librebound.so
+	@ln -s $(REB_DIR)/src/librebound.so .
+
+libreboundx.so: 
+	@echo "Compiling shared library libreboundx.so ..."
+	$(MAKE) -C $(REBX_DIR)/src/
+	@-rm -f libreboundx.so
+	@ln -s $(REBX_DIR)/src/libreboundx.so .
+
+clean:
+	@echo "Cleaning up shared library librebound.so ..."
+	@-rm -f librebound.so
+	$(MAKE) -C $(REB_DIR)/src/ clean
+	@echo "Cleaning up shared library libreboundx.so ..."
+	@-rm -f libreboundx.so
+	$(MAKE) -C $(REBX_DIR)/src/ clean
+	@echo "Cleaning up local directory ..."
+	@-rm -vf rebound

examples/gas_damping_timescale/problem.c

@@ -0,0 +1,74 @@
+/**
+ * Gas Damping Timescale Example
+ *
+ * This example shows how to add gas damping to a REBOUND simulation.
+ *
+ */
+#include <stdio.h>
+#include <stdlib.h>
+#include <unistd.h>
+#include <math.h>
+#include <string.h>
+#include "rebound.h"
+#include "reboundx.h"
+
+void heartbeat(struct reb_simulation* sim);
+
+int main(int argc, char* argv[]){
+
+    /* start the rebound simulation here */
+    struct reb_simulation* sim = reb_simulation_create();
+
+    // set units to use throughout the simulation
+    sim->G = 4*M_PI*M_PI;  //  Gravitational constant in AU^3/M_sun/yr^2
+
+    sim->heartbeat = heartbeat;
+
+    // add the host star
+    struct reb_particle star = {0}; 
+    star.m     = 1.;   
+    reb_simulation_add(sim, star); 
+
+    double m = 3.e-6;                // roughly 1 Earth Mass in Solar Masses
+    double a = 0.1;
+    double e = 0.05;
+    double inc = 5.*M_PI/180.;        // 5 degrees in radians
+    double Omega = 0.;
+    double omega = 0.;
+    double f = 0.;
+
+    struct reb_particle planet = reb_particle_from_orbit(sim->G, star, m, a, e, inc, Omega, omega, f);
+    reb_simulation_add(sim,planet);
+
+    // move simulation to center-of-mass frame
+    reb_simulation_move_to_com(sim);
+
+    // add in reboundx and load/add new force
+    struct rebx_extras* rebx = rebx_attach(sim);
+    struct rebx_force* gd = rebx_load_force(rebx, "gas_damping_timescale");
+    rebx_add_force(rebx, gd);
+
+    // Set the total simulation time  
+    double tmax = 5.e2;      // 500 years
+
+    // Set parameters for simulation
+    rebx_set_param_double(rebx, &gd->ap, "cs_coeff", 0.272);         // set speed sound coefficient to 0.272 AU^(3/4) yr^-1
+    rebx_set_param_double(rebx, &gd->ap, "tau_coeff", 0.003);        // set timescale coefficient to 0.003 yr AU^-2
+
+    // Set parameter for particle
+    rebx_set_param_double(rebx, &sim->particles[1].ap, "d_factor", 5.);         // set depletion factor to 5
+
+    // integrate simulation
+    reb_simulation_integrate(sim, tmax);
+    rebx_free(rebx);            // free all the memory allocated by reboundx
+    reb_simulation_free(sim);   // free all the memory allocated by rebound 
+}
+
+void heartbeat(struct reb_simulation* sim){
+    // output a e i of the planet
+    if(reb_simulation_output_check(sim, 50.)){
+        const struct reb_particle star = sim->particles[0];
+        const struct reb_orbit orbit = reb_orbit_from_particle(sim->G, sim->particles[1], star); // calculate orbit of planet
+        printf("%f\t%f\t%f\t%f\n",sim->t, orbit.a, orbit.e, orbit.inc);
+    }
+}

setup.py

@@ -57,7 +57,7 @@ def finalize_options(self):
         print("***", rebdir, "***", sitepackagesdir, "***", editable_rebdir, "***")
         self.include_dirs.append(rebdir)
         #self.include_dirs.append(editable_rebdir)
-        sources = [ 'src/modify_mass.c', 'src/integrator_euler.c', 'src/modify_orbits_forces.c', 'src/lense_thirring.c', 'src/integrator_rk2.c', 'src/track_min_distance.c', 'src/tides_spin.c', 'src/gas_dynamical_friction.c', 'src/rebxtools.c', 'src/inner_disk_edge.c', 'src/gravitational_harmonics.c', 'src/gr_potential.c', 'src/core.c', 'src/integrator_rk4.c', 'src/input.c', 'src/central_force.c', 'src/stochastic_forces.c', 'src/gr.c', 'src/modify_orbits_direct.c', 'src/tides_constant_time_lag.c', 'src/yarkovsky_effect.c', 'src/gr_full.c', 'src/steppers.c', 'src/integrate_force.c', 'src/interpolation.c', 'src/type_I_migration.c', 'src/output.c', 'src/radiation_forces.c', 'src/integrator_implicit_midpoint.c', 'src/exponential_migration.c', 'src/linkedlist.c'],
+        sources = [ 'src/central_force.c', 'src/core.c', 'src/exponential_migration.c', 'src/gas_damping_timescale.c', 'src/gas_dynamical_friction.c', 'src/gr.c', 'src/gr_full.c', 'src/gr_potential.c', 'src/gravitational_harmonics.c', 'src/inner_disk_edge.c', 'src/input.c', 'src/integrate_force.c', 'src/integrator_euler.c', 'src/integrator_implicit_midpoint.c', 'src/integrator_rk2.c', 'src/integrator_rk4.c', 'src/interpolation.c', 'src/lense_thirring.c', 'src/linkedlist.c', 'src/modify_mass.c', 'src/modify_orbits_direct.c', 'src/modify_orbits_forces.c', 'src/output.c', 'src/radiation_forces.c', 'src/rebxtools.c', 'src/steppers.c', 'src/stochastic_forces.c', 'src/tides_constant_time_lag.c', 'src/tides_spin.c', 'src/track_min_distance.c', 'src/type_I_migration.c', 'src/yarkovsky_effect.c'],
 
         self.library_dirs.append(rebdir+'/../')
         self.library_dirs.append(sitepackagesdir)
@@ -90,7 +90,7 @@ def finalize_options(self):
     extra_compile_args.append('-ffp-contract=off')
 
 libreboundxmodule = Extension('libreboundx',
-        sources = [ 'src/modify_mass.c', 'src/integrator_euler.c', 'src/modify_orbits_forces.c', 'src/lense_thirring.c', 'src/integrator_rk2.c', 'src/track_min_distance.c', 'src/tides_spin.c', 'src/gas_dynamical_friction.c', 'src/rebxtools.c', 'src/inner_disk_edge.c', 'src/gravitational_harmonics.c', 'src/gr_potential.c', 'src/core.c', 'src/integrator_rk4.c', 'src/input.c', 'src/central_force.c', 'src/stochastic_forces.c', 'src/gr.c', 'src/modify_orbits_direct.c', 'src/tides_constant_time_lag.c', 'src/yarkovsky_effect.c', 'src/gr_full.c', 'src/steppers.c', 'src/integrate_force.c', 'src/interpolation.c', 'src/type_I_migration.c', 'src/output.c', 'src/radiation_forces.c', 'src/integrator_implicit_midpoint.c', 'src/exponential_migration.c', 'src/linkedlist.c'],
+        sources = [ 'src/central_force.c', 'src/core.c', 'src/exponential_migration.c', 'src/gas_damping_timescale.c', 'src/gas_dynamical_friction.c', 'src/gr.c', 'src/gr_full.c', 'src/gr_potential.c', 'src/gravitational_harmonics.c', 'src/inner_disk_edge.c', 'src/input.c', 'src/integrate_force.c', 'src/integrator_euler.c', 'src/integrator_implicit_midpoint.c', 'src/integrator_rk2.c', 'src/integrator_rk4.c', 'src/interpolation.c', 'src/lense_thirring.c', 'src/linkedlist.c', 'src/modify_mass.c', 'src/modify_orbits_direct.c', 'src/modify_orbits_forces.c', 'src/output.c', 'src/radiation_forces.c', 'src/rebxtools.c', 'src/steppers.c', 'src/stochastic_forces.c', 'src/tides_constant_time_lag.c', 'src/tides_spin.c', 'src/track_min_distance.c', 'src/type_I_migration.c', 'src/yarkovsky_effect.c'],
                     include_dirs = ['src'],
                     library_dirs = [],
                     runtime_library_dirs = ["."],

src/Makefile

@@ -22,7 +22,7 @@ ifndef REBXGITHASH
 	PREDEF+= -DREBXGITHASH=$(REBXGITHASH)
 endif
 
-SOURCES=modify_mass.c integrator_euler.c modify_orbits_forces.c lense_thirring.c integrator_rk2.c track_min_distance.c tides_spin.c gas_dynamical_friction.c rebxtools.c inner_disk_edge.c gravitational_harmonics.c gr_potential.c core.c integrator_rk4.c input.c central_force.c stochastic_forces.c gr.c modify_orbits_direct.c tides_constant_time_lag.c yarkovsky_effect.c gr_full.c steppers.c integrate_force.c interpolation.c type_I_migration.c output.c radiation_forces.c integrator_implicit_midpoint.c exponential_migration.c linkedlist.c 
+SOURCES=central_force.c core.c exponential_migration.c gas_damping_timescale.c gas_dynamical_friction.c gr.c gr_full.c gr_potential.c gravitational_harmonics.c inner_disk_edge.c input.c integrate_force.c integrator_euler.c integrator_implicit_midpoint.c integrator_rk2.c integrator_rk4.c interpolation.c lense_thirring.c linkedlist.c modify_mass.c modify_orbits_direct.c modify_orbits_forces.c output.c radiation_forces.c rebxtools.c steppers.c stochastic_forces.c tides_constant_time_lag.c tides_spin.c track_min_distance.c type_I_migration.c yarkovsky_effect.c 
 
 OBJECTS=$(SOURCES:.c=.o)
 HEADERS=rebxtools.h reboundx.h linkedlist.h

src/core.c

@@ -62,6 +62,9 @@ void rebx_register_default_params(struct rebx_extras* rebx){
     rebx_register_param(rebx, "R_eq", REBX_TYPE_DOUBLE);
     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, "tau_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);
@@ -293,6 +296,10 @@ struct rebx_force* rebx_load_force(struct rebx_extras* const rebx, const char* n
         force->update_accelerations = rebx_modify_orbits_forces;
         force->force_type = REBX_FORCE_VEL;
     }
+    else if (strcmp(name, "gas_damping_timescale") == 0){
+        force->update_accelerations = rebx_gas_damping_timescale;
+        force->force_type = REBX_FORCE_VEL;
+    }
     else if (strcmp(name, "exponential_migration") == 0){
         force->update_accelerations = rebx_exponential_migration;
         force->force_type = REBX_FORCE_VEL;

src/core.h

@@ -70,6 +70,7 @@ void rebx_gr_potential(struct reb_simulation* const sim, struct rebx_force* cons
 void rebx_radiation_forces(struct reb_simulation* const sim, struct rebx_force* const force, struct reb_particle* const particles, const int N);
 void rebx_stochastic_forces(struct reb_simulation* const sim, struct rebx_force* const force, struct reb_particle* const particles, const int N);
 void rebx_modify_orbits_forces(struct reb_simulation* const sim, struct rebx_force* const force, struct reb_particle* const particles, const int N);
+void rebx_gas_damping_timescale(struct reb_simulation* const sim, struct rebx_force* const force, struct reb_particle* const particles, const int N);
 void rebx_exponential_migration(struct reb_simulation* const sim, struct rebx_force* const force, struct reb_particle* const particles, const int N);
 void rebx_tides_constant_time_lag(struct reb_simulation* const sim, struct rebx_force* const force, struct reb_particle* const particles, const int N);
 void rebx_central_force(struct reb_simulation* const sim, struct rebx_force* const force, struct reb_particle* const particles, const int N);

src/gas_damping_timescale.c

@@ -0,0 +1,141 @@
+/**
+ * @file    gas_damping_timescale.c
+ * @brief   Update orbits with prescribed timescales by directly changing orbital elements after each timestep
+ * @author  Phoebe Sandhaus <[email protected]>
+ * 
+ * @section     LICENSE
+ * Copyright (c) 2015 Dan Tamayo, Hanno Rein
+ *
+ * This file is part of reboundx.
+ *
+ * reboundx is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * reboundx is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with rebound.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ * The section after the dollar signs gets built into the documentation by a script.  All lines must start with space * space like below.
+ * Tables always must be preceded and followed by a blank line.  See http://docutils.sourceforge.net/docs/user/rst/quickstart.html for a primer on rst.
+ * $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+ *
+ * $Gas Effects$       // Effect category (must be the first non-blank line after dollar signs and between dollar signs to be detected by script).
+ *
+ * ======================= ===============================================
+ * Authors                 Phoebe Sandhaus
+ * Implementation Paper    `Sandhaus et al. in prep`
+ * Based on                `Dawson et al. 2016 <https://ui.adsabs.harvard.edu/abs/2016ApJ...822...54D/abstract>; Kominami & Ida 2002 <https://ui.adsabs.harvard.edu/abs/2002Icar..157...43K/abstract>`_
+ * C Example               :ref:`c_example_gas_damping_timescale`
+ * Python Example          `GasDampingTimescale.ipynb <https://github.com/PhoebeSandhaus/reboundx_gas_damping/tree/main/ipython_examples/GasDampingTimescale.ipynb>`_
+ * ======================= ===============================================
+ * 
+ * This updates particles' positions and velocities between timesteps by first calculating a damping timescale for each individual particle, and then applying the timescale to damp both the eccentricity and inclination of the particle. Note: The timescale of damping should be much greater than a particle's orbital period. The damping force should also be small as compared to the gravitational forces on the particle.
+ * 
+ * **Effect Parameters**
+ * 
+ * ============================ =========== ==================================================================
+ * Field (C type)               Required    Description
+ * ============================ =========== ==================================================================
+ * None                         -           -
+ * ============================ =========== ==================================================================
+ *
+ * **Particle Parameters**
+ *
+ * ============================ =========== ==============================================================================
+ * Field (C type)               Required    Description
+ * ============================ =========== ==============================================================================
+ * d_factor (double)            Yes         Depletion factor d in Equation 16 from Dawson et al. 2016; d=1 corresponds roughly to the full minimum mass solar nebula with Sigma_gas (surface gas density) = 1700 g cm^-2 at 1 AU [for d=1]; d>1 corresponds to a more depleted nebula
+ * cs_coeff (double)            Yes         Sound speed coefficient; Changing the value will change assumed units; Example: If you are using the units: AU, M_sun, yr -->  cs_coeff = 0.272  # AU^(3/4) yr^-1
+ * tau_coeff (double)           Yes         Timescale coefficient;   Changing the value will change assumed units; Example: If you are using the units: AU, M_sun, yr --> tau_coeff = 0.003  #  yr AU^-2
+ * ============================ =========== ==============================================================================
+ * 
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "rebound.h"
+#include "reboundx.h"
+#include "rebxtools.h"
+
+static struct reb_vec3d rebx_calculate_gas_damping_timescale(struct reb_simulation* const sim, struct rebx_force* const force, struct reb_particle* planet, struct reb_particle* star){
+    struct rebx_extras* const rebx = sim->extras;
+    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 tau_coeff = rebx_get_param(rebx, force->ap, "tau_coeff");
+
+    struct reb_vec3d a = {0};
+
+    if (d_factor == NULL || cs_coeff == NULL || tau_coeff == NULL){
+        rebx_error(rebx, "Need to set d_factor, cs_coeff, tau_coeff parameters.  See examples in documentation.\n");
+        return a;
+    }
+
+    // initialize positions and velocities
+    const double dvx = planet->vx - star->vx;
+    const double dvy = planet->vy - star->vy;
+    const double dvz = planet->vz - star->vz;
+    const double dx = planet->x-star->x;
+    const double dy = planet->y-star->y;
+    const double dz = planet->z-star->z;
+    const double r2 = dx*dx + dy*dy + dz*dz;
+
+    // initial semimajor axis, eccentricity, and inclination
+    const double a0 = o.a;
+    const double e0 = o.e;
+    const double inc0 = o.inc;
+    const double starMass = star->m;
+    const double planetMass = planet->m;
+
+    // eccentricity and inclination timescales from Dawson+16 Eqn 16
+    double coeff;
+
+    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/vk) {
+            coeff = v_over_cs*v_over_cs*v_over_cs;
+        }
+        else {
+            coeff = v_over_cs*v_over_cs*v_over_cs*v_over_cs;
+        }
+    }
+
+    double tau_e = -*tau_coeff*(*d_factor)*a0*a0*(starMass/planetMass)*coeff;
+    double tau_inc = 2.*tau_e;  // from Kominami & Ida 2002 [Eqs. 2.9 and 2.10]
+
+
+    if (tau_e < INFINITY || tau_inc < INFINITY){
+        const double vdotr = dx*dvx + dy*dvy + dz*dvz;
+        const double prefac = 2*vdotr/r2/tau_e;
+        a.x += prefac*dx;
+        a.y += prefac*dy;
+        a.z += prefac*dz + 2.*dvz/tau_inc;
+    }
+    return a;
+}
+
+void rebx_gas_damping_timescale(struct reb_simulation* const sim, struct rebx_force* const force, struct reb_particle* const particles, const int N){
+    int* ptr = rebx_get_param(sim->extras, force->ap, "coordinates");
+    enum REBX_COORDINATES coordinates = REBX_COORDINATES_JACOBI; // Default
+    if (ptr != NULL){
+        coordinates = *ptr;
+    }
+    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);
+}