Merge remote-tracking branch 'origin/master' into parmetis-perm

Conflicts:
	configure.ac

.gitignore

 *~
+*.hdf5
 
 Makefile
 Makefile.in
@@ -23,14 +24,12 @@ doc/Doxyfile
 examples/swift
 examples/swift_mpi
 examples/*/*.xmf
-examples/*/*.hdf5
 examples/*/*.h5
 examples/*/*.png
 examples/*/*.txt
 examples/*/*.dot
 examples/*/used_parameters.yml
 examples/*/*/*.xmf
-examples/*/*/*.hdf5
 examples/*/*/*.png
 examples/*/*/*.txt
 examples/*/*/*.dot
@@ -73,7 +72,6 @@ tests/test_nonsym_force_1_vec.dat
 tests/test_nonsym_force_2_vec.dat
 tests/testGreetings
 tests/testReading
-tests/input.hdf5
 tests/testSingle
 tests/testTimeIntegration
 tests/testSPHStep
@@ -107,6 +105,8 @@ tests/testDump
 tests/testLogger
 tests/benchmarkInteractions
 tests/testGravityDerivatives
+tests/testPotentialSelf
+tests/testPotentialPair
 
 theory/latex/swift.pdf
 theory/SPH/Kernels/kernels.pdf
@@ -277,9 +277,14 @@ _minted*
 *.sympy
 sympy-plots-for-*.tex/
 
+# python
+*.pyc
+
 # todonotes
 *.tdo
 
 # xindy
 *.xdy
 
+# macOS
+*.DS_Store
\ No newline at end of file

configure.ac

@@ -16,7 +16,7 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 # Init the project.
-AC_INIT([SWIFT],[0.6.0],[https://gitlab.cosma.dur.ac.uk/swift/swiftsim])
+AC_INIT([SWIFT],[0.7.0],[https://gitlab.cosma.dur.ac.uk/swift/swiftsim])
 swift_config_flags="$*"
 
 #  Need to define this, instead of using fifth argument of AC_INIT, until 2.64.
@@ -407,13 +407,47 @@ AC_HEADER_STDC
 # Check for the libraries we will need.
 AC_CHECK_LIB(m,sqrt,,AC_MSG_ERROR(something is wrong with the math library!))
 
-# Check for GSL
+# Check for GSL. We test for this in the standard directories by default,
+# and only disable if using --with-gsl=no or --without-gsl. When a value
+# is given GSL must be found.
 have_gsl="no"
-AC_CHECK_LIB([gslcblas], [cblas_dgemm])
-AC_CHECK_LIB([gsl], [gsl_integration_qag])
-if test "x$ac_cv_lib_gslcblas_cblas_dgemm" != "x"; then
-   have_gsl="yes"
+AC_ARG_WITH([gsl],
+    [AS_HELP_STRING([--with-gsl=PATH],
+       [root directory where GSL is installed @<:@yes/no@:>@]
+    )],
+    [with_gsl="$withval"],
+    [with_gsl="test"]
+)
+if test "x$with_gsl" != "xno"; then
+   if test "x$with_gsl" != "xyes" -a "x$with_gsl" != "xtest" -a "x$with_gsl" != "x"; then
+      GSL_LIBS="-L$with_gsl/lib -lgsl -lgslcblas"
+      GSL_INCS="-I$with_gsl/include"
+   else
+      GSL_LIBS="-lgsl -lgslcblas"
+      GSL_INCS=""
+   fi
+   #  GSL is not specified, so just check if we have it.
+   if test "x$with_gsl" = "xtest"; then
+      AC_CHECK_LIB([gslcblas],[cblas_dgemm],[have_gsl="yes"],[have_gsl="no"],$GSL_LIBS)
+      if test "x$have_gsl" != "xno"; then
+         AC_DEFINE([HAVE_LIBGSLCBLAS],1,[The GSL CBLAS library appears to be present.])
+         AC_CHECK_LIB([gsl],[gsl_integration_qag],
+            AC_DEFINE([HAVE_LIBGSL],1,[The GSL library appears to be present.]),
+            [have_gsl="no"],$GSL_LIBS)
+      fi
+   else
+      AC_CHECK_LIB([gslcblas],[cblas_dgemm],
+         AC_DEFINE([HAVE_LIBGSLCBLAS],1,[The GSL CBLAS library appears to be present.]),
+         AC_MSG_ERROR(something is wrong with the GSL CBLAS library!), $GSL_LIBS)
+      AC_CHECK_LIB([gsl],[gsl_integration_qag],
+         AC_DEFINE([HAVE_LIBGSL],1,[The GSL library appears to be present.]),
+         AC_MSG_ERROR(something is wrong with the GSL library!), $GSL_LIBS)
+      have_gsl="yes"
+   fi
 fi
+AC_SUBST([GSL_LIBS])
+AC_SUBST([GSL_INCS])
+AM_CONDITIONAL([HAVEGSL],[test -n "$GSL_LIBS"])
 
 # Check for pthreads.
 AX_PTHREAD([LIBS="$PTHREAD_LIBS $LIBS" CFLAGS="$CFLAGS $PTHREAD_CFLAGS"
@@ -443,7 +477,7 @@ AC_ARG_WITH([parmetis],
     [AS_HELP_STRING([--with-parmetis=PATH],
        [root directory where ParMETIS is installed @<:@yes/no@:>@]
     )],
-    [],
+    [with_parmetis="$withval"],
     [with_parmetis="no"]
 )
 
@@ -506,13 +540,13 @@ AH_VERBATIM([__STDC_FORMAT_MACROS],
 #define __STDC_FORMAT_MACROS 1
 #endif])
 
-# Check for grackle. 
+# Check for grackle.
 have_grackle="no"
 AC_ARG_WITH([grackle],
     [AS_HELP_STRING([--with-grackle=PATH],
        [root directory where grackle is installed @<:@yes/no@:>@]
     )],
-    [],
+    [with_grackle="$withval"],
     [with_grackle="no"]
 )
 if test "x$with_grackle" != "xno"; then
@@ -1013,7 +1047,7 @@ esac
 #  External potential
 AC_ARG_WITH([ext-potential],
    [AS_HELP_STRING([--with-ext-potential=<pot>],
-      [external potential @<:@none, point-mass, isothermal, softened-isothermal, disc-patch, sine-wave default: none@:>@]
+      [external potential @<:@none, point-mass, point-mass-ring, point-mass-softened, isothermal, softened-isothermal, disc-patch, sine-wave, default: none@:>@]
    )],
    [with_potential="$withval"],
    [with_potential="none"]
@@ -1034,6 +1068,12 @@ case "$with_potential" in
    sine-wave)
       AC_DEFINE([EXTERNAL_POTENTIAL_SINE_WAVE], [1], [Sine wave external potential in 1D])
    ;;
+   point-mass-ring)
+      AC_DEFINE([EXTERNAL_POTENTIAL_POINTMASS_RING], [1], [Point mass potential for Keplerian Ring (Hopkins 2015).])
+   ;;
+   point-mass-softened)
+      AC_DEFINE([EXTERNAL_POTENTIAL_POINTMASS_SOFT], [1], [Softened point-mass potential with form 1/(r^2 + softening^2).])
+   ;;
    *)
       AC_MSG_ERROR([Unknown external potential: $with_potential])
    ;;

examples/EAGLE_100/eagle_100.yml

@@ -6,6 +6,15 @@ InternalUnitSystem:
   UnitCurrent_in_cgs:  1             # Amperes
   UnitTemp_in_cgs:     1             # Kelvin
 
+# Cosmological parameters
+Cosmology:
+  h:              0.6777        # Reduced Hubble constant
+  a_begin:        0.9090909     # Initial scale-factor of the simulation
+  a_end:          1.0           # Final scale factor of the simulation
+  Omega_m:        0.307         # Matter density parameter
+  Omega_lambda:   0.693         # Dark-energy density parameter
+  Omega_b:        0.0455        # Baryon density parameter
+  
 # Parameters governing the time integration
 TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).
@@ -19,7 +28,7 @@ Scheduler:
 # Parameters governing the snapshots
 Snapshots:
   basename:            eagle # Common part of the name of output files
-  time_first:          0.    # Time of the first output (in internal units)
+  time_first:          1.    # Time of the first output (in internal units)
   delta_time:          1e-3  # Time difference between consecutive outputs (in internal units)
 
 # Parameters governing the conserved quantities statistics
@@ -29,8 +38,8 @@ Statistics:
 # Parameters for the self-gravity scheme
 Gravity:
   eta:                   0.025    # Constant dimensionless multiplier for time integration. 
-  epsilon:               0.0001   # Softening length (in internal units).
-  theta:                 0.7      # Opening angle (Multipole acceptance criterion)
+  epsilon:               0.001   # Softening length (in internal units).
+  theta:                 0.85      # Opening angle (Multipole acceptance criterion)
   
 # Parameters for the hydrodynamics scheme
 SPH:

examples/EAGLE_12/eagle_12.yml

@@ -6,13 +6,6 @@ InternalUnitSystem:
   UnitCurrent_in_cgs:  1             # Amperes
   UnitTemp_in_cgs:     1             # Kelvin
 
-# Parameters governing the time integration
-TimeIntegration:
-  time_begin: 0.    # The starting time of the simulation (in internal units).
-  time_end:   1e-2  # 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:     1e-4  # The maximal time-step size of the simulation (in internal units).
-
 # Cosmological parameters
 Cosmology:
   h:              0.6777        # Reduced Hubble constant
@@ -21,6 +14,13 @@ Cosmology:
   Omega_m:        0.307         # Matter density parameter
   Omega_lambda:   0.693         # Dark-energy density parameter
   Omega_b:        0.0455        # Baryon density parameter
+
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   1e-2  # 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:     1e-4  # The maximal time-step size of the simulation (in internal units).
   
 Scheduler:
   max_top_level_cells:    15
@@ -30,6 +30,7 @@ Snapshots:
   basename:            eagle # Common part of the name of output files
   time_first:          1.    # Time of the first output (in internal units)
   delta_time:          1e-3  # Time difference between consecutive outputs (in internal units)
+  compression:         4
 
 # Parameters governing the conserved quantities statistics
 Statistics:

examples/EAGLE_25/eagle_25.yml

@@ -6,6 +6,15 @@ InternalUnitSystem:
   UnitCurrent_in_cgs:  1             # Amperes
   UnitTemp_in_cgs:     1             # Kelvin
 
+# Cosmological parameters
+Cosmology:
+  h:              0.6777        # Reduced Hubble constant
+  a_begin:        0.9090909     # Initial scale-factor of the simulation
+  a_end:          1.0           # Final scale factor of the simulation
+  Omega_m:        0.307         # Matter density parameter
+  Omega_lambda:   0.693         # Dark-energy density parameter
+  Omega_b:        0.0455        # Baryon density parameter
+  
 # Parameters governing the time integration
 TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).
@@ -19,8 +28,9 @@ Scheduler:
 # Parameters governing the snapshots
 Snapshots:
   basename:            eagle # Common part of the name of output files
-  time_first:          0.    # Time of the first output (in internal units)
+  time_first:          1.    # Time of the first output (in internal units)
   delta_time:          1e-3  # Time difference between consecutive outputs (in internal units)
+  compression:         4
 
 # Parameters governing the conserved quantities statistics
 Statistics:

examples/EAGLE_50/eagle_50.yml

@@ -6,6 +6,15 @@ InternalUnitSystem:
   UnitCurrent_in_cgs:  1             # Amperes
   UnitTemp_in_cgs:     1             # Kelvin
 
+# Cosmological parameters
+Cosmology:
+  h:              0.6777        # Reduced Hubble constant
+  a_begin:        0.9090909     # Initial scale-factor of the simulation
+  a_end:          1.0           # Final scale factor of the simulation
+  Omega_m:        0.307         # Matter density parameter
+  Omega_lambda:   0.693         # Dark-energy density parameter
+  Omega_b:        0.0455        # Baryon density parameter
+  
 # Parameters governing the time integration
 TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).
@@ -19,7 +28,7 @@ Scheduler:
 # Parameters governing the snapshots
 Snapshots:
   basename:            eagle # Common part of the name of output files
-  time_first:          0.    # Time of the first output (in internal units)
+  time_first:          1.    # Time of the first output (in internal units)
   delta_time:          1e-3  # Time difference between consecutive outputs (in internal units)
 
 # Parameters governing the conserved quantities statistics

examples/EAGLE_6/eagle_6.yml

@@ -6,6 +6,15 @@ InternalUnitSystem:
   UnitCurrent_in_cgs:  1             # Amperes
   UnitTemp_in_cgs:     1             # Kelvin
 
+# Cosmological parameters
+Cosmology:
+  h:              0.6777        # Reduced Hubble constant
+  a_begin:        0.9090909     # Initial scale-factor of the simulation
+  a_end:          1.0           # Final scale factor of the simulation
+  Omega_m:        0.307         # Matter density parameter
+  Omega_lambda:   0.693         # Dark-energy density parameter
+  Omega_b:        0.0455        # Baryon density parameter
+  
 Scheduler:
   max_top_level_cells: 8
   
@@ -15,12 +24,13 @@ TimeIntegration:
   time_end:   1e-2  # 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:     1e-4  # The maximal time-step size of the simulation (in internal units).
-
+  
 # Parameters governing the snapshots
 Snapshots:
   basename:            eagle # Common part of the name of output files
   time_first:          0.    # Time of the first output (in internal units)
   delta_time:          1e-3  # Time difference between consecutive outputs (in internal units)
+  compression:         4
 
 # Parameters governing the conserved quantities statistics
 Statistics:

examples/EvrardCollapse_3D/evrard.yml

@@ -18,7 +18,8 @@ Snapshots:
   basename:            evrard # Common part of the name of output files
   time_first:          0.       # Time of the first output (in internal units)
   delta_time:          0.1     # Time difference between consecutive outputs (in internal units)
-
+  compression:         4
+  
 # Parameters governing the conserved quantities statistics
 Statistics:
   delta_time:          1e-2 # Time between statistics output

examples/GreshoVortex_3D/gresho.yml

@@ -21,7 +21,8 @@ Snapshots:
   basename:            gresho # Common part of the name of output files
   time_first:          0.     # Time of the first output (in internal units)
   delta_time:          1e-1   # Time difference between consecutive outputs (in internal units)
-
+  compression:         4
+  
 # Parameters governing the conserved quantities statistics
 Statistics:
   delta_time:          1e-2 # Time between statistics output

examples/KelvinHelmholtzGrowthRate_3D/kelvinHelmholtzGrowthRate.yml

@@ -18,7 +18,8 @@ Snapshots:
   basename:            kelvinHelmholtzGrowthRate  # Common part of the name of output files
   time_first:          0.               # Time of the first output (in internal units)
   delta_time:          0.04      # Time difference between consecutive outputs (in internal units)
-
+  compression:         4
+  
 # Parameters governing the conserved quantities statistics
 Statistics:
   delta_time:          1e-2 # Time between statistics output

examples/KeplerianRing/README.md

+Keplerian Ring Test
+===================
+
+This test uses the '2D Keplerian Ring' to assess the accuracy of SPH
+calculations. For more information on the test, please see Cullen & Dehnen 2009
+(http://arxiv.org/abs/1006.1524) Section 4.3.
+
+It is key that for this test in particular that the initial conditions are
+perfectly arranged (here we use the same methodology as described in the paper
+referenced above).
+
+The test uses:
+
++ $M = 10^10$ for a central point mass
++ $r$ up to 5 (particles created to edge of box when relevant)
++ Approximately 16000 particles
++ $c_s = 0.01 \ll v_\phi = 10$, enforced by giving them a mass of 1 unit and an
+  internal energy of 0.015.
+
+Please note that the initial condition generator **requires python3 rather than
+python2**, as well as the plotting script.
+
+
+Code Setup
+----------
+
+You will need to recompile SWIFT with an external potential. This can be done
+by using
+
+    ./configure --with-ext-potential=point-mass
+
+in the root directory of the project. We suggest leaving all other code options
+as their defaults.
+
+Please also consider using:
+
+    ./configure --with-ext-potential=point-mass-softened
+
+if you are running with the initial conditions generated with the script and
+using a nonzero softening.
+
+
+Initial Conditions Generation
+-----------------------------
+
+The initial condition generation is handled by ```makeIC.py```, for which
+the options are available with
+
+    python3 makeIC.py --help
+
+however the defaults are very sensible. If you wish to change the central mass,
+you will need to change the external point mass in ```keplerian_ring.yml``` as
+well.
+
+There are three main types of generation, handled through `--generationmethod`:
+
++ `spiral`: using the original spiral method with a density distribution
+            created using different mass particles.
++ `gaussian`: a guassian density profile with density distribution changed
+              through the _distribution_ of equal mass particles.
++ `grid`: a grid with a flat density profile (similar to Hopkins 2013) imposed
+          on it through different particle masses.
+
+All of them have their benefits and drawbacks, so give each a go.
+
+Plotting
+--------
+
+Once you have ran swift (we suggest that you use the following)
+
+    ../swift -g -S -s -t 16 keplerian_ring.yml 2>&1 | tee output.log
+
+there will be around 350 ```.hdf5``` files in your directory. To check out
+the results of the example use the plotting script:
+
+    python3 plotSolution.py
+
+which will produce a png file in the directory. Check out `plotSolution.py
+--help` for more options.
+
+You can also try the `make_movie.py` script which should make a movie.
+
+
+Theory
+------
+
+We use the 'spiral' technique to generate the initial conditions. To do this,
+we first calculate
+$$
+    m_i = \frac{i - \frac{1}{2}}{N}
+$$
+for each particle $i$ with $N$ the total number of required particles. Then the
+radius of each particle is given by
+$$
+    r_i = f^{-1}(m_i)
+$$
+where for us the PDF $f$ is a Gaussian.
+
+To choose the radial positions, we then choose an initial angle $\theta_1$, and
+from there
+$$
+    \delta \theta_i = \left[2\pi\left( 1 - \frac{r_{i-1}}{r_i} \right)\right]~.
+$$
+
+### Inverse CDF of Gaussian
+
+We need:
+$$
+    f^{-1}(m_i) = r_i = r_{central} + \sqrt{2}\sigma \mathrm{erf}^{-1}(2m_i - 1),
+$$
+and thankfully scipy has an inverse error function built in.
+
+
+Implementation Details
+----------------------
+
++ The $m_i$ are generated by ```generate_m_i```.
++ The inverse Gaussian PDF is implemented as ```inverse_gaussian```.
++ The $\theta_i$ are generated by ```generate_theta_i```.
+
+
+Useful Information
+------------------
+
+$G$ in simulation units (by default using this yaml) is given by:
+$$
+    4.300970\times10^{-6}
+$$
+This is used to convert the mass of the central potential (used in the
+parameterfile) to the $GM$ required in the initial conditions generator.

examples/KeplerianRing/keplerian_ring.yml

+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1   # Grams
+  UnitLength_in_cgs:   1   # Centimeters
+  UnitVelocity_in_cgs: 1   # Centimeters per second
+  UnitCurrent_in_cgs:  1   # Amperes
+  UnitTemp_in_cgs:     1   # Kelvin
+
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   50    # The end time of the simulation (in internal units).
+  dt_min:     1e-6  # 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:            keplerian_ring # Common part of the name of output files
+  time_first:          0.        # Time of the first output (in internal units)
+  delta_time:          0.2      # Time difference between consecutive outputs (in internal units)
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          1e-2 # Time between statistics output
+
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+  delta_neighbours:      0.1      # The tolerance for the targetted number of neighbours.
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  initial_conditions.hdf5        # The file to read
+  shift_x:    0.                   # A shift to apply to all particles read from the ICs (in internal units).
+  shift_y:    0.
+  shift_z:    0.
+
+# External potential parameters
+PointMassPotential:
+  position_x:      5.     # location of external point mass in internal units
+  position_y:      5.     # here just take this to be the centre of the ring
+  position_z:      5.	
+  mass:            1.498351e7     # mass of external point mass in internal units
+  timestep_mult:   0.03     # controls time step
+  softening:       0.05

examples/KeplerianRing/make_movie.py

+"""
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2017
+#
+# Josh Borrow (joshua.borrow@durham.ac.uk)
+#
+# 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/>.
+#
+#
+# -----------------------------------------------------------------------------
+#
+# This program creates test plots for the initial condition generator provided
+# for the Keplerian Ring example.
+#
+###############################################################################
+"""
+
+import matplotlib
+matplotlib.use("Agg")
+
+
+import matplotlib.pyplot as plt
+import matplotlib.animation as anim
+import numpy as np
+import h5py as h5
+
+from tqdm import tqdm
+
+
+def get_colours(numbers, norm=None, cm_name="viridis"):
+    if norm is None:
+        norm = matplotlib.colors.Normalize(vmax=numbers.max(), vmin=numbers.min())
+
+    cmap = matplotlib.cm.get_cmap(cm_name)
+
+    return cmap(norm(numbers)), norm
+        
+
+def load_data(filename, silent=True, extra=None, logextra=True, exclude=None):
+    if not silent:
+        print(f"Loading data from {filename}")
+
+    with h5.File(filename, "r") as file_handle:
+        coords = file_handle['PartType0']['Coordinates'][...]
+        time = float(file_handle['Header'].attrs['Time'])
+        boxsize = file_handle['Header'].attrs['BoxSize']
+
+        # For some old runs we have z=0
+        if np.sum(coords[:, 2]) == 0:
+            centre_of_box = np.array(list(boxsize[:2]/2) + [0])
+        else:
+            centre_of_box = boxsize/2
+
+        
+        if exclude is not None:
+            distance_from_centre = coords - centre_of_box
+            r2 = np.sum(distance_from_centre * distance_from_centre, 1)
+            mask = r2 < (exclude * exclude)
+
+            masked = np.ma.array(coords, mask=np.array([mask]*3))
+
+            coords = masked.compressed()
+        
+
+        if extra is not None:
+            other = file_handle['PartType0'][extra][...]
+            if exclude is not None:
+                masked = np.ma.array(other, mask=mask)
+
+                other = masked.compressed()
+
+            if logextra:
+                other = np.log(other)
+
+            return coords, time, other
+        else:
+            return coords, time
+
+
+def rms(x):
+    return np.sqrt(sum(x**2))
+
+
+def rotation_velocity_at_r(r, params):
+    """
+    Gets the rotation velocity at a given radius r by assuming it is keplerian.
+
+    Assumes we are in cgs units, which may one day be our downfall.
+    """
+
+    unit_length = float(params[r"InternalUnitSystem:UnitLength_in_cgs"])
+
+    if unit_length != 1.:
+        print(f"Your unit length: {unit_length}")
+        raise InternalUnitSystemError(
+            "This function is only created to handle CGS units."
+        )
+
+    central_mass = float(params["PointMassPotential:mass"])