diff options
Diffstat (limited to 'fg21sim/extragalactic')
| -rw-r--r-- | fg21sim/extragalactic/clusters/halo.py | 119 | ||||
| -rw-r--r-- | fg21sim/extragalactic/clusters/helper.py | 102 | ||||
| -rw-r--r-- | fg21sim/extragalactic/clusters/main.py | 7 |
3 files changed, 98 insertions, 130 deletions
diff --git a/fg21sim/extragalactic/clusters/halo.py b/fg21sim/extragalactic/clusters/halo.py index 0f8752a..beaff9e 100644 --- a/fg21sim/extragalactic/clusters/halo.py +++ b/fg21sim/extragalactic/clusters/halo.py @@ -36,6 +36,10 @@ References .. [sarazin1999] Sarazin 1999, ApJ, 520, 529 http://adsabs.harvard.edu/abs/1999ApJ...520..529S + +.. [hogg1999] + Hogg 1999, arXiv:astro-ph/9905116 + http://adsabs.harvard.edu/abs/1999astro.ph..5116H """ import logging @@ -48,6 +52,7 @@ from .emission import SynchrotronEmission from ...share import CONFIGS, COSMO from ...utils.units import (Units as AU, UnitConversions as AUC) +from ...utils.convert import Fnu_to_Tb logger = logging.getLogger(__name__) @@ -364,73 +369,137 @@ class RadioHalo: emissivity = syncem.emissivity(frequencies) return emissivity - def calc_power(self, emissivity): + def calc_power(self, frequencies, emissivity=None, **kwargs): """ Calculate the halo synchrotron power (i.e., power *emitted* per - unit frequency) from emissivity. + unit frequency) by assuming the emissivity is uniform throughout + the halo volume. + + NOTE + ---- + The calculated power (a.k.a. spectral luminosity) is in units of + [W/Hz] which is common in radio astronomy, instead of [erg/s/Hz]. + 1 [W] = 1e7 [erg/s] Parameters ---------- - emissivity : float, or 1D `~numpy.ndarray` - The synchrotron emissivity at multiple frequencies. + frequencies : float, or 1D `~numpy.ndarray` + The frequencies where to calculate the synchrotron power. + Unit: [MHz] + emissivity : float, or 1D `~numpy.ndarray`, optional + The synchrotron emissivity at the input frequencies. + If not provided, then invoke above ``calc_emissivity()`` + method to calculate them. Unit: [erg/s/cm^3/Hz] + **kwargs : optional arguments, i.e., ``n_e`` and ``gamma`` Returns ------- power : float, or 1D `~numpy.ndarray` - The calculated synchrotron power w.r.t. each input emissivity. + The calculated synchrotron power at each input frequency. Unit: [W/Hz] """ - return helper.calc_power(emissivity, volume=self.volume) + frequencies = np.asarray(frequencies) + if emissivity is None: + emissivity = self.calc_emissivity(frequencies=frequencies, + **kwargs) + else: + emissivity = np.asarray(emissivity) + if emissivity.shape != frequencies.shape: + raise ValueError("input 'frequencies' and 'emissivity' " + "do not match") + power = emissivity * (self.volume * AUC.kpc2cm**3) # [erg/s/Hz] + power *= 1e-7 # [erg/s/Hz] -> [W/Hz] + return power - def calc_flux(self, emissivity): + def calc_flux(self, frequencies, **kwargs): """ Calculate the synchrotron flux density (i.e., power *observed* - per unit frequency) from emissivity. + per unit frequency) of the halo, with k-correction considered. + + NOTE + ---- + The *k-correction* must be applied to the flux density (Sν) or + specific luminosity (Lν) because the redshifted object is emitting + flux in a different band than that in which you are observing. + And the k-correction depends on the spectrum of the object in + question. For any other spectrum (i.e., vLv != const.), the flux + density Sv is related to the specific luminosity Lv by: + Sv = (1+z) L_v(1+z) / (4π DL^2), + where + * L_v(1+z) is the specific luminosity emitting at frequency v(1+z), + * DL is the luminosity distance to the object at redshift z. + + Reference: Ref.[hogg1999],Eq.(22) Parameters ---------- - emissivity : float, or 1D `~numpy.ndarray` - The synchrotron emissivity at multiple frequencies. - Unit: [erg/s/cm^3/Hz] + frequencies : float, or 1D `~numpy.ndarray` + The frequencies where to calculate the flux density. + Unit: [MHz] + **kwargs : optional arguments, i.e., ``n_e`` and ``gamma`` Returns ------- flux : float, or 1D `~numpy.ndarray` - The calculated synchrotron flux w.r.t. each input emissivity. + The calculated flux density w.r.t. each input frequency. Unit: [Jy] = 1e-23 [erg/s/cm^2/Hz] = 1e-26 [W/m^2/Hz] """ - power = self.calc_power(emissivity) # [W/Hz] - return helper.calc_flux(power, z=self.z_obs) + z = self.z_obs + freqz = np.asarray(frequencies) * (1+z) + power = self.calc_power(freqz, **kwargs) # [W/Hz] + DL = COSMO.DL(self.z_obs) * AUC.Mpc2m # [m] + flux = 1e26 * (1+z) * power / (4*np.pi * DL*DL) # [Jy] + return flux - def calc_brightness_mean(self, emissivity, frequency, pixelsize=None): + def calc_brightness_mean(self, frequencies, flux=None, pixelsize=None, + **kwargs): """ Calculate the mean surface brightness (power observed per unit frequency and per unit solid angle) expressed in *brightness - temperature* at the specified frequencies from emissivity. + temperature* at the specified frequencies. + + NOTE + ---- + If the solid angle that the object extends is smaller than the + specified pixel area, then is is assumed to have size of 1 pixel. Parameters ---------- - emissivity : float, or 1D `~numpy.ndarray` - The synchrotron emissivity at multiple frequencies. - Unit: [erg/s/cm^3/Hz] - frequency : float, or 1D `~numpy.ndarray` - The frequencies where the synchrotron emissivity is calculated. + frequencies : float, or 1D `~numpy.ndarray` + The frequencies where to calculate the mean brightness temperature Unit: [MHz] + flux : float, or 1D `~numpy.ndarray`, optional + The flux density w.r.t. each input frequency. + Unit: [Jy] pixelsize : float, optional The pixel size of the output simulated sky image. + If not provided, then invoke above ``calc_flux()`` method to + calculate them. Unit: [arcsec] + **kwargs : optional arguments, i.e., ``n_e`` and ``gamma`` Returns ------- Tb : float, or 1D `~numpy.ndarray` - The mean surface brightness at each frequency. + The mean brightness temperature at each frequency. Unit: [K] <-> [Jy/pixel] """ + frequencies = np.asarray(frequencies) + if flux is None: + flux = self.calc_flux(frequencies=frequencies, **kwargs) # [Jy] + else: + flux = np.asarray(flux) + if flux.shape != frequencies.shape: + raise ValueError("input 'frequencies' and 'flux' do not match") + omega = np.pi * self.angular_radius**2 # [arcsec^2] - flux = self.calc_flux(emissivity) - return helper.calc_brightness_mean(flux, frequency=frequency, - omega=omega, pixelsize=pixelsize) + if pixelsize and (omega < pixelsize**2): + omega = pixelsize ** 2 # [arcsec^2] + logger.warning("Object size < 1 pixel; force to be 1 pixel!") + + Tb = Fnu_to_Tb(flux, omega, frequencies) # [K] + return Tb def fp_injection(self, gamma, t=None): """ diff --git a/fg21sim/extragalactic/clusters/helper.py b/fg21sim/extragalactic/clusters/helper.py index e553da9..6b835a0 100644 --- a/fg21sim/extragalactic/clusters/helper.py +++ b/fg21sim/extragalactic/clusters/helper.py @@ -40,7 +40,6 @@ from ...share import CONFIGS, COSMO from ...utils.units import (Units as AU, Constants as AC, UnitConversions as AUC) -from ...utils.convert import Fnu_to_Tb from ...utils.draw import circle from ...utils.transform import circle2ellipse @@ -304,107 +303,6 @@ def magnetic_field(mass): return B -def calc_power(emissivity, volume): - """ - Calculate the synchrotron power (i.e., power *emitted* per unit - frequency) from emissivity, which assumed to be uniform within - the volume. - - NOTE - ---- - The calculated power (a.k.a. spectral luminosity) is in units of - [W/Hz] which is common in radio astronomy, instead of [erg/s/Hz]. - 1 [W] = 1e7 [erg/s] - - Parameters - ---------- - emissivity : float, or 1D `~numpy.ndarray` - The synchrotron emissivity at multiple frequencies. - Unit: [erg/s/cm^3/Hz] - volume : float - The volume of the radio halo - Unit: [kpc^3] - - Returns - ------- - power : float, or 1D `~numpy.ndarray` - The calculated synchrotron power w.r.t. each input emissivity. - Unit: [W/Hz] - """ - emissivity = np.asarray(emissivity) - power = emissivity * (volume * AUC.kpc2cm**3) # [erg/s/Hz] - power *= 1e-7 # [erg/s/Hz] -> [W/Hz] - return power - - -def calc_flux(power, z): - """ - Calculate the synchrotron flux density (i.e., power *observed* - per unit frequency) from radio power at a certain redshift (i.e., - distance). - - Parameters - ---------- - power : float, or 1D `~numpy.ndarray` - The synchrotron power at multiple frequencies. - Unit: [W/Hz] - - Returns - ------- - flux : float, or 1D `~numpy.ndarray` - The calculated synchrotron flux w.r.t. each input power. - Unit: [Jy] = 1e-23 [erg/s/cm^2/Hz] = 1e-26 [W/m^2/Hz] - """ - DL = COSMO.DL(z) * AUC.Mpc2m # [m] - flux = 1e26 * power / (4*np.pi * DL*DL) # [Jy] - return flux - - -def calc_brightness_mean(flux, frequency, omega, pixelsize=None): - """ - Calculate the mean surface brightness (power observed per unit - frequency and per unit solid angle) expressed in *brightness - temperature* at the specified frequencies from flux. - - NOTE - ---- - If the solid angle that the object extends is smaller than the - specified pixel area, then is is assumed to have size of 1 pixel. - - Parameters - ---------- - flux : float, or 1D `~numpy.ndarray` - The synchrotron flux densities at multiple frequencies. - Unit: [Jy] - frequency : float, or 1D `~numpy.ndarray` - The frequencies where the above flux calculated. - Unit: [MHz] - omega : float - The sky coverage (angular size) of the object. - Unit: [arcsec^2] - pixelsize : float, optional - The pixel size of the output simulated sky image. - Unit: [arcsec] - - Returns - ------- - Tb : float, or 1D `~numpy.ndarray` - The mean surface brightness at each frequency. - Unit: [K] <-> [Jy/pixel] - """ - if pixelsize and (omega < pixelsize**2): - omega = pixelsize ** 2 # [arcsec^2] - logger.warning("Object sky coverage < 1 pixel; force to be 1 pixel") - - Tb = [Fnu_to_Tb(Fnu, omega, freq) - for Fnu, freq in zip(np.array(flux, ndmin=1), - np.array(frequency, ndmin=1))] - if len(Tb) == 1: - return Tb[0] - else: - return np.array(Tb) - - def halo_rprofile(re, num_re=5, I0=1.0): """ Generate the radial profile of a halo. diff --git a/fg21sim/extragalactic/clusters/main.py b/fg21sim/extragalactic/clusters/main.py index 70d5e6b..e6c6154 100644 --- a/fg21sim/extragalactic/clusters/main.py +++ b/fg21sim/extragalactic/clusters/main.py @@ -303,9 +303,10 @@ class GalaxyClusters: halo.set_electron_spectrum(hdict["n_e"]) emissivity = halo.calc_emissivity(frequencies=self.frequencies) - power = halo.calc_power(emissivity) - flux = halo.calc_flux(emissivity) - Tb_mean = halo.calc_brightness_mean(emissivity, self.frequencies, + power = halo.calc_power(self.frequencies, emissivity=emissivity) + # k-correction considered + flux = halo.calc_flux(self.frequencies) + Tb_mean = halo.calc_brightness_mean(self.frequencies, flux=flux, pixelsize=self.sky.pixelsize) # Update or add new items hdict["frequency"] = self.frequencies # [MHz] |