Add a new potential MWPotential2014

configure.ac

@@ -2624,7 +2624,7 @@ esac
 #  External potential
 AC_ARG_WITH([ext-potential],
    [AS_HELP_STRING([--with-ext-potential=<pot>],
-      [external potential @<:@none, point-mass, point-mass-softened, isothermal, nfw, nfw-mn, hernquist, hernquist-sdmh05, disc-patch, sine-wave, constant, default: none@:>@]
+      [external potential @<:@none, point-mass, point-mass-softened, isothermal, nfw, nfw-mn, hernquist, hernquist-sdmh05, disc-patch, sine-wave, MWPotential2014, constant, default: none@:>@]
    )],
    [with_potential="$withval"],
    [with_potential="none"]
@@ -2660,6 +2660,9 @@ case "$with_potential" in
    point-mass-softened)
       AC_DEFINE([EXTERNAL_POTENTIAL_POINTMASS_SOFT], [1], [Softened point-mass potential with form 1/(r^2 + softening^2).])
    ;;
+   MWPotential2014)
+      AC_DEFINE([EXTERNAL_POTENTIAL_MWPotential2014], [1], [Milky-Way like potential composed of a Navarro-Frenk-White + Miyamoto-Nagai disk + Power spherical cuttoff external potential.])
+   ;;
    constant)
       AC_DEFINE([EXTERNAL_POTENTIAL_CONSTANT], [1], [Constant gravitational acceleration.])
    ;;

doc/RTD/source/ExternalPotentials/index.rst

@@ -344,7 +344,46 @@ follows the definitions of `Creasey, Theuns & Bower (2013)
 <https://adsabs.harvard.edu/abs/2013MNRAS.429.1922C>`_ equations (16) and (17).
 The potential is implemented along the x-axis.
 
+12. MWPotential2014 (``MWPotential2014``)
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+This potential is based on ``galpy``'s ``MWPotential2014`` from `Jo Bovy (2015) <https://ui.adsabs.harvard.edu/abs/2015ApJS..216...29B>`_ and consists in a NFW potential for the halo, an axisymmetric Miyamoto-Nagai potential for the disk and a bulge modeled by a power spherical law with exponential cut-off. The bulge is given by the density:
+
+:math:`\rho(r) = A \left( \frac{r_1}{r} \right)^\alpha \exp \left( - \frac{r^2}{r_c^2} \right)`,
+
+where :math:`A` is an amplitude, :math:`r_1` is a reference radius for amplitude, :math:`\alpha` is the inner power and :math:`r_c` is the cut-off radius.
+
+The resulting potential is:
+
+:math:`\Phi_{\mathrm{MW}}(R, z) = f_1 \Phi_{\mathrm{NFW}} + f_2 \Phi_{\mathrm{MN}} + f_3 \Phi_{\text{bulge}}`,
 
+where :math:`R^2 = x^2 + y^2` is the projected radius and :math:`f_1`, :math:`f_2` and :math:`f_3` are three coefficients that adjust the strength of each individual component.
+
+The parameters of the model are:
+
+.. code:: YAML
+
+    MWPotential2014Potential:
+      useabspos:        0          # 0 -> positions based on centre, 1 -> absolute positions
+      position:         [0.,0.,0.] # Location of centre of potential with respect to centre of the box (if 0) otherwise absolute (if 1) (internal units)
+      timestep_mult:    0.005      # Dimensionless pre-factor for the time-step condition, basically determines the fraction of the orbital time we use to do the time integration
+      epsilon:          0.001      # Softening size (internal units)
+      concentration:    9.823403437774843      # concentration of the Halo
+      M_200_Msun:       147.41031542774076e10  # M200 of the galaxy disk (in M_sun)
+      H:                1.2778254614201471     # Hubble constant in units of km/s/Mpc
+      Mdisk_Msun:       6.8e10                 # Mass of the disk (in M_sun)
+      Rdisk_kpc:        3.0                    # Effective radius of the disk (in kpc)
+      Zdisk_kpc:        0.280                  # Scale-height of the disk (in kpc)
+      amplitude:        1.0                    # Amplitude of the bulge
+      r_1_kpc:          1.0                    # Reference radius for amplitude of the bulge (in kpc)
+      alpha:            1.8                    # Exponent of the power law of the bulge
+      r_c_kpc:          1.9                    # Cut-off radius of the bulge (in kpc)
+      potential_factors: [0.4367419745056084, 1.002641971008805, 0.022264787598364262] #Coefficients that adjust the strength of the halo (1st component), the disk (2nd component) and the bulge (3rd component)
+
+Note that the default value of the "Hubble constant" here seems odd. As it
+enters multiplicatively with the :math:`f_1` term, the absolute normalisation is
+actually not important.
+      
 How to implement your own potential
 -----------------------------------
 

examples/GravityTests/Hernquist_circularorbit/makeIC.py

@@ -16,9 +16,6 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 ################################################################################
-from galpy.potential import NFWPotential
-from galpy.orbit import Orbit
-from galpy.util import bovy_conversion
 import numpy as np
 import matplotlib.pyplot as plt
 from astropy import units
@@ -27,7 +24,7 @@ import h5py as h5
 C = 8.0
 M_200 = 2.0
 N_PARTICLES = 3
-print("Initial conditions written to 'test_nfw.hdf5'")
+print("Initial conditions written to 'circularorbitshernquist.hdf5'")
 
 pos = np.zeros((3, 3))
 pos[0, 2] = 50.0

examples/GravityTests/MWPotential2014_circularorbit/README

