aboutsummaryrefslogtreecommitdiffstats
path: root/fg21sim/extragalactic
diff options
context:
space:
mode:
Diffstat (limited to 'fg21sim/extragalactic')
-rw-r--r--fg21sim/extragalactic/clusters/halo.py119
-rw-r--r--fg21sim/extragalactic/clusters/helper.py102
-rw-r--r--fg21sim/extragalactic/clusters/main.py7
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]