calc_mass.py = calc_mass_potential.py - potential calculation

1 files changed, 6 insertions, 96 deletions

diff --git a/calc_mass_potential.py b/calc_mass.py
index bb49501..d283c7e 100755
--- a/calc_mass_potential.py
+++ b/calc_mass.py
@@ -2,9 +2,12 @@
 #
 # Aaron LI
 # Created: 2016-06-24
-# Updated: 2016-07-11
+# Updated: 2016-07-13
 #
 # Change logs:
+# 2016-07-13:
+#   * Rename from 'calc_mass_potential.py' to 'calc_mass.py'
+#   * Split out the potential profile calculation -> 'calc_potential.py'
 # 2016-07-11:
 #   * Use a default config to allow a minimal user config
 # 2016-07-10:
@@ -80,39 +83,9 @@ For example:
 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
 """
 
 
@@ -135,7 +108,7 @@ from spline import SmoothSpline
 plt.style.use("ggplot")
 
 config_default = """
-## Configuration for `calc_mass_potential.py`
+## Configuration for `calc_mass.py`
 
 # electron density profile
 ne_profile = ne_profile.txt
@@ -210,8 +183,6 @@ class DensityProfile:
     m_total = None
     # total mass density profile (required by potential calculation)
     rho_total = None
-    # potential profile (cut at the largest available radius)
-    potential = None
     # fitted spline to the profiles
     d_spline = None
     ne_spline = None
@@ -482,46 +453,6 @@ class DensityProfile:
         self.rho_total = rho_total
         return rho_total
 
-    def calc_potential(self, verbose=True):
-        """
-        Calculate the gravitational potential profile.
-
-        NOTE:
-        The integral in the potential formula can only be integrated
-        to the largest radius of availability.
-        This limitation will affect the absolute values of the calculated
-        potentials, but not the shape of the potential profile.
-        """
-        def _int_inner(x):
-            return x**2 * self.eval_spline(spline="rho_total", x=x)
-
-        def _int_outer(x):
-            return x * self.eval_spline(spline="rho_total", x=x)
-
-        if self.rho_total is None:
-            self.calc_density_total(verbose=verbose)
-        if self.rho_total_spline is None:
-            self.fit_spline(spline="rho_total", log10=["x", "y"])
-        potential = np.zeros(self.radius.shape)
-        if verbose:
-            print("Calculating the potential profile (#%d): ..." %
-                  len(self.radius), end="", flush=True)
-        r_max = max(self.radius)
-        c = - 4 * np.pi * au.kpc.to(au.cm)**2
-        for i, r in enumerate(self.radius):
-            if verbose and (i+1) % 10 == 0:
-                print("%d..." % (i+1), end="", flush=True)
-            # total gravitational potential [ cm^2/s^2 ]
-            potential[i] = c * ac.G.cgs.value * (
-                (1/r) * integrate.quad(_int_inner, a=0.0, b=r,
-                                       epsrel=1e-2)[0] +
-                integrate.quad(_int_outer, a=r, b=r_max,
-                               epsrel=1e-2)[0])
-        if verbose:
-            print("DONE!", flush=True)
-        self.potential = potential
-        return potential
-
     def plot(self, profile, with_spline=True, ax=None, fig=None):
         x = self.radius
         xlabel = "3D Radius"
@@ -549,11 +480,6 @@ class DensityProfile:
             y = self.rho_total
             ylabel = "Total mass density"
             yunit = "g/cm$^3$"
-        elif profile == "potential":
-            y = self.potential
-            ylabel = "Gravitational potential"
-            yunit = "cm$^2$/s$^2$"
-            yscale = "linear"
         else:
             raise ValueError("unknown profile name: %s" % profile)
         #
@@ -590,11 +516,6 @@ class DensityProfile:
                                     self.radius_err,
                                     self.rho_total])
             header = "radius[kpc]  radius_err[kpc]  density_total[g/cm^3]"
-        elif profile == "potential":
-            data = np.column_stack([self.radius,
-                                    self.radius_err,
-                                    self.potential])
-            header = "radius[kpc]  radius_err[kpc]  potential[???]"
         else:
             raise ValueError("unknown profile name: %s" % profile)
         np.savetxt(outfile, data, header=header)
@@ -602,7 +523,7 @@ class DensityProfile:
 
 def main():
     parser = argparse.ArgumentParser(
-            description="Calculate the mass and potential profiles")
+            description="Calculate the gas/total mass profiles")
     parser.add_argument("config", nargs="?",
                         help="additional config")
     args = parser.parse_args()
@@ -651,17 +572,6 @@ def main():
     density_profile.plot(profile="rho_total", ax=ax, fig=fig)
     fig.savefig(config["rho_total_profile_image"], dpi=150)
 
-    outfile_potential = config["potential_profile"]
-    if outfile_potential != "":
-        density_profile.calc_potential(verbose=True)
-        density_profile.save(profile="potential",
-                             outfile=outfile_potential)
-        fig = Figure(figsize=(10, 8))
-        FigureCanvas(fig)
-        ax = fig.add_subplot(1, 1, 1)
-        density_profile.plot(profile="potential", ax=ax, fig=fig)
-        fig.savefig(config["potential_profile_image"], dpi=150)
-
 
 if __name__ == "__main__":
     main()