-
Notifications
You must be signed in to change notification settings - Fork 63
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
1 parent
8df4f28
commit b2c0d72
Showing
8 changed files
with
86 additions
and
86 deletions.
There are no files selected for viewing
File renamed without changes.
134 changes: 67 additions & 67 deletions
134
examples/gas_damping_forces/problem.c → examples/gas_damping_timescale/problem.c
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -1,67 +1,67 @@ | ||
/** | ||
* Gas Damping Forces Example | ||
* | ||
* This example shows how to add gas damping forces 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(); | ||
|
||
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_forces"); | ||
rebx_add_force(rebx, gd); | ||
|
||
// Set the total simulation time | ||
double tmax = 5.e2*2*M_PI; // 500 years in yr/2pi | ||
|
||
// Set parameter for particle | ||
rebx_set_param_double(rebx, &sim->particles[1].ap, "d_factor", 10.); // set depletion factor to 10 | ||
|
||
// 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.*2*M_PI)){ | ||
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); | ||
} | ||
} | ||
/** | ||
* Gas Damping Timescale Example | ||
* | ||
* This example shows how to add gas damping forces 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(); | ||
|
||
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*2*M_PI; // 500 years in yr/2pi | ||
|
||
// Set parameter for particle | ||
rebx_set_param_double(rebx, &sim->particles[1].ap, "d_factor", 10.); // set depletion factor to 10 | ||
|
||
// 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.*2*M_PI)){ | ||
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); | ||
} | ||
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -1,5 +1,5 @@ | ||
/** | ||
* @file gas_damping_forces.c | ||
* @file gas_damping_timescale.c | ||
* @brief Update orbits with prescribed timescales by directly changing orbital elements after each timestep | ||
* @author Phoebe Sandhaus <[email protected]> | ||
* | ||
|
@@ -31,8 +31,8 @@ | |
* 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_forces` | ||
* Python Example `GasDampingForces.ipynb <https://github.com/PhoebeSandhaus/reboundx_gas_damping/tree/main/ipython_examples/GasDampingForces.ipynb>`_ | ||
* 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. | ||
|
@@ -62,7 +62,7 @@ | |
#include "reboundx.h" | ||
#include "rebxtools.h" | ||
|
||
static struct reb_vec3d rebx_calculate_gas_damping_forces(struct reb_simulation* const sim, struct rebx_force* const force, struct reb_particle* planet, struct reb_particle* star){ | ||
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); | ||
|
||
|
@@ -134,13 +134,13 @@ static struct reb_vec3d rebx_calculate_gas_damping_forces(struct reb_simulation* | |
return a; | ||
} | ||
|
||
void rebx_gas_damping_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){ | ||
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_forces, particles, N); | ||
rebx_com_force(sim, force, coordinates, back_reactions_inclusive, reference_name, rebx_calculate_gas_damping_timescale, particles, N); | ||
} |