calc_mass_potential.py: add gravitational pontential calculation document

1 files changed, 40 insertions, 1 deletions

diff --git a/calc_mass_potential.py b/calc_mass_potential.py
index 4197d51..c1300b9 100755
--- a/calc_mass_potential.py
+++ b/calc_mass_potential.py
@@ -2,9 +2,11 @@
 #
 # Weitian LI
 # Created: 2016-06-24
-# Updated: 2016-06-25
+# Updated: 2016-06-26
 #
 # Change logs:
+# 2016-06-26:
+#   * Add document on gravitational potential calculation
 # 2016-06-25:
 #   * Rename method 'interpolate()' to 'fit_spline()'
 #   * Use 'InterpolatedUnivariateSpline' instead of 'interp1d'
@@ -60,10 +62,41 @@ For example:
     (1.0 keV) * (1.0 kpc) / (0.6 * m_atom * G) ~= 3.7379e10 [ Msun ]
 which is consistent with the formula of (ref.[2], eq.(3))
 
+------------------------------------------------------------
+
+Gravitational potential profile calculation:
+
+Newton's theorems (ref.[3], Sec. 2.2.1):
+1. A body that is inside a spherical shell of matter experiences no net
+   gravitational force from that shell.
+2. The gravitational force on a body that lies outside a spherical shell
+   of matter is the same as it would be if all the shell's matter were
+   concentrated into a point at its center.
+
+Therefore, the gravitational potential produced by a spherical shell of
+mass 'M' is:
+    phi = (1) - G * M / R;  (r <= R, i.e., inside the shell)
+          (2) - G * M / r;  (r > R, i.e., outside the shell)
+
+The total gravitational potential may be considered to be the sum of the
+potentials of spherical shells of mass
+    dM(r) = 4 * pi * rho(r) * r^2 dr,
+so we may calculate the gravitational potential at 'R' generated by an
+arbitrary spherically symmetric density distribution 'rho(r)' by adding
+the contributions to the potential produced by shells
+    (1) with r < R,
+and
+    (2) with r > R.
+In this way, we obtain
+    phi(R) = - (G/R) * \int_0^R dM(r) - G * \int_R^{\inf} dM(r)/r
+           = - 4*pi*G * ((1/R) * \int_0^R r^2 * rho(r) dr +
+             \int_R^{\inf} r * rho(r) dr)
 
+------------------------------------------------------------
 References:
 [1] Ettori et al., 2013, Space Science Review, 177, 119-154
 [2] Walker et al., 2012, MNRAS, 422, 3503
+[3] Tremaine & Binney, Galactic Dynamics, 2nd edition, 2008
 
 
 Sample configuration file:
@@ -350,6 +383,12 @@ class DensityProfile:
         self.m_total = m_total
         return m_total
 
+    def calc_potential(self, verbose=True):
+        """
+        Calculate the gravitational potential profile.
+        """
+        pass
+
     def plot(self, profile, ax=None, fig=None):
         pass