From 39a2a068fd74bfde3ad758a01f44dd3c6b72807e Mon Sep 17 00:00:00 2001
From: Aaron LI <aly@aaronly.me>
Date: Wed, 3 Jan 2018 22:16:55 +0800
Subject: clusters/halo: split emissivity/power/flux calculations into
HaloEmission
Add emission.HaloEmission class to calculate halo emissivity, power, flux,
brightness etc.
---
fg21sim/extragalactic/clusters/emission.py | 201 ++++++++++++++++++++++++++++-
fg21sim/extragalactic/clusters/halo.py | 174 -------------------------
2 files changed, 200 insertions(+), 175 deletions(-)
diff --git a/fg21sim/extragalactic/clusters/emission.py b/fg21sim/extragalactic/clusters/emission.py
index ddecab3..a798470 100644
--- a/fg21sim/extragalactic/clusters/emission.py
+++ b/fg21sim/extragalactic/clusters/emission.py
@@ -31,7 +31,11 @@ import numpy as np
import scipy.special
from scipy import integrate, interpolate
-from ...utils.units import (Units as AU, Constants as AC)
+from ...share import COSMO
+from ...utils.convert import Fnu_to_Tb
+from ...utils.units import (Units as AU,
+ UnitConversions as AUC,
+ Constants as AC)
logger = logging.getLogger(__name__)
@@ -223,3 +227,198 @@ class SynchrotronEmission:
return syncem[0]
else:
return syncem
+
+
+class HaloEmission:
+ """
+ Calculate the synchrotron emission of a (giant) radio halo.
+
+ Parameters
+ ----------
+ gamma : 1D `~numpy.ndarray`
+ The Lorentz factors γ of the electron spectrum.
+ n_e : 1D `~numpy.ndarray`
+ The electron spectrum (w.r.t. Lorentz factors γ).
+ Unit: [cm^-3]
+ B : float
+ The magnetic field strength.
+ Unit: [uG]
+ radius : float, optional
+ The radio halo radius.
+ Required to calculate the power.
+ Unit: [kpc]
+ redshift : float, optional
+ The redshift to the radio halo.
+ Required to calculate the flux, which also requires ``radius``.
+ """
+ def __init__(self, gamma, n_e, B, radius=None, redshift=None):
+ self.gamma = np.asarray(gamma)
+ self.n_e = np.asarray(n_e)
+ self.B = B
+ self.radius = radius
+ self.redshift = redshift
+
+ @property
+ def angular_radius(self):
+ """
+ The angular radius of the radio halo.
+ Unit: [arcsec]
+ """
+ if self.redshift is None:
+ raise RuntimeError("parameter 'redshift' is required")
+ if self.radius is None:
+ raise RuntimeError("parameter 'radius' is required")
+
+ DA = COSMO.DA(self.redshift) * 1e3 # [Mpc] -> [kpc]
+ theta = self.radius / DA # [rad]
+ return theta * AUC.rad2arcsec
+
+ @property
+ def volume(self):
+ """
+ The halo volume.
+ Unit: [kpc^3]
+ """
+ if self.radius is None:
+ raise RuntimeError("parameter 'radius' is required")
+
+ return (4*np.pi/3) * self.radius**3
+
+ def calc_emissivity(self, frequencies):
+ """
+ Calculate the synchrotron emissivity for the derived electron
+ spectrum.
+
+ Parameters
+ ----------
+ frequencies : float, or 1D `~numpy.ndarray`
+ The frequencies where to calculate the synchrotron emissivity.
+ Unit: [MHz]
+
+ Returns
+ -------
+ emissivity : float, or 1D `~numpy.ndarray`
+ The calculated synchrotron emissivity at each specified
+ frequency.
+ Unit: [erg/s/cm^3/Hz]
+ """
+ syncem = SynchrotronEmission(gamma=self.gamma, n_e=self.n_e, B=self.B)
+ emissivity = syncem.emissivity(frequencies)
+ return emissivity
+
+ def calc_power(self, frequencies, emissivity=None):
+ """
+ Calculate the halo synchrotron power (i.e., power *emitted* per
+ 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
+ ----------
+ 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]
+
+ Returns
+ -------
+ power : float, or 1D `~numpy.ndarray`
+ The calculated synchrotron power at each input frequency.
+ Unit: [W/Hz]
+ """
+ frequencies = np.asarray(frequencies)
+ if emissivity is None:
+ emissivity = self.calc_emissivity(frequencies=frequencies)
+ else:
+ emissivity = np.asarray(emissivity)
+ power = emissivity * (self.volume * AUC.kpc2cm**3) # [erg/s/Hz]
+ power *= 1e-7 # [erg/s/Hz] -> [W/Hz]
+ return power
+
+ def calc_flux(self, frequencies):
+ """
+ Calculate the synchrotron flux density (i.e., power *observed*
+ 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)
+
+ Returns
+ -------
+ flux : float, or 1D `~numpy.ndarray`
+ 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]
+ """
+ if self.redshift is None:
+ raise RuntimeError("parameter 'redshift' is required")
+
+ freqz = np.asarray(frequencies) * (1+self.redshift)
+ power = self.calc_power(freqz) # [W/Hz]
+ DL = COSMO.DL(self.redshift) * AUC.Mpc2m # [m]
+ flux = 1e26 * (1+self.redshift) * power / (4*np.pi * DL*DL) # [Jy]
+ return flux
+
+ def calc_brightness_mean(self, frequencies, flux=None, 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.
+
+ 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
+ ----------
+ 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]
+
+ Returns
+ -------
+ Tb : float, or 1D `~numpy.ndarray`
+ 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) # [Jy]
+ else:
+ flux = np.asarray(flux)
+ omega = np.pi * self.angular_radius**2 # [arcsec^2]
+ if pixelsize and (omega < pixelsize**2):
+ omega = pixelsize ** 2 # [arcsec^2]
+ logger.warning("Halo size < 1 pixel; force to be 1 pixel!")
+
+ Tb = Fnu_to_Tb(flux, omega, frequencies) # [K]
+ return Tb
diff --git a/fg21sim/extragalactic/clusters/halo.py b/fg21sim/extragalactic/clusters/halo.py
index 4e5f30f..29262b8 100644
--- a/fg21sim/extragalactic/clusters/halo.py
+++ b/fg21sim/extragalactic/clusters/halo.py
@@ -62,12 +62,10 @@ import numpy as np
from . import helper
from .solver import FokkerPlanckSolver
-from .emission import SynchrotronEmission
from ...share import CONFIGS, COSMO
from ...utils.units import (Units as AU,
UnitConversions as AUC,
Constants as AC)
-from ...utils.convert import Fnu_to_Tb
logger = logging.getLogger(__name__)
@@ -494,178 +492,6 @@ class RadioHalo:
else:
raise ValueError("given electron spectrum has wrong shape!")
- def calc_emissivity(self, frequencies, n_e=None, gamma=None, B=None):
- """
- Calculate the synchrotron emissivity for the derived electron
- spectrum.
-
- Parameters
- ----------
- frequencies : float, or 1D `~numpy.ndarray`
- The frequencies where to calculate the synchrotron emissivity.
- Unit: [MHz]
- n_e : 1D `~numpy.ndarray`, optional
- The electron spectrum (w.r.t. Lorentz factors γ).
- If not provided, then use the cached ``self.electron_spec``
- that was solved at above.
- Unit: [cm^-3]
- gamma : 1D `~numpy.ndarray`, optional
- The Lorentz factors γ of the electron spectrum.
- If not provided, then use ``self.gamma``.
- B : float, optional
- The magnetic field strength.
- If not provided, then use ``self.B_obs``.
- Unit: [uG]
-
- Returns
- -------
- emissivity : float, or 1D `~numpy.ndarray`
- The calculated synchrotron emissivity at each specified
- frequency.
- Unit: [erg/s/cm^3/Hz]
- """
- if n_e is None:
- n_e = self.electron_spec
- if gamma is None:
- gamma = self.gamma
- if B is None:
- B = self.B_obs
- syncem = SynchrotronEmission(gamma=gamma, n_e=n_e, B=B)
- emissivity = syncem.emissivity(frequencies)
- return emissivity
-
- def calc_power(self, frequencies, emissivity=None, **kwargs):
- """
- Calculate the halo synchrotron power (i.e., power *emitted* per
- 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
- ----------
- 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``, ``gamma``, and ``B``.
-
- Returns
- -------
- power : float, or 1D `~numpy.ndarray`
- The calculated synchrotron power at each input frequency.
- Unit: [W/Hz]
- """
- 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, frequencies, **kwargs):
- """
- Calculate the synchrotron flux density (i.e., power *observed*
- 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
- ----------
- frequencies : float, or 1D `~numpy.ndarray`
- The frequencies where to calculate the flux density.
- Unit: [MHz]
- **kwargs : optional arguments, i.e., ``n_e``, ``gamma``, and ``B``.
-
- Returns
- -------
- flux : float, or 1D `~numpy.ndarray`
- 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]
- """
- 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, 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.
-
- 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
- ----------
- 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``, ``gamma``, and ``B``.
-
- Returns
- -------
- Tb : float, or 1D `~numpy.ndarray`
- 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]
- if pixelsize and (omega < pixelsize**2):
- omega = pixelsize ** 2 # [arcsec^2]
- logger.warning("Halo 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):
"""
Electron injection (rate) term for the Fokker-Planck equation.
--
cgit v1.2.2