Commit 6102664c authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Applied code formatting rule, fixed compilation bug, removed non-free legend.pro files.

parent b2345b3c
sets-up for a potential of a patch f disk, see Creasey, Theuns &
Bower,2013.
Setup for a potential of a patch disk, see Creasey, Theuns &
Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948
rho(z) = (Sigma/2b) / cosh^2(z/b)
where Sigma is the surface density, and b the scale height
......
This diff is collapsed.
......@@ -3,8 +3,8 @@
# Generate the initial conditions if they are not present.
if [ ! -e Isothermal.hdf5 ]
then
echo "Generating initial conditions for the point mass potential box example..."
echo "Generating initial conditions for the disk-patch example..."
python makeIC.py 1000
fi
../swift -g -t 2 externalGravity.yml
../../swift -g -t 2 disk-patch.yml
genertae and evolve a disk-patch, where gas is in hydrostatic
generate and evolve a disk-patch, where gas is in hydrostatic
equilibrium with an imposed external gravutation force, using teh
equations from Creasey, Theuns & Bower 2013.
equations from Creasey, Theuns & Bower, 2013, MNRAS, Volume 429,
Issue 3, p.1922-1948
Run the swift using setting
#define EXTERNAL_POTENTIAL_DISK_PATCH
......
This diff is collapsed.
This diff is collapsed.
......@@ -329,11 +329,12 @@ int main(int argc, char *argv[]) {
/* Initialise the hydro properties */
struct hydro_props hydro_properties;
hydro_props_init(&hydro_properties, params);
if (with_hydro && myrank == 0) eos_print();
if (with_hydro && myrank == 0) eos_print();
/* Initialise the external potential properties */
struct external_potential potential;
if (with_external_gravity) potential_init(params, &prog_const, &us, &potential);
if (with_external_gravity)
potential_init(params, &prog_const, &us, &potential);
if (with_external_gravity && myrank == 0) potential_print(&potential);
/* Read particles and space information from (GADGET) ICs */
......
......@@ -48,7 +48,8 @@
/* Equation of state choice */
#define EOS_IDEAL_GAS
/* if EOS_ISOTHERMAL_GAS is defined, keep thermal energy per unit mass equal to EOS_ISOTHERMAL_GAS in programme units */
/* if EOS_ISOTHERMAL_GAS is defined, keep thermal energy per unit mass equal to
* EOS_ISOTHERMAL_GAS in programme units */
//#define EOS_ISOTHERMAL_GAS (20.2615290634)
/* Kernel function to use */
......@@ -74,10 +75,10 @@
//#define EXTERNAL_POTENTIAL_POINTMASS
//#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
#define EXTERNAL_POTENTIAL_DISK_PATCH
/* Add viscuous force to gas particles to speed-up glass making for disk-patch ICs */
/* Add viscuous force to gas particles to speed-up glass making for disk-patch
* ICs */
//#define EXTERNAL_POTENTIAL_DISK_PATCH_ICS
/* Are we debugging ? */
//#define SWIFT_DEBUG_CHECKS
......
......@@ -41,7 +41,7 @@
/* ------------------------------------------------------------------------- */
#if defined(EOS_IDEAL_GAS)
#if defined(EOS_ISOTHERMAL_GAS)
#error: cannot have ideal and isothermal gas at the same time!
#error : cannot have ideal and isothermal gas at the same time!
#endif
/**
......@@ -128,15 +128,16 @@ gas_soundspeed_from_internal_energy(float density, float u) {
return sqrtf(u * hydro_gamma * hydro_gamma_minus_one);
}
__attribute__((always_inline)) INLINE static void eos_print(){
message(" assuming an equation of state of an ideal gas, with gamma = %e", hydro_gamma);
__attribute__((always_inline)) INLINE static void eos_print() {
message(" assuming an equation of state of an ideal gas, with gamma = %e",
hydro_gamma);
}
__attribute__((always_inline)) INLINE static float hydro_update_entropy(float density, float entropy)
{
__attribute__((always_inline)) INLINE static float hydro_update_entropy(
float density, float entropy) {
return entropy;
}
__attribute__((always_inline)) INLINE static float hydro_update_entropy_rate(float density, float entropy)
{
__attribute__((always_inline)) INLINE static float hydro_update_entropy_rate(
float density, float entropy) {
return hydro_gamma_minus_one * pow_minus_gamma_minus_one(density);
}
/* ------------------------------------------------------------------------- */
......@@ -174,7 +175,8 @@ __attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
__attribute__((always_inline)) INLINE static float
gas_entropy_from_internal_energy(float density, float u) {
return hydro_gamma_minus_one * EOS_ISOTHERMAL_GAS / pow(density, hydro_gamma-1);
return hydro_gamma_minus_one * EOS_ISOTHERMAL_GAS /
pow(density, hydro_gamma - 1);
}
/**
......@@ -202,16 +204,19 @@ gas_soundspeed_from_internal_energy(float density, float u) {
error("Missing definition !");
return 0.f;
}
__attribute__((always_inline)) INLINE static void eos_print(){
message(" assuming an equation of state for an isothermal gas with utherm= %e and gamma = %e", EOS_ISOTHERMAL_GAS, hydro_gamma);
__attribute__((always_inline)) INLINE static void eos_print() {
message(
" assuming an equation of state for an isothermal gas with utherm= %e "
"and gamma = %e",
EOS_ISOTHERMAL_GAS, hydro_gamma);
}
__attribute__((always_inline)) INLINE static float hydro_update_entropy(float density, float entropy)
{
__attribute__((always_inline)) INLINE static float hydro_update_entropy(
float density, float entropy) {
float thermal_energy = EOS_ISOTHERMAL_GAS;
return gas_entropy_from_internal_energy(density, thermal_energy);
}
__attribute__((always_inline)) INLINE static float hydro_update_entropy_rate(float density, float entropy)
{
__attribute__((always_inline)) INLINE static float hydro_update_entropy_rate(
float density, float entropy) {
return 0;
}
/* ------------------------------------------------------------------------- */
......
......@@ -47,7 +47,8 @@ gravity_compute_timestep_external(const struct external_potential* potential,
phys_const, gp));
#endif
#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
dt = fmin(dt, external_gravity_disk_patch_timestep(potential, phys_const, gp));
dt =
fmin(dt, external_gravity_disk_patch_timestep(potential, phys_const, gp));
#endif
return dt;
}
......@@ -74,7 +75,6 @@ gravity_compute_timestep_self(const struct phys_const* const phys_const,
return FLT_MAX;
}
/**
* @brief Initialises the g-particles for the first time
*
......@@ -117,7 +117,7 @@ __attribute__((always_inline)) INLINE static void gravity_init_gpart(
* @param const_G Newton's constant
*/
__attribute__((always_inline)) INLINE static void gravity_end_force(
struct gpart* gp, const double const_G) {
struct gpart* gp, const double const_G) {
/* Let's get physical... */
gp->a_grav[0] *= const_G;
......@@ -135,8 +135,7 @@ __attribute__((always_inline)) INLINE static void gravity_end_force(
* @param gp The particle to act upon.
*/
__attribute__((always_inline)) INLINE static void external_gravity(
const double time,
const struct external_potential* potential,
const double time, const struct external_potential* potential,
const struct phys_const* const phys_const, struct gpart* gp) {
#ifdef EXTERNAL_POTENTIAL_POINTMASS
......
......@@ -21,7 +21,6 @@
#include "hydro.h"
#include "io_properties.h"
#include "kernel_hydro.h"
#include "equation_of_state.h"
/**
* @brief Specifies which particle fields to read from a dataset
......
......@@ -64,7 +64,7 @@ struct part {
/* Derivative of the density with respect to smoothing length. */
float rho_dh;
// union {
union {
struct {
......@@ -100,6 +100,7 @@ struct part {
float h_dt;
} force;
};
/* Particle ID. */
long long id;
......
......@@ -33,7 +33,7 @@
* @param potential The external potential properties to initialize
*/
void potential_init(const struct swift_params* parameter_file,
const struct phys_const* const phys_const,
const struct phys_const* const phys_const,
struct UnitSystem* us,
struct external_potential* potential) {
......@@ -67,17 +67,19 @@ void potential_init(const struct swift_params* parameter_file,
#endif /* EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL */
#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
potential -> disk_patch_potential.surface_density =
parser_get_param_double(parameter_file, "Disk-PatchPotential:surface_density");
potential -> disk_patch_potential.scale_height =
parser_get_param_double(parameter_file, "Disk-PatchPotential:scale_height");
potential -> disk_patch_potential.z_disk =
parser_get_param_double(parameter_file, "Disk-PatchPotential:z_disk");
potential -> disk_patch_potential.timestep_mult =
parser_get_param_double(parameter_file, "Disk-PatchPotential:timestep_mult");
potential -> disk_patch_potential.dynamical_time = sqrt(potential->disk_patch_potential.scale_height / (phys_const->const_newton_G * potential->disk_patch_potential.surface_density));
potential->disk_patch_potential.surface_density = parser_get_param_double(
parameter_file, "Disk-PatchPotential:surface_density");
potential->disk_patch_potential.scale_height = parser_get_param_double(
parameter_file, "Disk-PatchPotential:scale_height");
potential->disk_patch_potential.z_disk =
parser_get_param_double(parameter_file, "Disk-PatchPotential:z_disk");
potential->disk_patch_potential.timestep_mult = parser_get_param_double(
parameter_file, "Disk-PatchPotential:timestep_mult");
potential->disk_patch_potential.dynamical_time =
sqrt(potential->disk_patch_potential.scale_height /
(phys_const->const_newton_G *
potential->disk_patch_potential.surface_density));
#endif /* EXTERNAL_POTENTIAL_DISK_PATCH */
}
/**
......@@ -108,13 +110,17 @@ void potential_print(const struct external_potential* potential) {
#endif /* EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL */
#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
message("Disk-patch potential properties are surface_density = %e disk height= %e scale height= %e timestep multiplier= %e",
potential->disk_patch_potential.surface_density,
potential->disk_patch_potential.z_disk,
potential->disk_patch_potential.scale_height,
potential->disk_patch_potential.timestep_mult);
message(
"Disk-patch potential properties are surface_density = %e disk height= "
"%e scale height= %e timestep multiplier= %e",
potential->disk_patch_potential.surface_density,
potential->disk_patch_potential.z_disk,
potential->disk_patch_potential.scale_height,
potential->disk_patch_potential.timestep_mult);
#ifdef EXTERNAL_POTENTIAL_DISK_PATCH_ICS
message("Disk-patch potential: imposing growth of gravity over time, and adding viscous force to gravity");
message(
"Disk-patch potential: imposing growth of gravity over time, and adding "
"viscous force to gravity");
#endif
#endif /* EXTERNAL_POTENTIAL_DISK_PATCH */
}
......@@ -25,8 +25,8 @@
#include "../config.h"
/* Some standard headers. */
#include <math.h>
#include <float.h>
#include <math.h>
/* Local includes. */
#include "const.h"
......@@ -35,7 +35,7 @@
#include "part.h"
#include "physical_constants.h"
#include "units.h"
#include "math.h"
/* External Potential Properties */
struct external_potential {
......@@ -56,11 +56,11 @@ struct external_potential {
#endif
#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
struct {
double surface_density;
double scale_height;
double z_disk;
double dynamical_time;
double timestep_mult;
double surface_density;
double scale_height;
double z_disk;
double dynamical_time;
double timestep_mult;
} disk_patch_potential;
#endif
};
......@@ -68,17 +68,17 @@ struct external_potential {
/* external potential: disk_patch */
#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
/**
* @brief Computes the time-step due to the acceleration from a hydrostatic disk (Creasy, Theuns & Bower, 2013)
* @brief Computes the time-step due to the acceleration from a hydrostatic disk
* (Creasy, Theuns & Bower, 2013)
*
* @param phys_cont The physical constants in internal units.
* @param gp Pointer to the g-particle data.
*/
__attribute__((always_inline)) INLINE static float external_gravity_disk_patch_timestep(
const struct external_potential* potential,
const struct phys_const* const phys_const,
const struct gpart* const g) {
__attribute__((always_inline)) INLINE static float
external_gravity_disk_patch_timestep(const struct external_potential* potential,
const struct phys_const* const phys_const,
const struct gpart* const g) {
/* initilize time step to disk dynamical time */
const float dt_dyn = potential->disk_patch_potential.dynamical_time;
float dt = dt_dyn;
......@@ -87,78 +87,85 @@ __attribute__((always_inline)) INLINE static float external_gravity_disk_patch_t
const float dz = fabs(g->x[2] - potential->disk_patch_potential.z_disk);
/* vertical cceleration */
const float z_accel = 2 * M_PI * phys_const->const_newton_G * potential->disk_patch_potential.surface_density
* tanh(dz / potential->disk_patch_potential.scale_height);
const float z_accel = 2 * M_PI * phys_const->const_newton_G *
potential->disk_patch_potential.surface_density *
tanh(dz / potential->disk_patch_potential.scale_height);
/* demand that dt * velocity < fraction of scale height of disk */
float dt1 = FLT_MAX;
if(fabs(g->v_full[2]) > 0)
{
dt1 = potential->disk_patch_potential.scale_height / fabs(g->v_full[2]);
if(dt1 < dt) dt = dt1;
}
if (fabs(g->v_full[2]) > 0) {
dt1 = potential->disk_patch_potential.scale_height / fabs(g->v_full[2]);
if (dt1 < dt) dt = dt1;
}
/* demand that dt^2 * acceleration < fraction of scale height of disk */
float dt2 = FLT_MAX;
if(fabs(z_accel) > 0)
{
dt2 = potential->disk_patch_potential.scale_height / fabs(z_accel);
if(dt2 < dt * dt) dt = sqrt(dt2);
}
if (fabs(z_accel) > 0) {
dt2 = potential->disk_patch_potential.scale_height / fabs(z_accel);
if (dt2 < dt * dt) dt = sqrt(dt2);
}
/* demand that dt^3 jerk < fraction of scale height of disk */
float dt3 = FLT_MAX;
if(abs(g->v_full[2]) > 0)
{
const float dz_accel_over_dt = 2 * M_PI * phys_const->const_newton_G * potential->disk_patch_potential.surface_density /
potential->disk_patch_potential.scale_height / cosh(dz / potential->disk_patch_potential.scale_height) / cosh(dz / potential->disk_patch_potential.scale_height) * fabs(g->v_full[2]);
dt3 = potential->disk_patch_potential.scale_height / fabs(dz_accel_over_dt);
if(dt3 < dt * dt * dt) dt = pow(dt3, 1./3.);
}
if (abs(g->v_full[2]) > 0) {
const float dz_accel_over_dt =
2 * M_PI * phys_const->const_newton_G *
potential->disk_patch_potential.surface_density /
potential->disk_patch_potential.scale_height /
cosh(dz / potential->disk_patch_potential.scale_height) /
cosh(dz / potential->disk_patch_potential.scale_height) *
fabs(g->v_full[2]);
dt3 = potential->disk_patch_potential.scale_height / fabs(dz_accel_over_dt);
if (dt3 < dt * dt * dt) dt = pow(dt3, 1. / 3.);
}
return potential->disk_patch_potential.timestep_mult * dt;
}
/**
* @brief Computes the gravitational acceleration from a hydrostatic disk (Creasy, Theuns & Bower, 2013)
* @brief Computes the gravitational acceleration from a hydrostatic disk
* (Creasy, Theuns & Bower, 2013)
*
* @param phys_cont The physical constants in internal units.
* @param gp Pointer to the g-particle data.
*/
__attribute__((always_inline)) INLINE static void external_gravity_disk_patch_potential(const double time, const struct external_potential* potential, const struct phys_const* const phys_const, struct gpart* g) {
__attribute__((always_inline)) INLINE static void
external_gravity_disk_patch_potential(
const double time, const struct external_potential* potential,
const struct phys_const* const phys_const, struct gpart* g) {
const float G_newton = phys_const->const_newton_G;
const float t_dyn = potential->disk_patch_potential.dynamical_time;
const float dz = g->x[2] - potential->disk_patch_potential.z_disk;
g->a_grav[0] += 0;
g->a_grav[1] += 0;
g->a_grav[0] += 0;
g->a_grav[1] += 0;
#ifdef EXTERNAL_POTENTIAL_DISK_PATCH_ICS
float reduction_factor = 1.;
if(time < 5 * t_dyn)
reduction_factor = time / (5 * t_dyn);
if (time < 5 * t_dyn) reduction_factor = time / (5 * t_dyn);
#else
const float reduction_factor = 1.;
#endif
const float z_accel = reduction_factor * 2 * G_newton * M_PI * potential->disk_patch_potential.surface_density
* tanh(fabs(dz) / potential->disk_patch_potential.scale_height);
const float z_accel =
reduction_factor * 2 * G_newton * M_PI *
potential->disk_patch_potential.surface_density *
tanh(fabs(dz) / potential->disk_patch_potential.scale_height);
if (dz > 0)
g->a_grav[2] -= z_accel / G_newton; /* returned acceleraton is multiplied by G later on */
if (dz < 0)
g->a_grav[2] += z_accel / G_newton;
g->a_grav[2] -=
z_accel /
G_newton; /* returned acceleraton is multiplied by G later on */
if (dz < 0) g->a_grav[2] += z_accel / G_newton;
if(abs(g->id_or_neg_offset) == 1)
message(" time= %e, rf= %e, az= %e", time, reduction_factor, g->a_grav[2]);
if (abs(g->id_or_neg_offset) == 1)
message(" time= %e, rf= %e, az= %e", time, reduction_factor, g->a_grav[2]);
#ifdef EXTERNAL_POTENTIAL_DISK_PATCH_ICS
/* TT: add viscous drag */
g->a_grav[0] -= g->v_full[0] / (2 * t_dyn) / G_newton;
g->a_grav[1] -= g->v_full[1] / (2 * t_dyn) / G_newton;
g->a_grav[2] -= g->v_full[2] / (2 * t_dyn) / G_newton;
#endif
g->a_grav[0] -= g->v_full[0] / (2 * t_dyn) / G_newton;
g->a_grav[1] -= g->v_full[1] / (2 * t_dyn) / G_newton;
g->a_grav[2] -= g->v_full[2] / (2 * t_dyn) / G_newton;
#endif
}
#endif /* EXTERNAL_POTENTIAL_DISK_PATCH */
......@@ -177,129 +184,145 @@ external_gravity_isothermalpotential_timestep(
const struct phys_const* const phys_const, const struct gpart* const g) {
const struct gpart* const g) {
const float dx = g->x[0] - potential->isothermal_potential.x;
const float dy = g->x[1] - potential->isothermal_potential.y;
const float dz = g->x[2] - potential->isothermal_potential.z;
const float rinv2 = 1.f / (dx * dx + dy * dy + dz * dz);
const float drdv =
dx * (g->v_full[0]) + dy * (g->v_full[1]) + dz * (g->v_full[2]);
const double vrot = potential->isothermal_potential.vrot;
const float dota_x =
vrot * vrot * rinv2 * (g->v_full[0] - 2 * drdv * dx * rinv2);
const float dota_y =
vrot * vrot * rinv2 * (g->v_full[1] - 2 * drdv * dy * rinv2);
const float dota_z =
vrot * vrot * rinv2 * (g->v_full[2] - 2 * drdv * dz * rinv2);
const float dota_2 = dota_x * dota_x + dota_y * dota_y + dota_z * dota_z;
const float a_2 = g->a_grav[0] * g->a_grav[0] + g->a_grav[1] * g->a_grav[1] +
g->a_grav[2] * g->a_grav[2];
return potential->isothermal_potential.timestep_mult * sqrtf(a_2 / dota_2);
}
/**
* @brief Computes the gravitational acceleration of a particle due to a point
* mass
*
* Note that the accelerations are multiplied by Newton's G constant later on.
*
* @param potential The #external_potential used in the run.
* @param phys_const The physical constants in internal units.
* @param g Pointer to the g-particle data.
*/
__attribute__((always_inline)) INLINE static void external_gravity_isothermalpotential(const struct external_potential* potential, const struct phys_const* const phys_const, struct gpart* g) {
__attribute__((always_inline)) INLINE static void
external_gravity_isothermalpotential(const struct external_potential* potential,
const struct phys_const* const phys_const,
struct gpart* g) {
const float G_newton = phys_const->const_newton_G;
const float dx = g->x[0] - potential->isothermal_potential.x;
const float dy = g->x[1] - potential->isothermal_potential.y;
const float dz = g->x[2] - potential->isothermal_potential.z;
const float rinv2 = 1.f / (dx * dx + dy * dy + dz * dz);
const double vrot = potential->isothermal_potential.vrot;
const double term = -vrot * vrot * rinv2 / G_newton;
g->a_grav[0] += term * dx;
g->a_grav[1] += term * dy;
g->a_grav[2] += term * dz;
}
const float dx = g->x[0] - potential->isothermal_potential.x;
const float dy = g->x[1] - potential->isothermal_potential.y;
const float dz = g->x[2] - potential->isothermal_potential.z;
const float rinv2 = 1.f / (dx * dx + dy * dy + dz * dz);
const float drdv =
dx * (g->v_full[0]) + dy * (g->v_full[1]) + dz * (g->v_full[2]);
const double vrot = potential->isothermal_potential.vrot;
const float dota_x =
vrot * vrot * rinv2 * (g->v_full[0] - 2 * drdv * dx * rinv2);
const float dota_y =
vrot * vrot * rinv2 * (g->v_full[1] - 2 * drdv * dy * rinv2);
const float dota_z =
vrot * vrot * rinv2 * (g->v_full[2] - 2 * drdv * dz * rinv2);
const float dota_2 =
dota_x * dota_x + dota_y * dota_y + dota_z * dota_z;
const float a_2 = g->a_grav[0] * g->a_grav[0] +
g->a_grav[1] * g->a_grav[1] +
g->a_grav[2] * g->a_grav[2];
return potential->isothermal_potential.timestep_mult *
sqrtf(a_2 / dota_2);
}
/**
* @brief Computes the gravitational acceleration of a particle due to a
* point
* mass
*
* Note that the accelerations are multiplied by Newton's G constant
* later on.
*
* @param potential The #external_potential used in the run.
* @param phys_const The physical constants in internal units.
* @param g Pointer to the g-particle data.
*/
__attribute__((always_inline)) INLINE static void
external_gravity_isothermalpotential(
const struct external_potential* potential,
const struct phys_const* const phys_const, struct gpart* g) {
__attribute__((always_inline)) INLINE static void
external_gravity_isothermalpotential(
const struct external_potential* potential,
const struct phys_const* const phys_const, struct gpart* g) {
const float G_newton = phys_const->const_newton_G;
const float dx = g->x[0] - potential->isothermal_potential.x;
const float dy = g->x[1] - potential->isothermal_potential.y;
const float dz = g->x[2] - potential->isothermal_potential.z;
const float rinv2 = 1.f / (dx * dx + dy * dy + dz * dz);
const double vrot = potential->isothermal_potential.vrot;
const double term = -vrot * vrot * rinv2 / G_newton;
g->a_grav[0] += term * dx;
g->a_grav[1] += term * dy;
g->a_grav[2] += term * dz;
}
#endif /* EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL */
/* Include exteral pointmass potential */
#ifdef EXTERNAL_POTENTIAL_POINTMASS
/**
* @brief Computes the time-step due to the acceleration from a point mass
*
* @param potential The properties of the externa potential.
* @param phys_const The physical constants in internal units.
* @param g Pointer to the g-particle data.
*/
__attribute__((always_inline)) INLINE static float external_gravity_pointmass_timestep(const struct external_potential* potential,
const struct phys_const* const phys_const,
const struct gpart* const g) {
const float G_newton = phys_const->const_newton_G;
const float dx = g->x[0] - potential->point_mass.x;
const float dy = g->x[1] - potential->point_mass.y;
const float dz = g->x[2] - potential->point_mass.z;
const float rinv = 1.f / sqrtf(dx * dx + dy * dy + dz * dz);
const float drdv = (g->x[0] - potential->point_mass.x) * (g->v_full[0]) +
(g->x[1] - potential->point_mass.y) * (g->v_full[1]) +
(g->x[2] - potential->point_mass.z) * (g->v_full[2]);
const float dota_x = G_newton * potential->point_mass.mass * rinv * rinv *
rinv * (-g->v_full[0] + 3.f * rinv * rinv * drdv * dx);
const float dota_y = G_newton * potential->point_mass.mass * rinv * rinv *
rinv * (-g->v_full[1] + 3.f * rinv * rinv * drdv * dy);
const float dota_z = G_newton * potential->point_mass.mass * rinv * rinv *
rinv * (-g->v_full[2] + 3.f * rinv * rinv * drdv * dz);
const float dota_2 = dota_x * dota_x + dota_y * dota_y + dota_z * dota_z;
const float a_2 = g->a_grav[0] * g->a_grav[0] + g->a_grav[1] * g->a_grav[1] +
g->a_grav[2] * g->a_grav[2];
return potential->point_mass.timestep_mult * sqrtf(a_2 / dota_2);