diff options
| -rwxr-xr-x | calc_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() |