diff --git a/doc/RTD/source/ExternalPotentials/index.rst b/doc/RTD/source/ExternalPotentials/index.rst
index 0cdfefdeba1350c4cd17fca941540747c44d6bf8..0454da53ab5a591452bb68bc1860859ecd83dc26 100644
--- a/doc/RTD/source/ExternalPotentials/index.rst
+++ b/doc/RTD/source/ExternalPotentials/index.rst
@@ -150,22 +150,22 @@ The parameters of the model are:
6. Hernquist potential (``hernquist``)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-A potential that is given by the Hernquist (1990)
-potential:
+We can set up a potential as given by Hernquist (1990):
-:math:`\Phi(r) = - \frac{G_{\rm N}M}{r+a}.`
+:math:`\Phi(r) = - \frac{G_{\rm N}M}{r+a},`
-The free parameters of the Hernquist potential are mass, scale length,
-and softening. The potential can be set at any position in the box.
-The potential add an additional time step constrain that limits the
-time step to a maximum of a specified fraction of the circular orbital
+where :math:`M` is the total Hernquist mass and :math: `a` is the Hernquist-
+equivalent scale radius of the potential. The potential can be set at any
+position in the box. It adds an additional time step constraint that limits
+the time step to a maximum of a specified fraction of the circular orbital
time at the current radius of the particle. The other criteria
-(CFL, self-gravity, ...) are applied on top of this criterion.
+(CFL, self-gravity, ...) are applied on top of this criterion. For example, a
+fraction of 0.01 means that an orbital period would be resolved by 100 time steps.
-The Hernquist potential can be run in the most basic version, then only the
-Hernquist mass, scale length, softening length and fraction of the
-orbital time for the time step limit are used, the parameters of the
-model in this case are:
+In the most basic version, the Hernquist potential can be run by providing
+only the Hernquist mass, scale length, softening length and fraction of the
+orbital time for the time stepping. In this case the model parameters may
+look something like:
.. code:: YAML
@@ -173,20 +173,23 @@ model in this case are:
useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions
position: [0.,0.,0.] # Location of centre of isothermal potential with respect to centre of the box (if 0) otherwise absolute (if 1) (internal units)
mass: 1e12 # Mass of the Hernquist potential
- scalelength: 2.0 # scale length a
- timestep_mult: 0.01 # Dimensionless pre-factor for the time-step condition, basically determines the fraction of the orbital time we use to do the time integration
+ scalelength: 2.0 # scale length a (internal units)
+ timestep_mult: 0.01 # Dimensionless pre-factor for the time-step condition, determines the fraction of the orbital time we use to do the time integration; fraction of 0.01 means we resolve an orbit with 100 timesteps
epsilon: 0.2 # Softening size (internal units)
Besides the basic version, it is also possible to run the Hernquist
potential for idealised disk galaxies that follow the approach of
-Springel+ 2005. The default Hernquist potential uses a corrected value
-for the formulation that improves the match with the NFW (below) with
-the same M200 (Nobels+ in prep). Contrary to above, the idealised disk
-setup runs with a specified M200, concentration and reduced Hubble
-constant that set both the mass and scale length parameter. The reduced
-Hubble constant is only used to determine R200.
-
-The parameters of the model are:
+`Springel, Di Matteo & Hernquist (2005)
+<https://ui.adsabs.harvard.edu/abs/2005MNRAS.361..776S/abstract>`_. This
+potential, however, uses a corrected value of the formulation that improves
+the match with the NFW profile (below) with the same M200 (Nobels+ in prep).
+Contrary to the above (idealizeddisk: 0 setup), the idealised disk setup runs
+by specifying one out of :math:`M_{200}`, :math:`V_{200}`, or :math:`R_{200}`,
+plus concentration and reduced Hubble constant.
+
+In this case, we don't provide the 'mass' and 'scalelength' parameters, but
+'M200' (or 'V200', or 'R200') and 'concentration' :math:`c`, as well as reduced Hubble
+constant :math:`h` to define the potential. The parameters of the model may look something like:
.. code:: YAML
@@ -199,21 +202,52 @@ The parameters of the model are:
concentration: 9.0 # concentration of the Halo
diskfraction: 0.040 # Disk mass fraction
bulgefraction: 0.0 # Bulge mass fraction
- timestep_mult: 0.01 # Dimensionless pre-factor for the time-step condition, basically determines the fraction of the orbital time we use to do the time integration
+ timestep_mult: 0.01 # Dimensionless pre-factor for the time-step condition, determines the fraction of the orbital time we use to do the time integration; fraction of 0.01 means we resolve an orbit with 100 timesteps
epsilon: 0.2 # Softening size (internal units)
+The user should specify one out of 'M200', 'V200', or 'R200' to define
+the potential. The reduced Hubble constant is then used to determine the
+other two. Then, the scale radius is calculated as :math:`R_s = R_{200}/c`,
+where :math:`c` is the concentration, and the Hernquist-equivalent scale-length
+is calculated as:
+
+:math:`a = \frac{b+\sqrt{b}}{1-b} R_{200}`,
+
+where :math:`b = \frac{2}{c^2}(\ln(1+c) - \frac{c}{1+c})`.
+
+Two examples using the Hernquist potential can be found in ``swiftsim/examples/GravityTests/``.
+The ``Hernquist_radialinfall`` example puts 5 particles with zero velocity in a Hernquist
+potential, resulting in radial orbits. The ``Hernquist_circularorbit``example puts three
+particles on a circular orbit in a Hernquist potential, one at the inner region, one at the
+maximal circular velocity, and one in the outer region. To run these examples, SWIFT must
+be configured with the flag ``--with-ext-potential=hernquist``, or ``hernquist-sdmh05``
+(see below).
-7. Hernquist potential (``hernquist-sdmh05``)
-^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+7. Hernquist SDMH05 potential (``hernquist-sdmh05``)
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
This is the same potential as Hernquist with the difference being the
way that the idealised disk part is calculated. This potential uses
exactly the same approach as `Springel, Di Matteo & Hernquist (2005)
<https://ui.adsabs.harvard.edu/abs/2005MNRAS.361..776S/abstract>`_,
-this means that ICs generated with the original `MakeNewDisk` code can
-be used when using this potential. Contrary to the updated Hernquist
-potential (above) it is not possible to have an identically matched
-NFW potential.
+which means that ICs generated with the original `MakeNewDisk` code can
+be used with this potential. Contrary to the updated Hernquist
+potential (above), it is not possible to have an identically matched
+NFW potential in this case.
+
+The parameters needed for this potential are the same set of variables as
+above, i.e. 'mass' and 'scalelength' when we don't use the idealised
+disk, and 'concentration' plus one out of 'M200', 'V200', or 'R200' if
+we do. As one of the three is provided, the reduced Hubble constant is
+used to calculate the other two. Then, the scale radius is calculated
+using the NFW definition, :math:`R_s = R_{200}/c`, and the Hernquist-
+equivalent scale length is given by
+
+:math:`a = R_s \sqrt{2(\ln(1+c) - \frac{c}{1+c})}.`
+
+Runs that provide e.g. M200 and c (using idealised disk) are thus equivalent
+to providing mass and scale length if calculated by the above equation (without
+idealized disk).
8. Navarro-Frenk-White potential (``nfw``):
diff --git a/examples/GravityTests/Hernquist_circularorbit/plotprog.py b/examples/GravityTests/Hernquist_circularorbit/plotprog.py
index a19c66e7f30e0e4012a23a4d38dd23045deea6e2..881c74301ec4a466f6984ee081b7ac6e3d555f5b 100755
--- a/examples/GravityTests/Hernquist_circularorbit/plotprog.py
+++ b/examples/GravityTests/Hernquist_circularorbit/plotprog.py
@@ -22,7 +22,7 @@ import h5py
import matplotlib.pyplot as plt
from scipy.integrate import odeint
-t = np.linspace(0, 40, 1e5)
+t = np.linspace(0, 40, 100000)
y0 = [0, 10]
a = 30.0
G = 4.300927e-06
diff --git a/src/potential/hernquist/potential.h b/src/potential/hernquist/potential.h
index 2ccfbe109015b50e4a63d16b0ab6edce068aaa19..144181f432f1a754e2934964c6b7b7e4bfda9cf4 100644
--- a/src/potential/hernquist/potential.h
+++ b/src/potential/hernquist/potential.h
@@ -73,6 +73,9 @@ struct external_potential {
* time to get the time steps */
double timestep_mult;
+ /*! Inverse of the sqrt of G*M, a common factor */
+ double sqrtgm_inv;
+
/*! Mode to use 0 for simplest form of potential purely 1 for idealized
* galaxies */
int usedisk;
@@ -92,21 +95,18 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep(
const struct phys_const* restrict phys_const,
const struct gpart* restrict g) {
- const float G_newton = phys_const->const_newton_G;
-
/* Calculate the relative potential with respect to the centre of the
* potential */
const float dx = g->x[0] - potential->x[0];
const float dy = g->x[1] - potential->x[1];
const float dz = g->x[2] - potential->x[2];
- /* calculate the radius */
+ /* Calculate the radius */
const float r = sqrtf(dx * dx + dy * dy + dz * dz + potential->epsilon2);
- const float sqrtgm_inv = 1.f / sqrtf(G_newton * potential->mass);
/* Calculate the circular orbital period */
- const float period = 2.f * M_PI * sqrtf(r) * potential->al *
- (1 + r / potential->al) * sqrtgm_inv;
+ const float period =
+ 2.f * M_PI * sqrtf(r) * (potential->al + r) * potential->sqrtgm_inv;
/* Time-step as a fraction of the cirecular orbital time */
const float time_step = potential->timestep_mult * period;
@@ -211,7 +211,7 @@ static INLINE void potential_init_backend(
potential->x[2] += s->dim[2] / 2.;
}
- /* check whether we use the more advanced idealized disk setting */
+ /* Check whether we use the more advanced idealized disk setting */
potential->usedisk = parser_get_opt_param_int(
parameter_file, "HernquistPotential:idealizeddisk",
idealized_disk_default);
@@ -269,7 +269,7 @@ static INLINE void potential_init_backend(
potential->M200 = M200;
potential->R200 = R200;
- /* get the concentration from the parameter file */
+ /* Get the concentration from the parameter file */
potential->c = parser_get_param_double(parameter_file,
"HernquistPotential:concentration");
@@ -286,7 +286,7 @@ static INLINE void potential_init_backend(
const double b = 2. * cc_inv * cc_inv * (log(1. + cc) - cc / (1. + cc));
/* Calculate the Hernquist equivalent scale length */
- potential->al = (b + sqrt(b)) / (1 - b) * R200;
+ potential->al = R200 * (b + sqrt(b)) / (1 - b);
/* Define R200 inv*/
const double R200_inv = 1. / R200;
@@ -296,9 +296,9 @@ static INLINE void potential_init_backend(
(potential->R200 + potential->al) * R200_inv *
R200_inv * potential->M200;
- /* Depending on the disk mass and and the bulge mass the halo
- * gets a different mass, because of this we read the fractions
- * from the parameter file and calculate the absolute mass*/
+ /* Depending on the disk mass and and the bulge mass, the halo
+ * gets a different mass. Because of this, we read the fractions
+ * from the parameter file and calculate the absolute mass */
const double diskfraction = parser_get_param_double(
parameter_file, "HernquistPotential:diskfraction");
const double bulgefraction = parser_get_param_double(
@@ -318,12 +318,14 @@ static INLINE void potential_init_backend(
parser_get_param_double(parameter_file, "HernquistPotential:epsilon");
potential->epsilon2 = epsilon * epsilon;
- /* Compute the minimal time-step. */
- /* This is the circular orbital time at the softened radius */
+ /* Calculate a common factor in the calculation, i.e. 1/sqrt(GM)*/
const float sqrtgm = sqrtf(phys_const->const_newton_G * potential->mass);
- potential->mintime = 2.f * sqrtf(epsilon) * potential->al * M_PI *
- (1. + epsilon / potential->al) / sqrtgm *
- potential->timestep_mult;
+ potential->sqrtgm_inv = 1. / sqrtgm;
+
+ /* Compute the minimal time-step. */
+ /* This is a fraction of the circular orbital time at the softened radius */
+ potential->mintime = potential->timestep_mult * 2.f * sqrtf(epsilon) * M_PI *
+ (potential->al + epsilon) * potential->sqrtgm_inv;
}
/**
diff --git a/src/potential/hernquist_sdmh05/potential.h b/src/potential/hernquist_sdmh05/potential.h
index b73ee208d442842d48d51a4129aa8df6933cc74c..8755d24ef939c70c012a1f00fd693e6fc5b8668b 100644
--- a/src/potential/hernquist_sdmh05/potential.h
+++ b/src/potential/hernquist_sdmh05/potential.h
@@ -63,6 +63,9 @@ struct external_potential {
/*! Time-step condition pre-factor, is multiplied times the circular orbital
* time to get the time steps */
double timestep_mult;
+
+ /* Additional common parameter inverse of sqrt(GM)*/
+ double sqrtgm_inv;
};
/**
@@ -79,8 +82,6 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep(
const struct phys_const* restrict phys_const,
const struct gpart* restrict g) {
- const float G_newton = phys_const->const_newton_G;
-
/* Calculate the relative potential with respect to the centre of the
* potential */
const float dx = g->x[0] - potential->x[0];
@@ -89,11 +90,10 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep(
/* calculate the radius */
const float r = sqrtf(dx * dx + dy * dy + dz * dz + potential->epsilon2);
- const float sqrtgm_inv = 1.f / sqrtf(G_newton * potential->mass);
/* Calculate the circular orbital period */
- const float period = 2.f * M_PI * sqrtf(r) * potential->al *
- (1 + r / potential->al) * sqrtgm_inv;
+ const float period =
+ 2.f * M_PI * sqrtf(r) * (potential->al + r) * potential->sqrtgm_inv;
/* Time-step as a fraction of the cirecular orbital time */
const float time_step = potential->timestep_mult * period;
@@ -254,7 +254,7 @@ static INLINE void potential_init_backend(
/* message("M200 = %g, R200 = %g, V200 = %g", M200, R200, V200); */
/* message("H0 = %g", H0); */
- /* get the concentration from the parameter file */
+ /* Get the concentration from the parameter file */
const double concentration = parser_get_param_double(
parameter_file, "HernquistPotential:concentration");
@@ -262,12 +262,12 @@ static INLINE void potential_init_backend(
const double RS = R200 / concentration;
/* Calculate the Hernquist equivalent scale length */
- potential->al = RS * sqrt(1. * (log(1. + concentration) -
+ potential->al = RS * sqrt(2. * (log(1. + concentration) -
concentration / (1. + concentration)));
- /* Depending on the disk mass and and the bulge mass the halo
- * gets a different mass, because of this we read the fractions
- * from the parameter file and calculate the absolute mass*/
+ /* Depending on the disk mass and the bulge mass, the halo
+ * gets a different mass. Because of this, we read the fractions
+ * from the parameter file and calculate the absolute mass */
const double diskfraction = parser_get_param_double(
parameter_file, "HernquistPotential:diskfraction");
const double bulgefraction = parser_get_param_double(
@@ -287,12 +287,14 @@ static INLINE void potential_init_backend(
parser_get_param_double(parameter_file, "HernquistPotential:epsilon");
potential->epsilon2 = epsilon * epsilon;
- /* Compute the minimal time-step. */
- /* This is the circular orbital time at the softened radius */
+ /* Calculate a common factor in the calculation, i.e. 1/sqrt(GM)*/
const float sqrtgm = sqrtf(phys_const->const_newton_G * potential->mass);
- potential->mintime = 2.f * sqrtf(epsilon) * potential->al * M_PI *
- (1. + epsilon / potential->al) / sqrtgm *
- potential->timestep_mult;
+ potential->sqrtgm_inv = 1. / sqrtgm;
+
+ /* Compute the minimal time-step. */
+ /* This is a fraction of the circular orbital time at the softened radius */
+ potential->mintime = potential->timestep_mult * 2.f * sqrtf(epsilon) * M_PI *
+ (potential->al + epsilon) * potential->sqrtgm_inv;
}
/**
@@ -304,7 +306,7 @@ static inline void potential_print_backend(
const struct external_potential* potential) {
message(
- "external potential is 'hernquist Springel, Di Matteo & Herquist 2005' "
+ "external potential is 'hernquist Springel, Di Matteo & Hernquist 2005' "
"with properties are (x,y,z) = (%e, %e, %e), mass = %e scale length = %e "
", minimum time = %e timestep multiplier = %e",
potential->x[0], potential->x[1], potential->x[2], potential->mass,