diff options
| -rwxr-xr-x | calc_mass_potential.py | 129 |
1 files changed, 59 insertions, 70 deletions
diff --git a/calc_mass_potential.py b/calc_mass_potential.py index b038050..8122aef 100755 --- a/calc_mass_potential.py +++ b/calc_mass_potential.py @@ -2,9 +2,12 @@ # # Aaron LI # Created: 2016-06-24 -# Updated: 2016-06-29 +# Updated: 2016-07-04 # # Change logs: +# 2016-07-04: +# * Remove unnecessary configuration options +# * Update radii unit to be "kpc", mass unit to be "Msun" # 2016-06-29: # * Update "plot()" to support electron number density profile # 2016-06-28: @@ -112,21 +115,14 @@ Sample configuration file: ## Configuration for `calc_mass_potential.py` ## Date: 2016-06-28 -# 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 profiles (interpolation) num_dp = 1000 @@ -180,8 +176,7 @@ class DensityProfile: * gravitational potential profile (cut at the largest available radius) NOTE: - * The radii (of density profile and cooling function profile) - should have unit [ cm ] + * The radii should have unit [ kpc ] * The density should have unit [ cm^-3 ] or [ g cm^-3 ] """ # available splines @@ -239,21 +234,35 @@ class DensityProfile: 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] + # 3-column density profile: r[kpc], r_err[kpc], 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() + if data.shape[1] == 2: + # 2-column cooling function profile: r[kpc], cf[flux/EM] + self.cf_radius = data[:, 0].copy() + self.cf_value = data[:, 1].copy() + elif data.shape[1] == 3: + # 3-column cooling function profile: r[kpc], r_err, cf[flux/EM] + self.cf_radius = data[:, 0].copy() + self.cf_value = data[:, 2].copy() + else: + raise ValueError("invalid cooling function profile data") 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() + if data.shape[1] == 2: + # 2-column temperature profile: r[kpc], T[keV] + self.t_radius = data[:, 0].copy() + self.t_value = data[:, 1].copy() + elif data.shape[1] == 3: + # 3-column temperature profile: r[kpc], r_err[kpc], T[keV] + self.t_radius = data[:, 0].copy() + self.t_value = data[:, 2].copy() + else: + raise ValueError("invalid temperature profile data") def calc_density_electron(self): """ @@ -295,8 +304,9 @@ class DensityProfile: # flux per unit volume cf_new = self.eval_spline(spline="cooling_function", x=self.r) flux = cf_new * ne ** 2 / AstroParams.ratio_ne_np - # project the 3D flux - projector = Projection(rout=self.r+self.r_err) + # project the 3D flux into 2D brightness + rout = (self.r + self.r_err) * au.kpc.to(au.cm) + projector = Projection(rout) brightness = projector.project(flux) return brightness @@ -503,13 +513,13 @@ class DensityProfile: if verbose: print("Calculating the gas mass profile (#%d): ..." % len(self.radius), end="", flush=True) + c = au.kpc.to(au.cm)**3 * au.g.to(au.solMass) for i, r in enumerate(self.radius): if verbose and (i+1) % 50 == 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=1e-5*au.kpc.to(au.cm), - epsrel=1e-3)[0] + # enclosed gas mass [ Msun ] + m_gas[i] = c * integrate.quad(_f_rho_gas, a=0.0, b=r, + epsrel=1e-3)[0] if verbose: print("DONE!", flush=True) self.m_gas = m_gas @@ -529,8 +539,8 @@ class DensityProfile: # # calculate the coefficient of the total mass formula # M0 = (k_B * T_gas(r) * r) / (mu * m_atom * G) - M0 = ((1.0*au.keV) * (1.0*au.cm) / - (AstroParams.mu * ac.u * ac.G)).to(au.g).value + M0 = ((1.0*au.keV) * (1.0*au.kpc) / + (AstroParams.mu * ac.u * ac.G)).to(au.solMass).value m_total = np.zeros(self.radius.shape) if verbose: print("Calculating the total mass profile (#%d): ..." % @@ -540,11 +550,11 @@ class DensityProfile: print("%d..." % (i+1), end="", flush=True) T = self.eval_spline(spline="temperature", x=r) dT_dr = derivative(lambda r: self.eval_spline("temperature", r), - r, dx=0.01*au.kpc.to(au.cm)) + r, dx=0.01) ne = self.eval_spline(spline="electron", x=r) dne_dr = derivative(lambda r: self.eval_spline("electron", r), - r, dx=0.01*au.kpc.to(au.cm)) - # enclosed total mass [ g ] + r, dx=0.01) + # enclosed total mass [ Msun ] m_total[i] = - M0 * T * r * (((r / T) * dT_dr) + ((r / ne) * dne_dr)) if verbose: @@ -563,10 +573,12 @@ class DensityProfile: if verbose: print("Calculating the total mass density profile ...") rho_total = np.zeros(self.radius.shape) + # unit conversion: Msun/kpc^3 -> g/cm^3 + c = au.solMass.to(au.g) / au.kpc.to(au.cm)**3 for i, r in enumerate(self.radius): dM_dr = derivative(lambda r: self.eval_spline("mass_total", r), - r, dx=0.01*au.kpc.to(au.cm)) - rho_total[i] = (dM_dr / (4 * np.pi * r**2)) + r, dx=0.01) + rho_total[i] = (dM_dr / (4 * np.pi * r**2)) * c self.rho_total = rho_total return rho_total @@ -595,15 +607,15 @@ class DensityProfile: 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) - potential[i] = - 4 * np.pi * ac.G.cgs.value * ( + # 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, - epsabs=1e-3*au.kpc.to(au.cm), epsrel=1e-2)[0] + integrate.quad(_int_outer, a=r, b=r_max, - epsabs=1e-3*au.kpc.to(au.cm), epsrel=1e-2)[0]) if verbose: print("DONE!", flush=True) @@ -611,26 +623,26 @@ class DensityProfile: return potential def plot(self, profile, with_spline=True, ax=None, fig=None): - x = self.radius * au.cm.to(au.kpc) + x = self.radius xlabel = "3D Radius" xunit = "kpc" xscale = "log" yscale = "log" x_spl, y_spl = None, None if profile == "electron": - x = self.r * au.cm.to(au.kpc) + x = self.r y = self.ne ylabel = "Electron number density" yunit = "cm$^{-3}$" if with_spline: - x_spl = self.radius * au.cm.to(au.kpc) + x_spl = self.radius y_spl = self.eval_spline(spline="electron", x=self.radius) elif profile == "mass_gas": - y = self.m_gas * au.g.to(au.solMass) + y = self.m_gas ylabel = "Gas mass" yunit = "M$_{\odot}$" elif profile == "mass_total": - y = self.m_total * au.g.to(au.solMass) + y = self.m_total ylabel = "Total mass" yunit = "M$_{\odot}$" elif profile == "rho_total": @@ -640,7 +652,7 @@ class DensityProfile: elif profile == "potential": y = self.potential ylabel = "Gravitational potential" - yunit = "???" + yunit = "cm$^2$/s$^2$" yscale = "linear" else: raise ValueError("unknown profile name: %s" % profile) @@ -653,7 +665,10 @@ class DensityProfile: ax.set_xscale(xscale) ax.set_yscale(yscale) ax.set_xlim(min(x), max(x)) - ax.set_ylim(min(y)/1.2, max(y)*1.2) + 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() @@ -664,22 +679,22 @@ class DensityProfile: data = np.column_stack([self.radius, self.radius_err, self.m_gas]) - header = "radius[cm] radius_err[cm] mass_gas(<r)[g]" + header = "radius[kpc] radius_err[kpc] mass_gas(<r)[Msun]" elif profile == "mass_total": data = np.column_stack([self.radius, self.radius_err, self.m_total]) - header = "radius[cm] radius_err[cm] mass_total(<r)[g]" + header = "radius[kpc] radius_err[kpc] mass_total(<r)[Msun]" 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]" + 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[cm] radius_err[cm] potential[???]" + header = "radius[kpc] radius_err[kpc] potential[???]" else: raise ValueError("unknown profile name: %s" % profile) np.savetxt(outfile, data, header=header) @@ -692,37 +707,11 @@ def main(): 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") |