From a5303c868402bec6b1a4b46dccb4c7fef4c829d7 Mon Sep 17 00:00:00 2001
From: Aaron LI <aaronly.me@outlook.com>
Date: Mon, 27 Jun 2016 11:09:43 +0800
Subject: calc_mass_potential.py: implement gravitational pontential
 calculation

---
 calc_mass_potential.py | 97 ++++++++++++++++++++++++++++++++++++++++++++------
 1 file changed, 87 insertions(+), 10 deletions(-)

diff --git a/calc_mass_potential.py b/calc_mass_potential.py
index c1300b9..eca425c 100755
--- a/calc_mass_potential.py
+++ b/calc_mass_potential.py
@@ -2,9 +2,12 @@
 #
 # Weitian LI
 # Created: 2016-06-24
-# Updated: 2016-06-26
+# Updated: 2016-06-27
 #
 # Change logs:
+# 2016-06-27:
+#   * Implement potential profile calculation:
+#     methods 'calc_density_total()' and 'calc_potential()'
 # 2016-06-26:
 #   * Add document on gravitational potential calculation
 # 2016-06-25:
@@ -129,6 +132,9 @@ m_gas_profile = mass_gas_profile.txt
 # output total (gravitational) mass profile
 m_total_profile = mass_total_profile.txt
 
+# output total mass density profile
+rho_total_profile = rho_total_profile.txt
+
 # output gravitational potential profile
 potential_profile = potential_profile.txt
 ------------------------------------------------------------
@@ -189,6 +195,8 @@ class DensityProfile:
     m_gas = None
     # total (gravitational) mass profile: M_total(<r)
     m_total = None
+    # total mass density profile (required by potential calculation)
+    rho_total = None
     # potential profile (cut at the largest available radius)
     potential = None
 
@@ -221,7 +229,7 @@ class DensityProfile:
         """
         if self.cf_radius is None or self.cf_value is None:
             raise ValueError("cooling function profile missing")
-        ne = self.calc_electron_density()
+        ne = self.calc_density_electron()
         # flux per unit volume
         flux = self.cf_spline(self.r) * ne ** 2 / AstroParams.ratio_ne_np
         # project the 3D flux
@@ -229,7 +237,7 @@ class DensityProfile:
         brightness = projector.project(flux)
         return brightness
 
-    def calc_electron_density(self):
+    def calc_density_electron(self):
         """
         Calculate the electron number density from the gas mass density
         if necessary.
@@ -240,7 +248,7 @@ class DensityProfile:
             self.ne = self.d / AstroParams.m_atom / AstroParams.mu_e
         return self.ne
 
-    def calc_gas_density(self):
+    def calc_density_gas(self):
         """
         Calculate the gas mass density from the electron number density
         if necessary.
@@ -337,8 +345,8 @@ class DensityProfile:
         #
         m_gas = np.zeros(self.radius.shape)
         if verbose:
-            print("Calculating the gas mass profile (#%d): ..." % len(m_gas),
-                  end="", flush=True)
+            print("Calculating the gas mass profile (#%d): ..." %
+                  len(self.radius), end="", flush=True)
         for i, r in enumerate(self.radius):
             if verbose and (i+1) % 50 == 0:
                 print("%d..." % (i+1), end="", flush=True)
@@ -368,7 +376,7 @@ class DensityProfile:
         m_total = np.zeros(self.radius.shape)
         if verbose:
             print("Calculating the total mass profile (#%d): ..." %
-                  len(m_total), end="", flush=True)
+                  len(self.radius), end="", flush=True)
         for i, r in enumerate(self.radius):
             if verbose and (i+1) % 100 == 0:
                 print("%d..." % (i+1), end="", flush=True)
@@ -383,11 +391,64 @@ class DensityProfile:
         self.m_total = m_total
         return m_total
 
+    def calc_density_total(self, verbose=True):
+        """
+        Calculate the total mass density profile, which is required to
+        calculate the following gravitational potential profile.
+        """
+        rho_total = np.zeros(self.radius.shape)
+        if verbose:
+            print("Calculating the total mass density profile ...")
+            print(">>> Fitting spline to total mass profile ...")
+        self.m_total_spline = interpolate.InterpolatedUnivariateSpline(
+            x=self.radius, y=self.m_total, ext="const")
+        for i, r in enumerate(self.radius):
+            rho_total[i] = (derivative(self.m_total_spline, r,
+                                       dx=0.01*au.kpc.to(au.cm)) /
+                            (4 * np.pi * r**2))
+        self.rho_total = rho_total
+        if verbose:
+            print(">>> Fitting spline to total mass density profile ...")
+        self.rho_total_spline = interpolate.InterpolatedUnivariateSpline(
+            x=self.radius, y=rho_total, ext="const")
+
     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.
         """
-        pass
+        def _int_inner(x):
+            return x**2 * self.rho_total_spline(x)
+
+        def _int_outer(x):
+            return x * self.rho_total_spline(x)
+
+        if self.rho_total is None:
+            self.calc_density_total(verbose=verbose)
+        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)
+        for i, r in enumerate(self.radius):
+            if verbose and (i+1) % 50 == 0:
+                print("%d..." % (i+1), end="", flush=True)
+            potential[i] = - 4 * np.pi * ac.G.cgs.value * (
+                (1/r) * integrate.quad(_int_inner, a=0.0, b=r,
+                                       epsabs=1e-5*au.kpc.to(au.cm),
+                                       epsrel=1e-3)[0] +
+                integrate.quad(_int_outer, a=r, b=r_max,
+                               epsabs=1e-5*au.kpc.to(au.cm),
+                               epsrel=1e-3)[0])
+        if verbose:
+            print("DONE!", flush=True)
+        self.potential = potential
+        return potential
 
     def plot(self, profile, ax=None, fig=None):
         pass
@@ -403,6 +464,16 @@ class DensityProfile:
                                     self.radius_err,
                                     self.m_total])
             header = "radius[cm]  radius_err[cm]  mass_total(<r)[g]"
+        elif profile == "rho_total":
+            data = np.column_stack([self.radius,
+                                    self.radius_err,
+                                    self.rho_total])
+            header = "radius[cm]  radius_err[cm]  density_total[g/cm^3]"
+        elif profile == "potential":
+            data = np.column_stack([self.radius,
+                                    self.radius_err,
+                                    self.potential])
+            header = "radius[cm]  radius_err[cm]  potential[???]"
         else:
             raise ValueError("unknown profile name: %s" % profile)
         np.savetxt(outfile, data, header=header)
@@ -451,8 +522,8 @@ def main():
                                      density_type="electron")
     density_profile.load_cf_profile(cf_profile)
     density_profile.load_t_profile(t_profile)
-    density_profile.calc_electron_density()
-    density_profile.calc_gas_density()
+    density_profile.calc_density_electron()
+    density_profile.calc_density_gas()
     density_profile.fit_spline(verbose=True)
     density_profile.gen_radius(num=config.as_int("num_dp"))
     density_profile.calc_mass_gas(verbose=True)
@@ -461,6 +532,12 @@ def main():
     density_profile.calc_mass_total(verbose=True)
     density_profile.save(profile="mass_total",
                          outfile=config["m_total_profile"])
+    density_profile.calc_density_total(verbose=True)
+    density_profile.save(profile="rho_total",
+                         outfile=config["rho_total_profile"])
+    density_profile.calc_potential(verbose=True)
+    density_profile.save(profile="potential",
+                         outfile=config["potential_profile"])
 
 
 if __name__ == "__main__":
-- 
cgit v1.2.2