summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rwxr-xr-xcalc_mass.py (renamed from calc_mass_potential.py)102
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()