diff options
| -rwxr-xr-x | deproject_sbp.py | 71 |
1 files changed, 62 insertions, 9 deletions
diff --git a/deproject_sbp.py b/deproject_sbp.py index acae430..8f28062 100755 --- a/deproject_sbp.py +++ b/deproject_sbp.py @@ -6,6 +6,7 @@ # # Change logs: # 2016-06-22: +# * Add class 'DensityProfile', the inversion to 'BrightnessProfile' # * Add classes 'AstroParams' and 'BrightnessProfile' # * Add class 'ChandraPixel' # * Update documentation @@ -160,7 +161,7 @@ class AstroParams: # molecular weight per electron (0.3 solar abundance; grsa table) [1] mu_e = 1.155 # atomic mass unit - m_atom = 1.660539040e-24 # [ g ] + m_atom = au.u.to(au.g) # [ g ] class Projection: @@ -1013,14 +1014,11 @@ class BrightnessProfile: * The radii should have unit [ pixel ] (Chandra pixel) * The brightness should have unit [ photon s^-1 cm^-2 pixel^-2 ], i.e., [ Flux pixel^-2 ] (radial profile column `SUR_FLUX`) - - Reference: ref. [4], eq.(9) below """ - # input SBP data: [r, r_err, s, s_err] + # input SBP data: [r, r_err, s] r = None r_err = None s = None - s_err = None # redshift of the source z = None # `ChandraPixel` instance for unit conversion @@ -1035,11 +1033,10 @@ class BrightnessProfile: self.pixel = ChandraPixel(z) def load_data(self, data): - # 4-column SBP + # 3-column SBP: [r, r_err, brightness] self.r = data[:, 0].copy() self.r_err = data[:, 1].copy() self.s = data[:, 2].copy() - self.s_err = data[:, 3].copy() def load_cf_data(self, data): # 2-column cooling function profile @@ -1058,11 +1055,10 @@ class BrightnessProfile: self.r_err *= cm_per_pixel self.cf_radius *= cm_per_pixel self.s /= cm_per_pixel**2 - self.s_err /= cm_per_pixel**2 self.units_converted = True def get_radius(self): - return self.r.copy() + return (self.r.copy(), self.r_err.copy()) def calc_electron_density(self): """ @@ -1094,6 +1090,63 @@ class BrightnessProfile: rho = ne * AstroParams.mu_e * AstroParams.m_atom return rho + +class DensityProfile: + """ + Calculate the 2D surface brightness profile by projecting the + electron number density profile / gas mass density profile and + incorporating the cooling function profile. + + NOTE: + * The radii 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 + + def __init__(self, cf_data): + self.load_cf_data(data=cf_data) + + 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_data(self, data): + # 2-column cooling function profile + self.cf_radius = data[:, 0].copy() + self.cf_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.density_type == "gas": + # convert gas mass to electron number density + ne = self.d / AstroParams.m_atom / AstroParams.mu_e + else: + ne = self.d + # + cf_interp = scipy.interpolate.interp1d(x=self.cf_radius, + y=self.cf_value, + kind="linear", + assume_sorted=True) + # flux per unit volume + flux = 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 + ###################################################################### |