diff options
| -rwxr-xr-x | calc_potential.py | 234 |
1 files changed, 234 insertions, 0 deletions
diff --git a/calc_potential.py b/calc_potential.py new file mode 100755 index 0000000..f2e3b66 --- /dev/null +++ b/calc_potential.py @@ -0,0 +1,234 @@ +#!/usr/bin/env python3 +# +# Aaron LI +# Created: 2016-07-13 +# Updated: 2016-07-13 +# +# Change logs: +# + +""" +Calculate the gravitational potential profile through the +total mass density profile. + +Newton's theorems (ref.[1], 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] Tremaine & Binney, Galactic Dynamics, 2nd edition, 2008 +""" + + +import argparse + +import numpy as np +import astropy.units as au +import astropy.constants as ac +import scipy.integrate as integrate +from configobj import ConfigObj +import matplotlib.pyplot as plt +from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas +from matplotlib.figure import Figure + +from spline import SmoothSpline + +plt.style.use("ggplot") + + +config_default = """ +## Configuration for `calc_potential.py` + +# total mass density profile +rho_total_profile = rho_total_profile.txt + +# output gravitational potential profile +potential_profile = potential_profile.txt +potential_profile_image = potential_profile.png +""" + + +class DensityProfile: + """ + Cluster's gas/total mass density profile. + """ + # supported types of density profile + DENSITY_TYPES = ["total", "gas"] + # available splines + SPLINES = ["rho_total"] + # input density data: [r, r_err, m] + r = None + r_err = None + d = None + # redshift of the object + redshift = None + # fitted smoothing spline to total density profile + rho_total_spline = None + # calculated potential profile + potential = None + + def __init__(self, data, density_type="total"): + self.load_data(data=data, density_type=density_type) + + def load_data(self, data, density_type="total"): + if density_type not in self.DENSITY_TYPES: + raise ValueError("invalid density_types: %s" % density_type) + # 3-column density profile: r[kpc], r_err[kpc], rho[g/cm^3] + self.r = data[:, 0].copy() + self.r_err = data[:, 1].copy() + self.d = data[:, 2].copy() + self.density_type = density_type + + def calc_potential(self, verbose=True): + """ + Calculate the gravitational potential profile from the + total mass density 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.density_type != "total": + raise ValueError("total mass density profile is required") + + if self.rho_total_spline is None: + self.fit_spline(spline="rho_total", log10=["x", "y"]) + potential = np.zeros(self.r.shape) + if verbose: + print("Calculating the potential profile (#%d): ..." % + len(self.r), end="", flush=True) + r_max = max(self.r) + c = - 4 * np.pi * au.kpc.to(au.cm)**2 + for i, r in enumerate(self.r): + 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, ax=None, fig=None): + x = self.r + xlabel = "3D Radius" + xunit = "kpc" + xscale = "log" + yscale = "log" + y = self.potential + ylabel = "Gravitational potential" + yunit = "cm$^2$/s$^2$" + yscale = "linear" + # + if ax is None: + fig, ax = plt.subplots(1, 1) + ax.plot(x, y, linewidth=2) + ax.set_xscale(xscale) + ax.set_yscale(yscale) + ax.set_xlim(min(x), max(x)) + y_min, y_max = min(y), max(y) + y_min = y_min/1.2 if y_min > 0 else y_min*1.2 + y_max = y_max*1.2 if y_max > 0 else y_max/1.2 + ax.set_ylim(y_min, y_max) + ax.set_xlabel(r"%s (%s)" % (xlabel, xunit)) + ax.set_ylabel(r"%s (%s)" % (ylabel, yunit)) + fig.tight_layout() + return (fig, ax) + + def save(self, outfile): + """ + Save calculated overdensity profile. + """ + data = np.column_stack([self.r, + self.r_err, + self.potential]) + header = "radius[kpc] radius_err[kpc] potential[cm^2/s^2]" + np.savetxt(outfile, data, header=header) + + def fit_spline(self, spline, log10=[]): + if spline not in self.SPLINES: + raise ValueError("invalid spline: %s" % spline) + # + if spline == "rho_total": + # input total mass density profile + x = self.r + y = self.d + spl = "rho_total_spline" + else: + raise ValueError("invalid spline: %s" % spline) + setattr(self, spl, SmoothSpline(x=x, y=y)) + getattr(self, spl).fit(log10=log10) + + def eval_spline(self, spline, x): + """ + Evaluate the specified spline at the supplied positions. + Also check whether the spline was fitted in the log-scale space, + and transform the evaluated values if necessary. + """ + if spline == "rho_total": + spl = self.rho_total_spline + else: + raise ValueError("invalid spline: %s" % spline) + # + return spl.eval(x) + + +def main(): + parser = argparse.ArgumentParser( + description="Calculate the gravitational potential profile") + parser.add_argument("config", nargs="?", + help="additional config") + args = parser.parse_args() + + config = ConfigObj(config_default.splitlines()) + if args.config is not None: + config_user = ConfigObj(open(args.config)) + config.merge(config_user) + + density_data = np.loadtxt(config["rho_total_profile"]) + density_profile = DensityProfile(data=density_data, density_type="total") + density_profile.calc_potential(verbose=True) + density_profile.save(outfile=config["potential_profile"]) + fig = Figure(figsize=(10, 8)) + FigureCanvas(fig) + ax = fig.add_subplot(1, 1, 1) + density_profile.plot(ax=ax, fig=fig) + fig.savefig(config["potential_profile_image"], dpi=150) + + +if __name__ == "__main__": + main() |