+# Context
+
+This example tests the MWPotential2014, a reference Milky-Way potential implemented on galpy (https://docs.galpy.org). 
+Details of the parameters can be found in galpy's paper, p.7:
+galpy: A Python Library for Galactic Dynamics, Jo Bovy (2015), Astrophys. J. Supp., 216, 29 (<arXiv/1412.3451>)
+
+# How to run this example
+
+In a terminal at the root of the "swiftsim" directory, type:
+
+`./autogen.sh`
+
+Then, configure swift to take into account external potentials. Type:
+
+`./configure --with-ext-potential=MWPotential2014`
+
+Feel free to adapt other configurations parameters depending on your system. Then, build the code by typing:
+
+`make`
+
+Finally, to run this example, open a terminal in the present directory and type:
+
+`./run.sh`
+

examples/GravityTests/MWPotential2014_circularorbit/makeIC.py

+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2023 Darwin Roduit (darwin.roduit@epfl.ch)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program 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 Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+################################################################################
+import numpy as np
+import h5py as h5
+
+box_size = 1000.0
+N_PARTICLES = 3
+print("Initial conditions written to 'circular_orbits_MW.hdf5'")
+
+pos = np.zeros((3, 3))
+pos[0, 2] = 5.0
+pos[1, 1] = 5.0
+pos[2, 0] = 30.0
+pos += np.array(
+    [box_size / 2, box_size / 2, box_size / 2]
+)  # shifts the particles to the center of the box
+vel = np.zeros((3, 3))
+vel[0, 0] = 198.5424557586175
+vel[1, 0] = 225.55900735974072
+vel[2, 1] = 188.5272441374569
+
+ids = np.array([1.0, 2.0, 3.0])
+mass = np.array([1.0, 1.0, 1.0]) * 1e-10
+
+# File
+file = h5.File("circular_orbits_MW.hdf5", "w")
+
+# Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 3.086e21
+grp.attrs["Unit mass in cgs (U_M)"] = 1.98848e43
+grp.attrs["Unit time in cgs (U_t)"] = 3.086e16
+grp.attrs["Unit current in cgs (U_I)"] = 1.0
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.0
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = box_size
+grp.attrs["NumPart_Total"] = [0, N_PARTICLES, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [0, N_PARTICLES, 0, 0, 0, 0]
+grp.attrs["Time"] = 0.0
+grp.attrs["NumFilesPerSnapshot"] = 1
+grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["Dimension"] = 3
+
+# Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = 1
+
+# Particle group
+grp1 = file.create_group("/PartType1")
+ds = grp1.create_dataset("Velocities", (N_PARTICLES, 3), "f", data=vel)
+ds = grp1.create_dataset("Masses", (N_PARTICLES,), "f", data=mass)
+ds = grp1.create_dataset("ParticleIDs", (N_PARTICLES,), "L", data=ids)
+ds = grp1.create_dataset("Coordinates", (N_PARTICLES, 3), "d", data=pos)
+
+file.close()

examples/GravityTests/MWPotential2014_circularorbit/makePlots.py

+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2023 Darwin Roduit (darwin.roduit@epfl.ch)
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program 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 Lesser General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+################################################################################
+import numpy as np
+import h5py
+import matplotlib.pyplot as plt
+
+
+#%%Functions
+
+
+def get_positions_and_time(N_snapshots, N_part, output_dir, boxsize):
+    xx = np.zeros((N_part, N_snapshots))
+    yy = np.zeros((N_part, N_snapshots))
+    zz = np.zeros((N_part, N_snapshots))
+    time = np.zeros(N_snapshots)
+    pot = np.zeros((N_part, N_snapshots))
+    for i in range(0, N_snapshots):
+        Data = h5py.File(output_dir + "/output_%04d.hdf5" % i, "r")
+        header = Data["Header"]
+        time[i] = header.attrs["Time"][0]
+        particles = Data["PartType1"]
+        positions = particles["Coordinates"]
+        xx[:, i] = positions[:, 0] - boxsize / 2.0
+        yy[:, i] = positions[:, 1] - boxsize / 2.0
+        zz[:, i] = positions[:, 2] - boxsize / 2.0
+        pot[:, i] = particles["Potentials"][:]
+    return xx, yy, zz, time, pot
+
+
+def plot_orbits(x, y, z, color, save_fig_name_suffix):
+    # Plots the orbits
+    fig, ax = plt.subplots(nrows=1, ncols=3, num=1, figsize=(12, 4.1))
+    fig.suptitle("Orbits", fontsize=15)
+    ax[0].clear()
+    ax[1].clear()
+
+    for i in range(0, N_part):
+        ax[0].plot(x[i, :], y[i, :], color[i])
+
+    ax[0].set_aspect("equal", "box")
+    ax[0].set_xlim([-35, 35])
+    ax[0].set_ylim([-35, 35])
+    ax[0].set_ylabel("y (kpc)")
+    ax[0].set_xlabel("x (kpc)")
+
+    for i in range(0, N_part):
+        ax[1].plot(x[i, :], z[i, :], col[i])
+
+    ax[1].set_aspect("equal", "box")
+    ax[1].set_xlim([-35, 35])
+    ax[1].set_ylim([-35, 35])
+    ax[1].set_ylabel("z (kpc)")
+    ax[1].set_xlabel("x (kpc)")
+    ax[1].legend(
+        [
+            "Particule 1, $R = 5$ kpc",
+            "Particule 2, $R = 5$ kpc",
+            "Particule 3, $R = 30$ kpc",
+        ]
+    )
+
+    for i in range(0, N_part):
+        ax[2].plot(y[i, :], z[i, :], col[i])
+
+    ax[2].set_aspect("equal", "box")
+    ax[2].set_xlim([-35, 35])
+    ax[2].set_ylim([-35, 35])
+    ax[2].set_ylabel("z (kpc)")
+    ax[2].set_xlabel("y (kpc)")
+    plt.tight_layout()
+    plt.savefig("circular_orbits" + save_fig_name_suffix + ".png", bbox_inches="tight")
+    plt.close()
+
+
+def plot_deviation_from_circular_orbits(x, y, z, time, color, save_fig_name_suffix):
+    # Plot of the deviation from circular orbit
+    R_1_th = 5.0  # kpc
+    R_2_th = 5.0  # kpc
+    R_3_th = 30.0  # kpc
+
+    fig2, ax2 = plt.subplots(nrows=1, ncols=2, num=2, figsize=(12, 4.5))
+    fig2.suptitle("Deviation from circular orbit", fontsize=15)
+    ax2[0].clear()
+    ax2[1].clear()
+
+    # Gather the x,y and z components of each particule into one array
+    pos_1 = np.array([x[0, :], y[0, :], z[0, :]])
+    pos_2 = np.array([x[1, :], y[1, :], z[1, :]])
+    pos_3 = np.array([x[2, :], y[2, :], z[2, :]])
+
+    # Compute the radii
+    r_1 = np.linalg.norm(pos_1, axis=0)
+    error_1 = np.abs(r_1 - R_1_th) / R_1_th * 100
+    r_2 = np.linalg.norm(pos_2, axis=0)
+    error_2 = np.abs(r_2 - R_2_th) / R_2_th * 100
+    r_3 = np.linalg.norm(pos_3, axis=0)
+    error_3 = np.abs(r_3 - R_3_th) / R_3_th * 100
+
+    ax2[0].plot(time, error_1, color[0])
+    ax2[1].plot(time, error_2, color[1])
+    ax2[1].plot(time, error_3, color[2])
+    ax2[0].set_ylabel("Deviation (\%)")
+    ax2[0].set_xlabel("Time (Gyr)")
+    ax2[1].set_ylabel("Deviation (\%)")
+    ax2[1].set_xlabel("Time (Gyr)")
+    ax2[0].legend(["Particule 1, $R = 5$ kpc"])
+    ax2[1].legend(["Particule 2, $R = 5$ kpc", "Particule 3, $R = 30$ kpc"])
+
+    plt.tight_layout()
+    plt.savefig("deviation" + save_fig_name_suffix + ".png", bbox_inches="tight")
+    plt.close()
+    return r_1, r_2, r_3
+
+
+def plot_deviation_from_original_data(r_1, r_2, r_3, time, color, save_fig_name_suffix):
+    """Make a comparison with the obtained data and ours to check nothing is broken."""
+    filename = "original_radii.txt"
+    r_1_original, r_2_original, r_3_original = np.loadtxt(filename)
+
+    # Plots the deviation wrt the original data
+    fig3, ax3 = plt.subplots(nrows=1, ncols=3, num=3, figsize=(12, 4.3))
+    fig3.suptitle("Deviation from the original data", fontsize=15)
+    ax3[0].clear()
+    ax3[1].clear()
+    ax3[2].clear()
+
+    error_1 = np.abs(r_1 - r_1_original) / r_1_original * 100
+    error_2 = np.abs(r_2 - r_2_original) / r_2_original * 100
+    error_3 = np.abs(r_3 - r_3_original) / r_3_original * 100
+
+    ax3[0].plot(time, error_1, col[0])
+    ax3[1].plot(time, error_2, col[1])
+    ax3[2].plot(time, error_3, col[2])
+    ax3[0].set_ylabel("Deviation (\%)")
+    ax3[0].set_xlabel("Time (Gyr)")
+    ax3[1].set_ylabel("Deviation (\%)")
+    ax3[1].set_xlabel("Time (Gyr)")
+    ax3[2].set_ylabel("Deviation (\%)")
+    ax3[2].set_xlabel("Time (Gyr)")
+    ax3[0].legend(["Particule 1, $R = 5$ kpc"])
+    ax3[1].legend(["Particule 2, $R = 5$ kpc"])
+    ax3[2].legend(["Particule 3, $R = 30$ kpc"])
+    plt.tight_layout()
+    plt.savefig(
+        "deviation_from_original_data" + save_fig_name_suffix + ".png",
+        bbox_inches="tight",
+    )
+    plt.close()
+
+
+#%%Plots the orbits, the deviation from the circular orbit and the deviation from the original precomputed data
+# Notice that in this examples, the ouputs are set in suitable units in the parameters files.
+
+# General parameters
+N_snapshots = 201
+N_part = 3
+boxsize = 1000.0  # kpc
+col = ["b", "r", "c", "y", "k"]
+
+# First type of units (kpc)
+output_dir = "output_1"
+save_fig_name_suffix = "_simulation_kpc"
+x_1, y_1, z_1, time_1, pot_1 = get_positions_and_time(
+    N_snapshots, N_part, output_dir, boxsize
+)
+plot_orbits(x_1, y_1, z_1, col, save_fig_name_suffix)
+r_11, r_21, r_31 = plot_deviation_from_circular_orbits(
+    x_1, y_1, z_1, time_1, col, save_fig_name_suffix
+)
+plot_deviation_from_original_data(r_11, r_21, r_31, time_1, col, save_fig_name_suffix)
+
+# Second type of units (Mpc) (no need for units conversion to kpc, this is already done by swift in the snapshots)
+output_dir = "output_2"
+save_fig_name_suffix = "_simulation_Mpc"
+x_2, y_2, z_2, time_2, pot_2 = get_positions_and_time(
+    N_snapshots, N_part, output_dir, boxsize
+)
+plot_orbits(x_2, y_2, z_2, col, save_fig_name_suffix)
+r_12, r_22, r_32 = plot_deviation_from_circular_orbits(
+    x_2, y_2, z_2, time_2, col, save_fig_name_suffix
+)
+# plot_deviation_from_original_data(r_12, r_22, r_32, time_2, col, save_fig_name_suffix) #does not make sense since the original data are in kpc, not in Mpc
+
+#%%Saves our data to be the reference ones (precomputed)
+# Uncomment only if corrections of the precomputed data must occur !
+# Original data :  If some corrections occur in the potential default parameters, allows to correct
+# the data
+# filename = "original_radii.txt"
+# np.savetxt(filename, (r_1, r_2, r_3))

examples/GravityTests/MWPotential2014_circularorbit/original_radii.txt

+5.000000000000000000e+00 5.001166599459545559e+00 5.020579887129778207e+00 5.100189116397097600e+00 5.299839849332817820e+00 5.578710433997593476e+00 5.783030230255659987e+00 5.827598846292087131e+00 5.682565987909128147e+00 5.357447130085647657e+00 4.912578397139871988e+00 4.497394792994665380e+00 4.397021356282944105e+00 4.775431657567319910e+00 5.300932384346391579e+00 5.731042038011933570e+00 5.965091504370781728e+00 5.970489901937201971e+00 5.752272156326732500e+00 5.356957312236751534e+00 4.906750290982262008e+00 4.657121571444537089e+00 4.714457558888792477e+00 4.933301629890139317e+00 5.176061235968707486e+00 5.362349356536657119e+00 5.461122477067911873e+00 5.485087743207845534e+00 5.485339856546893600e+00 5.523077935189706800e+00 5.537107938749801228e+00 5.464804174270136095e+00 5.285227013585691580e+00 5.020384810243945672e+00 4.743263219255731578e+00 4.592936367545448206e+00 4.749715571294080618e+00 5.194699542807357240e+00 5.642069433457655769e+00 5.933223989683256150e+00 6.005807399351317244e+00 5.844582655179585196e+00 5.470114662333674715e+00 4.956958096079127962e+00 4.496399606097051027e+00 4.436328410208022710e+00 4.769513345769095025e+00 5.203928883967520846e+00 5.560443291869900051e+00 5.758575039733128342e+00 5.777233497154392161e+00 5.634900936928413984e+00 5.399734327667077949e+00 5.204348571913461008e+00 5.108383765639208818e+00 5.057693052077008034e+00 5.015035997691488667e+00 4.970823689765934361e+00 4.947419029227536846e+00 4.995207789090853723e+00 5.182270748184618192e+00 5.494607608495452489e+00 5.759878945928849525e+00 5.871757751107428369e+00 5.788023974259131066e+00 5.507664932242153810e+00 5.073360594415489366e+00 4.604596875155709590e+00 4.356042545916934294e+00 4.609781734172386791e+00 5.121533352837495556e+00 5.594559699714237055e+00 5.895503032313782477e+00 5.976985446160555604e+00 5.833832600072278218e+00 5.499788789299549840e+00 5.068750025511791435e+00 4.754624075584575671e+00 4.723949943832726817e+00 4.875445028387068724e+00 5.081882737105429371e+00 5.259925739770850761e+00 5.372637274737933843e+00 5.425242400507194418e+00 5.460118287161413342e+00 5.536220897651173090e+00 5.595059701109168948e+00 5.560241134615987235e+00 5.401506769780430872e+00 5.130567420218448582e+00 4.809371855699721365e+00 4.569652033714569406e+00 4.612229578661045437e+00 5.019851467780982901e+00 5.505814358798329700e+00 5.865468786224973918e+00 6.016174230579117932e+00 5.931623634176135695e+00 5.622177861123924814e+00 5.141870044592171674e+00 4.638547645705167533e+00 4.428307819955180413e+00 4.654500832829794099e+00 5.055593990419128048e+00 5.427492129376623176e+00 5.668492207620567491e+00 5.744061565586670248e+00 5.662423276901259683e+00 5.476861486048790084e+00 5.299864737421422056e+00 5.202470997364078364e+00 5.134308457492519295e+00 5.058899583759844276e+00 4.970512274614167225e+00 4.896917303931519250e+00 4.898286265130379569e+00 5.054476498797095374e+00 5.384729626924356083e+00 5.707038692608169761e+00 5.888268525975762557e+00 5.871970152607735471e+00 5.647189263247745394e+00 5.241550304616231060e+00 4.747405841925699477e+00 4.378927045831982667e+00 4.474186584819125123e+00 4.940610149752312275e+00 5.439852850044650800e+00 5.799438267533770563e+00 5.953178490303065118e+00 5.885318802304978725e+00 5.618141372733455263e+00 5.224255162563557597e+00 4.875392411458151720e+00 4.768956166135601471e+00 4.850718648821461976e+00 5.009868149351704325e+00 5.165493314050841889e+00 5.278258318859305298e+00 5.348100423409463566e+00 5.410955922599304913e+00 5.522902505871598144e+00 5.630663735706390227e+00 5.642294107886619337e+00 5.517118052537815842e+00 5.256074072268789088e+00 4.908246395925647043e+00 4.590870106450615218e+00 4.510269017110958956e+00 4.841676564404478356e+00 5.347528554683028013e+00 5.768526167155091144e+00 5.995637396520894669e+00 5.989832219840335092e+00 5.751871797307715362e+00 5.320044178608156926e+00 4.806704151149128634e+00 4.474934801273977136e+00 4.576629640621881379e+00 4.921538717217767811e+00 5.290768217179389055e+00 5.562191457746648915e+00 5.686534351951245014e+00 5.661795715468069368e+00 5.528388603164072101e+00 5.380407491001560238e+00 5.294789296099782661e+00 5.223565412511718797e+00 5.127422326009080322e+00 5.001198016563644266e+00 4.875299081804920220e+00 4.818639828265833813e+00 4.925709984592579360e+00 5.252798132458202041e+00 5.624850406729384744e+00 5.875379086292169539e+00 5.930742494670786513e+00 5.769927451901470938e+00 5.407923458724494203e+00 4.913709661908932880e+00 4.459744560941122593e+00 4.382798389554843688e+00 4.768537003581841738e+00 5.272907504565521819e+00 5.680248961594461576e+00 5.900602226243958803e+00 5.906378074467261818e+00 5.709143330590744192e+00 5.364799048397715708e+00 5.008887298284076195e+00 4.843398709646222322e+00 4.859874698091035228e+00 4.965584565534517658e+00 5.086695693063617796e+00 5.185431020385737000e+00 5.259341369321015591e+00 5.340599814904859244e+00 5.482731832086808588e+00 5.640753798791949158e+00 5.705643537489685002e+00 5.624938857189243357e+00
+5.000000000000000000e+00 4.999888359284162753e+00 5.000641684297063350e+00 5.000994447130620024e+00 5.001675693602675388e+00 5.001624794631814197e+00 5.001558065482559989e+00 5.000888414115478575e+00 5.000403424123103235e+00 4.999857363290737489e+00 4.999815117935768605e+00 5.000030078010763823e+00 5.000590428170696988e+00 5.001185604320395939e+00 5.001578829707598928e+00 5.001700786476619420e+00 5.001329335974549650e+00 5.000875735587428927e+00 5.000172435947884964e+00 4.999938404239019363e+00 4.999745387639388028e+00 5.000283699222444866e+00 5.000643388144556667e+00 5.001470334121859374e+00 5.001565799279720004e+00 5.001856736927319247e+00 5.001184047933382892e+00 5.000858585776699172e+00 5.000045524839611844e+00 4.999938201977301766e+00 4.999789233051740389e+00 5.000353429170261599e+00 5.000792119848084027e+00 5.001452348355899247e+00 5.001632422736415506e+00 5.001608535971686642e+00 5.001127158676252904e+00 5.000553905860360970e+00 5.000029329624915420e+00 4.999762162118537745e+00 4.999946677627926661e+00 5.000335136347326070e+00 5.001035783553454017e+00 5.001422330001844330e+00 5.001779305392438424e+00 5.001442290734069651e+00 5.001159754497918719e+00 5.000348863737741567e+00 5.000122347527026889e+00 4.999700132208437786e+00 5.000215353627103454e+00 5.000418901608895794e+00 5.001274128363529137e+00 5.001475223091539135e+00 5.001811379870214402e+00 5.001359505875552891e+00 5.000995546297854588e+00 5.000250370590015869e+00 4.999942233798068791e+00 4.999755726773103426e+00 5.000113010258659507e+00 5.000609015959638803e+00 5.001219673802144605e+00 5.001614354202936319e+00 5.001636066331053954e+00 5.001364334692166125e+00 5.000733606126373409e+00 5.000258441419644484e+00 4.999780701257404480e+00 4.999929042417343261e+00 5.000122602854585097e+00 5.000898089680752889e+00 5.001250786342685828e+00 5.001830340890685456e+00 5.001536553836007926e+00 5.001450960015414893e+00 5.000564462723501968e+00 5.000273841477640424e+00 4.999729946197625452e+00 5.000019708256767004e+00 5.000238197444725863e+00 5.000993609020794928e+00 5.001370698158823735e+00 5.001741598978786563e+00 5.001520071329959016e+00 5.001141309344722785e+00 5.000492961311684859e+00 5.000007097885235119e+00 4.999789471734320756e+00 4.999926978056772775e+00 5.000452501234716074e+00 5.000982517444656494e+00 5.001568938123119423e+00 5.001630494678826366e+00 5.001581127018102535e+00 5.000918624733289519e+00 5.000520523906761383e+00 4.999849635612912735e+00 4.999963116585512957e+00 4.999945908278064621e+00 5.000768878036669918e+00 5.001060851264059082e+00 5.001749243021976099e+00 5.001594604383085674e+00 5.001528201965431464e+00 5.000787308820795474e+00 5.000362783584799686e+00 4.999809889944675234e+00 4.999870251224309392e+00 5.000081273895164991e+00 5.000706925098050704e+00 5.001234607158520085e+00 5.001628277065568895e+00 5.001642902217743547e+00 5.001272369050412614e+00 5.000752316183323387e+00 5.000115144777662834e+00 4.999877063130510280e+00 4.999785278748652750e+00 5.000314975961055453e+00 5.000732923535549546e+00 5.001486668371692978e+00 5.001580864064437471e+00 5.001767526649143925e+00 5.001102524724919896e+00 5.000815405699214899e+00 4.999977098680566279e+00 4.999965430706085456e+00 4.999820211656420987e+00 5.000468706648866224e+00 5.000862674405269992e+00 5.001540997553139967e+00 5.001619721607925229e+00 5.001589058085179396e+00 5.001022878428104335e+00 5.000494347648318794e+00 4.999950437180141094e+00 4.999782011985283603e+00 4.999968771211793950e+00 5.000436397880152484e+00 5.001086192651295725e+00 5.001485312783913173e+00 5.001737933662662172e+00 5.001395316031677751e+00 5.001032550774532126e+00 5.000272753974341455e+00 5.000032632149568279e+00 4.999711660403833413e+00 5.000227368341616518e+00 5.000504073344837153e+00 5.001396382624039738e+00 5.001505342442887248e+00 5.001831185827574799e+00 5.001279357006295001e+00 5.000939059512067075e+00 5.000156498511878489e+00 4.999939419656115547e+00 4.999760777523306388e+00 5.000213242684194981e+00 5.000677422077679068e+00 5.001316266444979952e+00 5.001612758519516255e+00 5.001623006669117188e+00 5.001255205432448392e+00 5.000655856805036770e+00 5.000151079135440213e+00 4.999768706765014059e+00 4.999923501590752828e+00 5.000205665288950385e+00 5.000941154591838078e+00 5.001315826587424240e+00 5.001794596942993110e+00 5.001491317612160259e+00 5.001316595758519057e+00 5.000470387348008749e+00 5.000254580857065534e+00 4.999712075477516393e+00 5.000102976947414568e+00 5.000306391555606744e+00 5.001111012410059509e+00 5.001407461804255661e+00 5.001773556648744012e+00 5.001449382725613013e+00 5.001085296333972252e+00 5.000387772118549456e+00 4.999982482345314061e+00 4.999767868097886314e+00 5.000003634425123522e+00 5.000507888369177145e+00 5.001079803140623170e+00 5.001579362536846318e+00 5.001633342805087423e+00 5.001483416394396642e+00 5.000841451375783286e+00 5.000401977832478195e+00 4.999818075316714072e+00 4.999936405692012364e+00 5.000014185319412441e+00 5.000809312551328212e+00 5.001135437752519231e+00 5.001827481510867202e+00
+3.000000000000000000e+01 2.999779458340697147e+01 2.999913108927551875e+01 2.999822513999371054e+01 2.999869101241313629e+01 2.999908157608862069e+01 2.999867255062164162e+01 3.000035432115671696e+01 2.999906436031323764e+01 3.000130726868057707e+01 2.999985124246509116e+01 3.000120011581939750e+01 3.000101355345713472e+01 3.000145740135981853e+01 3.000252851377030439e+01 3.000205476548366335e+01 3.000437014149219905e+01 3.000296575532702192e+01 3.000506656450175669e+01 3.000416345420759257e+01 3.000531256555992954e+01 3.000561911665013071e+01 3.000580235611645818e+01 3.000730436832045811e+01 3.000650795996002174e+01 3.000919171647245420e+01 3.000740253616102748e+01 3.000908722559980291e+01 3.000846124470268705e+01 3.000913539394623442e+01 3.000966175811804604e+01 3.000931550053518393e+01 3.001098523842369659e+01 3.000961070226721716e+01 3.001169411489794570e+01 3.001000798292042404e+01 3.001105493839978067e+01 3.001049866425505996e+01 3.001050661892987037e+01 3.001107862320722219e+01 3.001004719415356803e+01 3.001174860804977840e+01 3.000967995094207907e+01 3.001106890886319079e+01 3.000941297065721614e+01 3.000977242153797064e+01 3.000925845687551785e+01 3.000859589580250741e+01 3.000923251594434049e+01 3.000755748209316565e+01 3.000935520671618661e+01 3.000667892784183266e+01 3.000748164415950114e+01 3.000598473874775962e+01 3.000580528773574329e+01 3.000550129646441633e+01 3.000435340448430566e+01 3.000525642811016169e+01 3.000315424345881254e+01 3.000455578094075548e+01 3.000223634049040911e+01 3.000270501649730548e+01 3.000162796403852283e+01 3.000117720180458747e+01 3.000135582863928008e+01 2.999999786615864394e+01 3.000144392723029085e+01 2.999919002924884381e+01 3.000046804618323293e+01 2.999877353580113848e+01 2.999916917874825728e+01 2.999876482042754233e+01 2.999828485654238719e+01 2.999917646486512979e+01 2.999782544062704659e+01 3.000001627767003498e+01 2.999779633743827034e+01 2.999911844063658251e+01 2.999819830488571526e+01 2.999865030911109187e+01 2.999902751315119787e+01 2.999860581120435299e+01 3.000027557069040540e+01 2.999897424318053041e+01 3.000120647142659891e+01 2.999974052798087598e+01 3.000108036609949735e+01 3.000088565466134227e+01 3.000132237211635555e+01 3.000238750303375213e+01 3.000190931347448497e+01 3.000422185867911296e+01 3.000281580647389745e+01 3.000491608599069338e+01 3.000401361861365146e+01 3.000516478912791385e+01 3.000547462841453950e+01 3.000566249461711266e+01 3.000717035022314505e+01 3.000638099269628611e+01 3.000907299236889969e+01 3.000729313447100211e+01 3.000898818458806971e+01 3.000837351407666986e+01 3.000905988756514375e+01 3.000959933442729266e+01 3.000926695882656148e+01 3.001095125287651655e+01 3.000959180150701044e+01 3.001169059343623147e+01 3.001002003559164066e+01 3.001108278198048041e+01 3.001054248304392047e+01 3.001056653531087548e+01 3.001115458351709364e+01 3.001013901847990439e+01 3.001185589737429282e+01 3.000980246935861828e+01 3.001120595492282561e+01 3.000956393020903334e+01 3.000993653574142783e+01 3.000943488286652894e+01 3.000878378812352665e+01 3.000943101649924571e+01 3.000776555887522434e+01 3.000957175115255993e+01 3.000690282240585205e+01 3.000771163940339292e+01 3.000621963347996513e+01 3.000604381912500784e+01 3.000574221323055824e+01 3.000459546464529481e+01 3.000549825809542170e+01 3.000339446894059847e+01 3.000479306306118232e+01 3.000246913202072463e+01 3.000293187427373454e+01 3.000184761572847592e+01 3.000138830507564691e+01 3.000155700520485880e+01 3.000018799611447662e+01 3.000162189224056064e+01 2.999935478803817901e+01 3.000061862764304976e+01 2.999890904656334101e+01 2.999928887416733758e+01 2.999886807861885174e+01 2.999837114075570099e+01 2.999924530142823542e+01 2.999787658246655297e+01 3.000004954612104058e+01 2.999781168046782653e+01 2.999911595907233774e+01 2.999817825072545574e+01 2.999861297335620591e+01 2.999897334516337466e+01 2.999853524782260195e+01 3.000018922740121141e+01 2.999887294236322077e+01 3.000109111179866161e+01 2.999961196637686811e+01 3.000093950505366180e+01 3.000073363009094862e+01 3.000116029442814991e+01 3.000221660273533431e+01 3.000173085754470392e+01 3.000403703825239532e+01 3.000262587834090766e+01 3.000472258970311401e+01 3.000381811956868461e+01 3.000496875334294344e+01 3.000527961996388271e+01 3.000546996417365975e+01 3.000698169725789910e+01 3.000619758041846197e+01 3.000889621429398701e+01 3.000712430259497765e+01 3.000882861345186825e+01 3.000822447317817065e+01 3.000892257912635230e+01 3.000947488798068008e+01 3.000915634131698795e+01 3.001085531632530490e+01 3.000951134895363381e+01 3.001162660663587900e+01 3.000997326126640274e+01 3.001105384355043526e+01 3.001053154124043232e+01 3.001057371847517885e+01 3.001118006937515403e+01 3.001018279679330192e+01 3.001191804078371206e+01 3.000988258949455556e+01 3.001130365826893609e+01 3.000967875120074169e+01 3.001006775551087458e+01 3.000958182273279817e+01 3.000894558530716338e+01 3.000960665058126864e+01 3.000795400716476991e+01 3.000977178013717150e+01

examples/GravityTests/MWPotential2014_circularorbit/params_unit_1.yml

+# Define the system of units to use internally.
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.98848e+43 # 10^10 Solar masses
+  UnitLength_in_cgs:   3.086e+21 # kpc
+  UnitVelocity_in_cgs: 1e5       # 1 km / s in cm/s
+  UnitCurrent_in_cgs:  1         # Amperes
+  UnitTemp_in_cgs:     1         # Kelvin
+
+# Parameters governing the time integration (Set dt_min and dt_max to the same value for a fixed time-step run.)
+TimeIntegration:
+  time_begin:          0.      # The starting time of the simulation (in internal units).
+  time_end:            2.0     # The end time of the simulation (in internal units).
+  dt_min:              1e-10    # The minimal time-step size of the simulation (in internal units).
+  dt_max:              1e0    # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:            output_1/output  # Common part of the name of output files
+  time_first:          0.      # Time of the first output (in internal units)
+  delta_time:          1e-2    # Time difference between consecutive outputs (in internal units)
+  UnitMass_in_cgs:     1.98848e+43
+  UnitLength_in_cgs:   3.086e+21
+  UnitVelocity_in_cgs: 1e5
+  UnitCurrent_in_cgs:  1
+  UnitTemp_in_cgs:     1
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          1.0    # Time between statistics output
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:          circular_orbits_MW.hdf5 # The file to read
+  periodic:           1
+
+# NFW_MN_PSC potential parameters
+MWPotential2014Potential:
+  useabspos:       0        # 0 -> positions based on centre, 1 -> absolute positions 
+  position:        [0.,0.,0.]    #Centre of the potential with respect to centre of the box
+  timestep_mult:   0.005     # Dimensionless pre-factor for the time-step condition
+  epsilon:         0.001      # Softening size (internal units)
+  #The following parameters are the default ones.
+  # concentration:    9.823403437774843
+  # M_200_Msun:       147.41031542774076e10
+  # H:                127.78254614201471e-2
+  # Mdisk_Msun:       6.8e10
+  # Rdisk_kpc:        3.0
+  # Zdisk_kpc:        0.280
+  # amplitude:        1.0
+  # r_1_kpc:          1.0
+  # alpha:            1.8
+  # r_c_kpc:          1.9
+  # potential_factors: [0.4367419745056084, 1.002641971008805, 0.022264787598364262]

examples/GravityTests/MWPotential2014_circularorbit/params_unit_2.yml

+# Define the system of units to use internally.
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.98848e+43 # 10^10 Solar masses
+  UnitLength_in_cgs:   3.086e+24 # Mpc
+  UnitVelocity_in_cgs: 1e5       # 1 km / s in cm/s
+  UnitCurrent_in_cgs:  1         # Amperes
+  UnitTemp_in_cgs:     1         # Kelvin
+
+# Parameters governing the time integration (Set dt_min and dt_max to the same value for a fixed time-step run.)
+TimeIntegration:
+  time_begin:          0.      # The starting time of the simulation (in internal units).
+  time_end:            2.0e-3     # The end time of the simulation (in internal units).
+  dt_min:              1e-13    # The minimal time-step size of the simulation (in internal units).
+  dt_max:              1e-3    # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:            output_2/output  # Common part of the name of output files
+  time_first:          0.      # Time of the first output (in internal units)
+  delta_time:          1e-5    # Time difference between consecutive outputs (in internal units)
+  UnitMass_in_cgs:     1.98848e+43
+  UnitLength_in_cgs:   3.086e+21
+  UnitVelocity_in_cgs: 1e5
+  UnitCurrent_in_cgs:  1
+  UnitTemp_in_cgs:     1
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          1.0    # Time between statistics output
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:          circular_orbits_MW.hdf5 # The file to read
+  periodic:           1
+
+# NFW_MN_PSC potential parameters
+MWPotential2014Potential:
+  useabspos:       0        # 0 -> positions based on centre, 1 -> absolute positions 
+  position:        [0.,0.,0.]    #Centre of the potential with respect to centre of the box
+  timestep_mult:   0.0005     # Dimensionless pre-factor for the time-step condition
+  epsilon:         0.001e-3      # Softening size (internal units)
+  #The following parameters are the default ones.
+  # concentration:    9.823403437774843
+  # M_200_Msun:       147.41031542774076e10
+  # H:                127.78254614201471e-2
+  # Mdisk_Msun:       6.8e10
+  # Rdisk_kpc:        3.0
+  # Zdisk_kpc:        0.280
+  # amplitude:        1.0
+  # r_1_kpc:          1.0
+  # alpha:            1.8
+  # r_c_kpc:          1.9
+  # potential_factors: [0.4367419745056084, 1.002641971008805, 0.022264787598364262]

examples/GravityTests/MWPotential2014_circularorbit/run.sh

+#!/bin/bash
+
+#Creates a directory for the outputs
+DIR=output_1 #First test of units conversion
+if [ -d "$DIR" ];
+then
+    echo "$DIR directory exists. Its content will be removed."
+    rm $DIR/output_*
+else
+    echo "$DIR directory does not exists. It will be created."
+    mkdir $DIR
+fi
+
+DIR=output_2 #Second test of units conversion
+if [ -d "$DIR" ];
+then
+    echo "$DIR directory exists. Its content will be removed."
+    rm $DIR/output_*
+else
+    echo "$DIR directory does not exists. It will be created."
+    mkdir $DIR
+fi
+
+#Clears the previous figures
+echo "Clearing existing figures."
+if [ -f "circular_orbits_simulation_kpc.png" ];
+then
+    rm circular_orbits_simulation_kpc.png
+fi
+
+if [ -f "circular_orbits_simulation_Mpc.png" ];
+then
+    rm circular_orbits_simulation_Mpc.png
+fi
+
+if [ -f "deviation_simulation_kpc.png" ];
+then
+    rm deviation_simulation_kpc.png
+fi
+if [ -f "deviation_simulation_Mpc.png" ];
+then
+    rm deviation_simulation_Mpc.png
+fi
+
+if [ -f "deviation_from_original_data_simulation_kpc.png" ];
+then
+    rm deviation_from_original_data_simulation_kpc.png
+fi
+
+if [ -f "deviation_from_original_data_simulation_Mpc.png" ];
+then
+    rm deviation_from_original_data_simulation_Mpc.png
+fi
+
+#Clears the IC file
+if [ -f "circular_orbits_MW.hdf5" ];
+then
+    rm circular_orbits_MW.hdf5
+fi
+
+
+#Generates the initial conditions
+echo "Generate initial conditions for circular orbits"
+if command -v python3 &>/dev/null; then
+    python3 makeIC.py
+else
+    python3 makeIC.py
+fi
+
+#Runs the simulation
+# self gravity G, external potential g, hydro s, threads t and high verbosity v
+echo "Starting the simulation with the first type of units (kpc)... Look at output_1.log for the simulation details."
+../../../swift --external-gravity --threads=8 params_unit_1.yml 2>&1 > output_1.log
+echo "Simulation ended."
+
+echo "Starting the simulation with the second type of units (Mpc)... Look at output_2.log for the simulation details."
+../../../swift --external-gravity --threads=8 params_unit_2.yml 2>&1 > output_2.log
+echo "Simulation ended."
+
+#Saves the plots
+echo "Save plots of the circular orbits and of the errors"
+if command -v python3 &>/dev/null; then
+    python3 makePlots.py
+else 
+    python3 makePlots.py
+fi

examples/parameter_example.yml

@@ -420,6 +420,23 @@ SineWavePotential:
   timestep_limit:   1.      # Time-step dimensionless pre-factor.
   growth_time:      0.      # (Optional) Time for the potential to grow to its final size.
 
+# MWPotential2014 potential
+MWPotential2014Potential:
+  useabspos:        0          # 0 -> positions based on centre, 1 -> absolute positions
+  position:         [0.,0.,0.] # Location of centre of potential with respect to centre of the box (if 0) otherwise absolute (if 1) (internal units)
+  timestep_mult:    0.005      # Dimensionless pre-factor for the time-step condition, basically determines the fraction of the orbital time we use to do the time integration
+  epsilon:          0.001      # Softening size (internal units)
+  concentration:    9.823403437774843      # concentration of the Halo
+  M_200_Msun:       147.41031542774076e10  # M200 of the galaxy disk (in M_sun)
+  H:                1.2778254614201471     # Hubble constant in units of km/s/Mpc
+  Mdisk_Msun:       6.8e10                 # Mass of the disk (in M_sun)
+  Rdisk_kpc:        3.0                    # Effective radius of the disk (in kpc)
+  Zdisk_kpc:        0.280                  # Scale-height of the disk (in kpc)
+  amplitude:        1.0                    # Amplitude of the bulge
+  r_1_kpc:          1.0                    # Reference radius for amplitude of the bulge (in kpc)
+  alpha:            1.8                    # Exponent of the power law of the bulge
+  r_c_kpc:          1.9                    # Cut-off radius of the bulge (in kpc)
+  potential_factors: [0.4367419745056084, 1.002641971008805, 0.022264787598364262] # Coefficients that adjust the strength of the halo (1st component), the disk (2nd component) and the bulge (3rd component)
 
 # Parameters related to entropy floors    ----------------------------------------------
 

src/potential.h

@@ -42,6 +42,8 @@
 #include "./potential/nfw/potential.h"
 #elif defined(EXTERNAL_POTENTIAL_NFW_MN)
 #include "./potential/nfw_mn/potential.h"
+#elif defined(EXTERNAL_POTENTIAL_MWPotential2014)
+#include "./potential/MWPotential2014/potential.h"
 #elif defined(EXTERNAL_POTENTIAL_DISC_PATCH)
 #include "./potential/disc_patch/potential.h"
 #elif defined(EXTERNAL_POTENTIAL_SINE_WAVE)

src/potential/MWPotential2014/potential.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2023  Darwin Roduit
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program 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 Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#ifndef SWIFT_POTENTIAL_MWPotential2014_H
+#define SWIFT_POTENTIAL_MWPotential2014_H
+
+/* Config parameters. */
+#include <config.h>
+
+/* Some standard headers. */
+#include <float.h>
+#include <math.h>
+
+/* Local includes. */
+#include "error.h"
+#include "gravity.h"
+#include "parser.h"
+#include "part.h"
+#include "physical_constants.h"
+#include "space.h"
+#include "units.h"
+
+#ifdef HAVE_LIBGSL
+#include <gsl/gsl_sf_gamma.h>
+#endif
+
+/**
+ * @brief External Potential Properties - MWPotential2014 composed by
+ * NFW + Miyamoto-Nagai + Power Spherical cut-off potentials
+ *
+ * halo --> rho_NFW(r) = rho_0 / ( (r/R_s)*(1+r/R_s)^2 )
+ * disk --> Phi_MN(R,z) = -G * Mdisk / (R^2 + (Rdisk +
+ * (z^2+Zdisk^2)^1/2)^2)^(1/2) bulge --> rho_PSC(r) =
+ * amplitude*(r_1/r)^alpha*exp(-(r/r_c)^2)
+ *
+ * We however parametrise this in terms of c and virial_mass, Mdisk, Rdisk
+ * and Zdisk. Also, each potential is given a contribution amplitude such that
+ * the resulting potential is:
+ *      Phi_tot = f_1 * Phi_NFW + f_2 * Phi_MN + f_3 * Phi_PSC,
+ * with f_1, f_2 and f_3 contained in the array f.
+ *
+ * This potential is inspired by the following article:
+ * galpy: A Python Library for Galactic Dynamics, Jo Bovy (2015),
+ * Astrophys. J. Supp., 216, 29 (arXiv/1412.3451).
+ */
+struct external_potential {
+
+  /*! Position of the centre of potential */
+  double x[3];
+
+  /*! The scale radius of the NFW potential */
+  double r_s;
+
+  /*! The pre-factor \f$ 4 \pi G \rho_0 \r_s^3 \f$ */
+  double pre_factor;
+
+  /*! Hubble parameter */
+  double H;
+
+  /*! The concentration parameter */
+  double c_200;
+
+  /*! The virial mass */
+  double M_200;
+
+  /*! Disk Size */
+  double Rdisk;
+
+  /*! Disk height */
+  double Zdisk;
+
+  /*! Disk Mass */
+  double Mdisk;
+
+  /*! Amplitude for the PSC potential */
+  double amplitude;
+
+  /*! Reference radius for amplitude */
+  double r_1;
+
+  /*! Inner power */
+  double alpha;
+
+  /*! Cut-off radius */
+  double r_c;
+
+  /*! Contribution of each potential : f[0]*NFW + f[1]*MN + f[2]*PSP */
+  double f[3];
+
+  /*! Prefactor \f$ 2 \pi amplitude r_1^\alpha r_c^(3-\alpha) \f$ */
+  double prefactor_psc_1;
+
+  /*! Prefactor \f$ 2 \pi amplitude r_1^\alpha r_c^(2-\alpha) \f$ */
+  double prefactor_psc_2;
+
+  /*! Gamma function evaluation \f$ \Gamma((3-\alpha)/2 \f$ */
+  double gamma_psc;
+
+  /*! Time-step condition pre_factor, this factor is used to multiply times the
+   * orbital time, so in the case of 0.01 we take 1% of the orbital time as
+   * the time integration steps */
+  double timestep_mult;
+
+  /*! Minimum time step based on the orbital time at the softening times
+   * the timestep_mult */
+  double mintime;
+
+  /*! Common log term \f$ \ln(1+c_{200}) - \frac{c_{200}}{1 + c_{200}} \f$ */
+  double log_c200_term;
+
+  /*! Softening length */
+  double eps;
+};
+
+/**
+ * @brief Computes the time-step due to the acceleration from the NFW + MN + PSC
+ * potential as a fraction (timestep_mult) of the circular orbital time of that
+ * particle.
+ *
+ * @param time The current time.
+ * @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 float external_gravity_timestep(
+    double time, const struct external_potential* restrict potential,
+    const struct phys_const* restrict phys_const,
+    const struct gpart* restrict g) {
+
+#ifdef HAVE_LIBGSL
+
+  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];
+
+  const float R2 = dx * dx + dy * dy;
+  const float r = sqrtf(R2 + dz * dz + potential->eps * potential->eps);
+
+  /* Vcirc for NFW */
+  const float M_NFW = potential->pre_factor * (logf(1.0f + r / potential->r_s) -
+                                               r / (r + potential->r_s));
+  const float Vcirc_NFW = sqrtf((phys_const->const_newton_G * M_NFW) / r);
+
+  /* Now for MN */
+  const float R = sqrtf(R2);
+  const float sqrt_term = sqrtf(dz * dz + potential->Zdisk * potential->Zdisk);
+  const float MN_denominator =
+      powf(R2 + powf(potential->Rdisk + sqrt_term, 2.0f), 1.5f);
+  const float dPhi_dR_MN = potential->Mdisk * R / MN_denominator;
+  const float dPhi_dz_MN = potential->Mdisk * dz *
+                           (potential->Rdisk + sqrt_term) /
+                           (sqrt_term * MN_denominator);
+  const float Vcirc_MN = sqrtf(phys_const->const_newton_G * R * dPhi_dR_MN +
+                               phys_const->const_newton_G * dz * dPhi_dz_MN);
+
+  /* Now for PSC */
+  const float r2 = r * r;
+  const float M_psc =
+      potential->prefactor_psc_1 *
+      (potential->gamma_psc -
+       gsl_sf_gamma_inc(1.5f - 0.5f * potential->alpha,
+                        r2 / (potential->r_c * potential->r_c)));
+  const float Vcirc_PSC = sqrtf(phys_const->const_newton_G * M_psc / r);
+
+  /* Total circular velocity */
+  const float Vcirc = sqrtf(potential->f[0] * Vcirc_NFW * Vcirc_NFW +
+                            potential->f[1] * Vcirc_MN * Vcirc_MN +
+                            potential->f[2] * Vcirc_PSC * Vcirc_PSC);
+
+  const float period = 2.0f * M_PI * r / Vcirc;
+
+  /* Time-step as a fraction of the circular period */
+  const float time_step = potential->timestep_mult * period;
+
+  return max(time_step, potential->mintime);
+
+#else
+  error("Code not compiled with GSL. Can't compute MWPotential2014.");
+  return 0.0;
+#endif
+}
+
+/**
+ * @brief Computes the gravitational acceleration from an NFW Halo potential +
+ * MN disk + PSC bulge.
+ *
+ * Note that the accelerations are multiplied by Newton's G constant
+ * later on.
+ *
+ * a_x = 4 pi \rho_0 r_s^3 ( 1/((r+rs)*r^2) - log(1+r/rs)/r^3) * x
+ *      - dphi_PSC/dr*x/r
+ * a_y = 4 pi \rho_0 r_s^3 ( 1/((r+rs)*r^2) - log(1+r/rs)/r^3) * y
+ *      - dphi_PSC/dr*y/r
+ * a_z = 4 pi \rho_0 r_s^3 ( 1/((r+rs)*r^2) - log(1+r/rs)/r^3) * z
+ *      - dphi_PSC/dr*z/r
+ *
+ *
+ * @param time The current time.
+ * @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_acceleration(
+    double time, const struct external_potential* restrict potential,
+    const struct phys_const* restrict phys_const, struct gpart* restrict g) {
+
+#ifdef HAVE_LIBGSL
+
+  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];
+
+  /* First for the NFW part */
+  const float R2 = dx * dx + dy * dy;
+  const float r = sqrtf(R2 + dz * dz + potential->eps * potential->eps);
+
+  const float r_inv = 1.0f / r;
+  const float M_NFW = potential->pre_factor * (logf(1.0f + r / potential->r_s) -
+                                               r / (r + potential->r_s));
+  const float dpot_dr_NFW = M_NFW * r_inv * r_inv;
+  const float pot_nfw =
+      -potential->pre_factor * logf(1.0f + r / potential->r_s) * r_inv;
+  g->a_grav[0] -= potential->f[0] * dpot_dr_NFW * dx * r_inv;
+  g->a_grav[1] -= potential->f[0] * dpot_dr_NFW * dy * r_inv;
+  g->a_grav[2] -= potential->f[0] * dpot_dr_NFW * dz * r_inv;
+  gravity_add_comoving_potential(g, potential->f[0] * pot_nfw);
+
+  /* Now the the MN disk */
+  const float f1 = sqrtf(potential->Zdisk * potential->Zdisk + dz * dz);
+  const float f2 = potential->Rdisk + f1;
+  const float f3 = powf(R2 + f2 * f2, -1.5f);
+  const float mn_term = potential->Rdisk + sqrtf(potential->Zdisk + dz * dz);
+  const float pot_mn = -potential->Mdisk / sqrtf(R2 + mn_term * mn_term);
+
+  g->a_grav[0] -= potential->f[1] * potential->Mdisk * f3 * dx;
+  g->a_grav[1] -= potential->f[1] * potential->Mdisk * f3 * dy;
+  g->a_grav[2] -= potential->f[1] * potential->Mdisk * f3 * (f2 / f1) * dz;
+  gravity_add_comoving_potential(g, potential->f[1] * pot_mn);
+
+  /* Now the the PSC bulge */
+  const float r2 = r * r;
+  const float M_psc =
+      potential->prefactor_psc_1 *
+      (potential->gamma_psc -
+       gsl_sf_gamma_inc(1.5f - 0.5f * potential->alpha,
+                        r2 / (potential->r_c * potential->r_c)));
+  const float dpot_dr = M_psc / r2;
+  const float pot_psc =
+      -M_psc / r - potential->prefactor_psc_2 *
+                       gsl_sf_gamma_inc(1.0f - 0.5f * potential->alpha,
+                                        r2 / (potential->r_c * potential->r_c));
+
+  g->a_grav[0] -= potential->f[2] * dpot_dr * dx * r_inv;
+  g->a_grav[1] -= potential->f[2] * dpot_dr * dy * r_inv;
+  g->a_grav[2] -= potential->f[2] * dpot_dr * dz * r_inv;
+  gravity_add_comoving_potential(g, potential->f[2] * pot_psc);
+
+#else
+  error("Code not compiled with GSL. Can't compute MWPotential2014.");
+#endif
+}
+
+/**
+ * @brief Computes the gravitational potential energy of a particle in an
+ * NFW potential + MN potential.
+ *
+ * phi = f[0] * (-4 * pi * G * rho_0 * r_s^3 * ln(1+r/r_s)) - f[1] * (G * Mdisk
+ * / sqrt(R^2 + (Rdisk + sqrt(z^2 + Zdisk^2))^2)) + f[2] * [- G / r * (2 * pi *
+ * amplitude * r_1^alpha r_c^(3 - alpha) * gamma_inf((3 - alpha)/2, r^2 / r_c^2)
+ * ) - 2 * pi * G * amplitude * r_1^alpha * r_c^(2 - alpha)
+ * Gamma_sup((2-alpha)/2, r^2 / r_c^2 ) ]
+ *
+ * @param time The current time (unused here).
+ * @param potential The #external_potential used in the run.
+ * @param phys_const Physical constants in internal units.
+ * @param g Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float
+external_gravity_get_potential_energy(
+    double time, const struct external_potential* potential,
+    const struct phys_const* const phys_const, const struct gpart* g) {
+
+#ifdef HAVE_LIBGSL
+
+  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];
+
+  /* First for the NFW profile */
+  const float R2 = dx * dx + dy * dy;
+  const float r = sqrtf(R2 + dz * dz + potential->eps * potential->eps);
+  const float r_inv = 1.0f / r;
+  const float pot_nfw =
+      -potential->pre_factor * logf(1.0f + r / potential->r_s) * r_inv;
+
+  /* Now for the MN disk */
+  const float sqrt_term = sqrtf(dz * dz + potential->Zdisk * potential->Zdisk);
+  const float MN_denominator =
+      sqrtf(R2 + powf(potential->Rdisk + sqrt_term, 2.0f));
+  const float mn_pot = -potential->Mdisk / MN_denominator;
+
+  /* Now for PSC bulge */
+  const float r2 = r * r;
+  const float M_psc =
+      potential->prefactor_psc_1 *
+      (potential->gamma_psc -
+       gsl_sf_gamma_inc(1.5f - 0.5f * potential->alpha,
+                        r2 / (potential->r_c * potential->r_c)));
+  const float psc_pot =
+      -M_psc / r - potential->prefactor_psc_2 *
+                       gsl_sf_gamma_inc(1.0f - 0.5f * potential->alpha,
+                                        r2 / (potential->r_c * potential->r_c));
+
+  return phys_const->const_newton_G *
+         (potential->f[0] * pot_nfw + potential->f[1] * mn_pot +
+          potential->f[2] * psc_pot);
+
+#else
+  error("Code not compiled with GSL. Can't compute MWPotential2014.");
+  return 0.0;
+#endif
+}
+
+/**
+ * @brief Initialises the external potential properties in the internal system
+ * of units.
+ *
+ * @param parameter_file The parsed parameter file
+ * @param phys_const Physical constants in internal units
+ * @param us The current internal system of units
+ * @param potential The external potential properties to initialize
+ */
+static INLINE void potential_init_backend(
+    struct swift_params* parameter_file, const struct phys_const* phys_const,
+    const struct unit_system* us, const struct space* s,
+    struct external_potential* potential) {
+
+#ifdef HAVE_LIBGSL
+
+  /* Read in the position of the centre of potential */
+  parser_get_param_double_array(
+      parameter_file, "MWPotential2014Potential:position", 3, potential->x);
+
+  /* Is the position absolute or relative to the centre of the box? */
+  const int useabspos = parser_get_param_int(
+      parameter_file, "MWPotential2014Potential:useabspos");
+
+  /* Define the default value in the above system of units*/
+  const double c_200_default = 9.823403437774843;
+  const double M_200_Msun_default = 147.41031542774076e10; /* M_sun  */
+  const double H_default = 127.78254614201471e-2;          /* no unit  */
+  const double Mdisk_Msun_default = 6.8e10;                /* M_sun  */
+  const double Rdisk_kpc_default = 3.0;                    /* kpc  */
+  const double Zdisk_kpc_default = 0.280;                  /* kpc  */
+  const double amplitude_default = 1.0;                    /* no unit  */
+  const double r_1_kpc_default = 1.0;                      /* kpc  */
+  const double alpha_default = 1.8;                        /* no unit  */
+  const double r_c_kpc_default = 1.9;                      /* kpc  */
+  potential->f[0] = 0.4367419745056084;                    /* no unit  */
+  potential->f[1] = 1.002641971008805;                     /* no unit */
+  potential->f[2] = 0.022264787598364262;                  /* no unit */
+
+  if (!useabspos) {
+    potential->x[0] += s->dim[0] / 2.;
+    potential->x[1] += s->dim[1] / 2.;
+    potential->x[2] += s->dim[2] / 2.;
+  }
+
+  /* Read the other parameters of the model */
+  potential->timestep_mult = parser_get_param_double(
+      parameter_file, "MWPotential2014Potential:timestep_mult");
+
+  /* Bug fix : Read the softening length from the params file */
+  potential->eps = parser_get_param_double(parameter_file,
+                                           "MWPotential2014Potential:epsilon");
+
+  potential->c_200 = parser_get_opt_param_double(
+      parameter_file, "MWPotential2014Potential:concentration", c_200_default);
+  potential->M_200 = parser_get_opt_param_double(
+      parameter_file, "MWPotential2014Potential:M_200_Msun",
+      M_200_Msun_default);
+  potential->H = parser_get_opt_param_double(
+      parameter_file, "MWPotential2014Potential:H", H_default);
+  potential->Mdisk = parser_get_opt_param_double(
+      parameter_file, "MWPotential2014Potential:Mdisk_kpc", Mdisk_Msun_default);
+  potential->Rdisk = parser_get_opt_param_double(
+      parameter_file, "MWPotential2014Potential:Rdisk_kpc", Rdisk_kpc_default);
+  potential->Zdisk = parser_get_opt_param_double(
+      parameter_file, "MWPotential2014Potential:Zdisk_kpc", Zdisk_kpc_default);
+  potential->amplitude = parser_get_opt_param_double(
+      parameter_file, "MWPotential2014Potential:amplitude", amplitude_default);
+  potential->r_1 = parser_get_opt_param_double(
+      parameter_file, "MWPotential2014Potential:r_1_kpc", r_1_kpc_default);
+  potential->alpha = parser_get_opt_param_double(
+      parameter_file, "MWPotential2014Potential:alpha", alpha_default);
+  potential->r_c = parser_get_opt_param_double(
+      parameter_file, "MWPotential2014Potential:r_c_kpc", r_c_kpc_default);
+  parser_get_opt_param_double_array(
+      parameter_file, "MWPotential2014Potential:potential_factors", 3,
+      potential->f);
+
+  /* Convert to internal system of units by using the
+   * physical constants defined in this system */
+  potential->M_200 *= phys_const->const_solar_mass;
+  potential->H *= phys_const->const_reduced_hubble;
+  potential->Mdisk *= phys_const->const_solar_mass;
+  potential->Rdisk *= 1000. * phys_const->const_parsec;
+  potential->Zdisk *= 1000. * phys_const->const_parsec;
+  potential->r_1 *= 1000. * phys_const->const_parsec;
+  potential->r_c *= 1000. * phys_const->const_parsec;
+
+  /* Compute rho_c */
+  const double rho_c = 3.0 * potential->H * potential->H /
+                       (8.0 * M_PI * phys_const->const_newton_G);
+
+  /* Compute R_200 */
+  const double R_200 =
+      cbrtf(3.0 * potential->M_200 / (4. * M_PI * 200.0 * rho_c));
+
+  /* NFW scale-radius */
+  potential->r_s = R_200 / potential->c_200;
+  const double r_s3 = potential->r_s * potential->r_s * potential->r_s;
+
+  /* Log(c_200) term appearing in many expressions */
+  potential->log_c200_term =
+      log(1. + potential->c_200) - potential->c_200 / (1. + potential->c_200);
+
+  const double rho_0 =
+      potential->M_200 / (4.f * M_PI * r_s3 * potential->log_c200_term);
+
+  /* Pre-factor for the accelerations (note G is multiplied in later on) */
+  potential->pre_factor = 4.0f * M_PI * rho_0 * r_s3;
+
+  /* Prefactor for the mass of the PSC profile */
+  potential->prefactor_psc_1 = 2.0 * M_PI * potential->amplitude *
+                               pow(potential->r_1, potential->alpha) *
+                               pow(potential->r_c, 3.0 - potential->alpha);
+
+  /* Gamma function value for the mass of the PSC profile */
+  potential->gamma_psc = gsl_sf_gamma(1.5 - 0.5 * potential->alpha);
+
+  /* Prefactor for the potential of the PSC profile */
+  potential->prefactor_psc_2 = 2.0 * M_PI * potential->amplitude *
+                               pow(potential->r_1, potential->alpha) *
+                               pow(potential->r_c, 2.0 - potential->alpha);
+
+  /* Compute the orbital time at the softening radius */
+  const double sqrtgm = sqrt(phys_const->const_newton_G * potential->M_200);
+  const double epslnthing = log(1.0 + potential->eps / potential->r_s) -
+                            potential->eps / (potential->eps + potential->r_s);
+
+  potential->mintime = 2. * M_PI * potential->eps * sqrtf(potential->eps) *
+                       sqrtf(potential->log_c200_term / epslnthing) / sqrtgm *
+                       potential->timestep_mult;
+#else
+  error("Code not compiled with GSL. Can't compute MWPotential2014.");
+#endif
+}
+
+/**
+ * @brief Prints the properties of the external potential to stdout.
+ *
+ * @param  potential The external potential properties.
+ */
+static INLINE void potential_print_backend(
+    const struct external_potential* potential) {
+
+  message(
+      "External potential is 'MWPotential2014' "
+      "with properties (in internal units) are "
+      "(x,y,z) = "
+      "(%e, %e, %e), c_200 = %e, M_200 = %e, H = %e, M_disk = %e, R_disk = %e, "
+      "z_disk = %e, amplitude = %e, r_1 = %e, alpha = %e, r_c = %e, timestep "
+      "multiplier = %e mintime = %e",
+      potential->x[0], potential->x[1], potential->x[2], potential->c_200,
+      potential->M_200, potential->H, potential->Mdisk, potential->Rdisk,
+      potential->Zdisk, potential->amplitude, potential->r_1, potential->alpha,
+      potential->r_c, potential->timestep_mult, potential->mintime);
+}
+
+#endif /* SWIFT_POTENTIAL_MWPotential2014_H */