diff options
| -rwxr-xr-x | calc_mass_potential.py | 329 |
1 files changed, 329 insertions, 0 deletions
diff --git a/calc_mass_potential.py b/calc_mass_potential.py new file mode 100755 index 0000000..22c3ccd --- /dev/null +++ b/calc_mass_potential.py @@ -0,0 +1,329 @@ +#!/usr/bin/env python3 +# +# Weitian LI +# Created: 2016-06-24 +# Updated: 2016-06-24 +# + +""" +Calculate the (gas and gravitational) mass profile and gravitational +potential profile from the electron number density profile. +The temperature profile is required. + +References: +[1] Ettori et al, 2013, Space Science Review, 177, 119-154 + + +Sample configuration file: +------------------------------------------------------------ +## Configuration for `calc_mass_potential.py` +## Date: 2016-06-24 + +# redshift used for pixel to distance conversion +redshift = <REDSHIFT> + +# electron density profile +ne_profile = ne_profile.txt + +# cooling function profile +cf_profile = coolfunc_profile.txt +# unit of the CF profile radius (default: pixel) +cf_unit = "pixel" + +# temperature profile +t_profile = t_profile.txt +# unit of the T profile radius (default: pixel) +t_unit = "pixel" + +# number of data points for the output profile calculation +num_dp = 100 + +# output gas mass profile +m_gas_profile = mass_gas_profile.txt + +------------------------------------------------------------ +""" + +import argparse + +import numpy as np +import astropy.units as au +import scipy.interpolate as interpolate +import scipy.integrate as integrate +from configobj import ConfigObj + +from astro_params import AstroParams, ChandraPixel +from projection import Projection + + +class DensityProfile: + """ + Utilize the 3D (electron number or gas mass) density profile to + calculate the following quantities: + * 2D projected surface brightness (requires cooling function profile) + * gas mass profile (integrated, M_gas(<r)) + * gravitational mass profile (M(<r); requires temperature profile) + * gravitational potential profile (cut at the largest available radius) + + NOTE: + * The radii (of density profile and cooling function profile) + should have unit [ cm ] + * The density should have unit [ cm^-3 ] or [ g cm^-3 ] + """ + # allowed density profile types + DENSITY_TYPES = ["electron", "gas"] + # input density data: [r, r_err, d] + r = None + r_err = None + d = None + # electron number density + ne = None + # gas mass density + rho_gas = None + # cooling function profile + cf_radius = None + cf_value = None + # temperature profile + t_radius = None + t_value = None + # interpolated profiles + d_interp = None + cf_interp = None + t_interp = None + # generated radial data points for profile calculation + radius = None + radius_err = None + # gas mass profile: M_gas(<r); same length as the above 'radius' + m_gas = None + # total (gravitational) mass profile: M_total(<r) + m_total = None + # potential profile (cut at the largest available radius) + potential = None + + def __init__(self, density, density_type="electron"): + self.load_data(data=density, density_type=density_type) + + def load_data(self, data, density_type="electron"): + if density_type not in self.DENSITY_TYPES: + raise ValueError("invalid density_types: %s" % density_type) + # 3-column density profile: [r, r_err, density] + self.r = data[:, 0].copy() + self.r_err = data[:, 1].copy() + self.d = data[:, 2].copy() + self.density_type = density_type + + def load_cf_profile(self, data): + # 2-column cooling function profile: r[cm], cf[flux/EM] + self.cf_radius = data[:, 0].copy() + self.cf_value = data[:, 1].copy() + + def load_t_profile(self, data): + # 2-column temperature profile: r[cm], T[keV] + self.t_radius = data[:, 0].copy() + self.t_value = data[:, 1].copy() + + def calc_brightness(self): + """ + Project the electron number density or gas mass density profile + to calculate the 2D surface brightness profile. + """ + if self.cf_radius is None or self.cf_value is None: + raise ValueError("cooling function profile missing") + ne = self.calc_electron_density() + # flux per unit volume + flux = self.cf_interp(self.r) * ne ** 2 / AstroParams.ratio_ne_np + # project the 3D flux + projector = Projection(rout=self.r+self.r_err) + brightness = projector.project(flux) + return brightness + + def calc_electron_density(self): + """ + Calculate the electron number density from the gas mass density + if necessary. + """ + if self.density_type == "electron": + self.ne = self.d.copy() + elif self.density_type == "gas": + self.ne = self.d / AstroParams.m_atom / AstroParams.mu_e + return self.ne + + def calc_gas_density(self): + """ + Calculate the gas mass density from the electron number density + if necessary. + """ + if self.density_type == "electron": + self.rho_gas = self.d * AstroParams.mu_e * AstroParams.m_atom + elif self.density_type == "gas": + self.rho_gas = self.d.copy() + return self.rho_gas + + def interpolate(self): + """ + Interpolate the density profile, cooling function profile, + and temperature profile. + + NOTE: + * Linear interpolation may be bad because the total mass calculation + needs to take the derivative of electron density profile and + temperature profile. Therefore, smooth spline is a better choice. + * Allow cooling function profile and temperature profile to be + extrapolated by filling with the last value if necessary. + + XXX: + * How to extrapolate the smooth spline if necessary ?? + """ + # density profile + # insert a data point at radius of zero + self.d_interp = interpolate.interp1d( + x=np.concatenate([[0.0], self.r]), + y=np.concatenate([[self.d[0]], self.d]), + kind="linear", + bounds_error=False, fill_value=self.d[-1], + assume_sorted=True) + if self.ne is not None: + self.ne_interp = interpolate.interp1d( + x=np.concatenate([[0.0], self.r]), + y=np.concatenate([[self.ne[0]], self.ne]), + kind="linear", + bounds_error=False, fill_value=self.ne[-1], + assume_sorted=True) + if self.rho_gas is not None: + self.rho_gas_interp = interpolate.interp1d( + x=np.concatenate([[0.0], self.r]), + y=np.concatenate([[self.rho_gas[0]], self.rho_gas]), + kind="linear", + bounds_error=False, fill_value=self.rho_gas[-1], + assume_sorted=True) + # cooling function profile + self.cf_interp = interpolate.interp1d( + x=self.cf_radius, y=self.cf_value, kind="linear", + bounds_error=False, fill_value=self.cf_value[-1], + assume_sorted=True) + # temperature profile + self.t_interp = interpolate.interp1d( + x=self.t_radius, y=self.t_value, kind="linear", + bounds_error=False, fill_value=self.t_value[-1], + assume_sorted=True) + + def gen_radius(self, num=100): + """ + Generate radial points for following mass and potential calculation. + + The generated radial points are logarithmic-evenly distributed. + """ + rout = self.r + self.r_err + rout_new = np.logspace(np.log10(rout[0]), np.log10(rout[-1]), + num=num, base=10.0) + rin_new = np.concatenate([[0.0], rout_new[:-1]]) + self.radius = (rout_new + rin_new) / 2.0 + self.radius_err = (rout_new - rin_new) / 2.0 + + def calc_mass_gas(self, verbose=False): + """ + Calculate the gas mass profile, i.e., the mass of the gas within + each radius. + + Reference: ref.[1], eq.(9) + """ + def _f_rho_gas(r): + return self.rho_gas_interp(r) * 4*np.pi * r**2 + # + m_gas = np.zeros(self.radius.shape) + if verbose: + print("Calculating the gas mass profile (#%d) ... " % len(m_gas)) + for i, r in enumerate(self.radius): + if verbose and (i+1) % 10 == 0: + print("%d..." % (i+1), end="", flush=True) + # integrated gas mass [ g ] + m_gas[i] = integrate.quad(_f_rho_gas, a=0.0, b=r, + epsabs=1.0e5, epsrel=1.e-2)[0] + if verbose: + print("DONE!", flush=True) + self.m_gas = m_gas + return m_gas + + def calc_mass_total(self, verbose=True): + """ + Calculate the total mass (i.e., gravitational mass) profile, + under the assumption of hydrostatic equilibrium (HE). + + References: ref.[1], eq.(5,6,7) + """ + if self.cf_radius is None or self.cf_value is None: + raise ValueError("cooling function profile required") + if self.t_radius is None or self.t_value is None: + raise ValueError("temperature profile required") + # + m_total = np.zeros(self.radius.shape) + if verbose: + print("Calculating the total mass profile (#%d) ... " % + len(m_total)) + # TODO: + self.m_total = m_total + return m_total + + def save(self, profile, outfile): + if profile == "mass_gas": + data = np.column_stack([self.radius, + self.radius_err, + self.m_gas]) + header = "radius[cm] radius_err[cm] mass_gas(<r)[g]" + else: + raise ValueError("unknown profile name: %s" % profile) + np.savetxt(outfile, data, header=header) + + +def main(): + parser = argparse.ArgumentParser( + description="Calculate the mass and potential profiles") + parser.add_argument("config", nargs="?", default="mass_potential.conf", + help="config for mass and potential profiles " + + "calculation (default: mass_potential.conf)") + args = parser.parse_args() + + config = ConfigObj(args.config) + + redshift = config.as_float("redshift") + pixel = ChandraPixel(z=redshift) + + ne_profile = np.loadtxt(config["ne_profile"]) + + cf_profile = np.loadtxt(config["cf_profile"]) + cf_unit = "pixel" + try: + cf_unit = config["cf_unit"] + except KeyError: + pass + if cf_unit == "pixel": + conv_factor = pixel.get_length().to(au.cm).value + else: + conv_factor = au.Unit(cf_unit).to(au.cm) + cf_profile[:, 0] *= conv_factor + + t_profile = np.loadtxt(config["t_profile"]) + t_unit = "pixel" + try: + t_unit = config["t_unit"] + except KeyError: + pass + if t_unit == "pixel": + conv_factor = pixel.get_length().to(au.cm).value + else: + conv_factor = au.Unit(t_unit).to(au.cm) + t_profile[:, 0] *= conv_factor + + density_profile = DensityProfile(density=ne_profile, + density_type="electron") + density_profile.load_cf_profile(cf_profile) + density_profile.load_t_profile(t_profile) + density_profile.calc_gas_density() + density_profile.interpolate() + density_profile.gen_radius(num=config.as_int("num_dp")) + density_profile.calc_mass_gas(verbose=True) + density_profile.save(profile="mass_gas", outfile=config["m_gas_profile"]) + + +if __name__ == "__main__": + main() |