#!/usr/bin/env python3 # # Copyright (c) 2016 Aaron LI # MIT license # # Created: 2016-07-04 # # Change logs: # 2016-07-13: # * Error on missing user config file # * Force the input temperature data to be 2D # 2016-07-11: # * Use a default config to allow a minimal user config # 2016-07-04: # * Set default "rcut=3000" for TemperatureProfile.extrapolate() # """ Fit the deprojected ICM temperature data points with a self-proposed temperature profile model, i.e., the *wang2012* model: T(r) = A * (pow(x,n)+xi*a2) / (pow(x,n)+a2) / pow(1+x*x/a3/a3, beta) + T0 With the fitted temperature profile model, we can interpolate and extrapolate the temperature profile for later mass profile calculation. """ import argparse import json from collections import OrderedDict import numpy as np import astropy.units as au import lmfit from configobj import ConfigObj import matplotlib.pyplot as plt from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas from matplotlib.figure import Figure from fitting_models import FittingModel from astro_params import ChandraPixel plt.style.use("ggplot") config_default = """ ## Configuration for `fit_tprofile.py` # redshift of the object (for pixel to distance conversion) redshift = # input temperature profile data file t_profile_data = t_profile_data.txt # cut radius to which stop the extrapolation (unit: kpc) rcut_extrap = 3000 # number of data points for the output temperature profile num_dp = 1000 # output json file to save the fitting results t_profile_json = t_profile.json # output interpolated and extrapolated temperature profile t_profile = t_profile.txt t_profile_image = t_profile.png [model_params] # name = initial, lower, upper, variable (FIXED/False to fix the parameter) A = 5.0, 1.0, 500 n = 5.0, 0.1, 10 xi = 0.3, 0.1, 1.0 a2 = 2000, 1.0, 1e+05 a3 = 1000, 100, 3000 #beta = 0.5, 0.1, 1.0, FIXED beta = 0.5, 0.1, 1.0 T0 = 2.0, 1.0, 5.0 """ class Wang2012Model(FittingModel): """ *wang2012* model proposed to fit the ICM temperature profile. """ name = "Wang2012 Temperature Profile Model" # model parameters params = lmfit.Parameters() params.add_many( # (name, value, vary, min, max, expr) ("A", 5.0, True, 1.0, 500, None), ("n", 5.0, True, 0.1, 10, None), ("xi", 0.3, True, 0.1, 1.0, None), ("a2", 2000, True, 1.0, 1.0e5, None), ("a3", 1000, True, 100, 3000, None), ("beta", 0.5, True, 0.1, 1.0, None), ("T0", 2.0, True, 1.0, 5.0, None)) def __init__(self, fit_method="lbfgsb", params=None): super().__init__(fit_method=fit_method, params=params, scale=False) @staticmethod def model(x, params): parvals = params.valuesdict() A = parvals["A"] n = parvals["n"] xi = parvals["xi"] a2 = parvals["a2"] a3 = parvals["a3"] beta = parvals["beta"] T0 = parvals["T0"] return (A * (x**n + xi*a2) / (x**n + a2) / ((1 + (x/a3)**2) ** beta) + T0) class TemperatureProfile: """ Fit the deprojected ICM temperature data points with a temperature profile model, and output the interpolated and extrapolated temperature profile for later mass profile calculation. The input radii have unit "pixel", which are first converted to "kpc" and then fitted with the model. The output temperature profile also has unit "kpc" for radii. """ # input temperature profile data: [r, r_err, t, t_err] r = None r_err = None t = None t_err = None # redshift of the source z = None # `ChandraPixel` instance for unit conversion pixel = None # flag to indicate whether the units are converted units_converted = False # model to be fitted model = None def __init__(self, data, z): self.load_data(data) self.z = z self.pixel = ChandraPixel(z) self.model = Wang2012Model() def load_data(self, data): # 4-column t profile: [r, r_err, temperature, temperature_err] self.r = data[:, 0].copy() self.r_err = data[:, 1].copy() self.t = data[:, 2].copy() self.t_err = data[:, 3].copy() def convert_units(self): """ Convert the units of input data: radius: pixel -> kpc """ if not self.units_converted: kpc_per_pixel = self.pixel.get_length().to(au.kpc).value self.r *= kpc_per_pixel self.r_err *= kpc_per_pixel self.units_converted = True def fit(self): self.model.load_data(xdata=self.r, ydata=self.t, xerr=self.r_err, yerr=self.t_err) self.model.fit() def extrapolate(self, rcut=3000, num=1000): """ Interpolate and extrapolate the fitted temperature profile. The output radii are generated to be linear-evenly distributed. """ self.rcut_extrap = rcut self.num_dp = num radius = np.linspace(0.0, rcut, num+1) rin = radius[:-1] rout = radius[1:] self.r_extrapolated = (rout + rin) / 2.0 self.r_err_extrapolated = (rout - rin) / 2.0 self.t_extrapolated = self.model.f(self.r_extrapolated) def report(self, outfile=None): """ Report the temperature profile model fitting results. """ results = OrderedDict([ ("redshift", self.z), ("rcut_extrap", self.rcut_extrap), ("num_dp", self.num_dp), ("model", self.model.name), ("fitting", self.model.report(rtype="fitting")), ("params", self.model.report(rtype="parameters")), ]) results_json = json.dumps(results, indent=2) if outfile is None: print(results_json) else: open(outfile, "w").write(results_json+"\n") def save(self, outfile): data = np.column_stack([self.r_extrapolated, self.r_err_extrapolated, self.t_extrapolated]) header = "radius[kpc] radius_err[kpc] temperature[keV]" np.savetxt(outfile, data, header=header) def plot(self, ax=None, fig=None): if ax is None: fig, ax = plt.subplots(1, 1) ax.errorbar(self.r, self.t, xerr=self.r_err, yerr=self.t_err, fmt="none", elinewidth=2, capthick=2) # fitted model ax.plot(self.r_extrapolated, self.t_extrapolated, color="black", linestyle="solid", linewidth=2) ax.set_xlabel("Radius (kpc)") ax.set_ylabel("Temperature (keV)") fig.tight_layout() return (fig, ax) def main(): parser = argparse.ArgumentParser( description="temperature profile fit, interpolate and extrapolate") parser.add_argument("config", nargs="?", default="tprofile.conf", help="configuration (default: tprofile.conf)") args = parser.parse_args() config = ConfigObj(config_default.splitlines()) config_user = ConfigObj(open(args.config)) config.merge(config_user) tprofile_data = np.loadtxt(config["t_profile_data"], ndmin=2) redshift = config.as_float("redshift") tprofile = TemperatureProfile(tprofile_data, redshift) tprofile.convert_units() # Load parameters settings from config params = config["model_params"] for p, value in params.items(): variable = True if len(value) == 4 and value[3].upper() in ["FIXED", "FALSE"]: variable = False tprofile.model.set_param(name=p, value=float(value[0]), min=float(value[1]), max=float(value[2]), vary=variable) tprofile.fit() tprofile.extrapolate(rcut=config.as_float("rcut_extrap"), num=config.as_int("num_dp")) tprofile.report(outfile=config["t_profile_json"]) tprofile.save(outfile=config["t_profile"]) fig = Figure(figsize=(10, 8)) FigureCanvas(fig) ax = fig.add_subplot(1, 1, 1) tprofile.plot(ax=ax, fig=fig) fig.savefig(config["t_profile_image"], dpi=150) if __name__ == "__main__": main()