From 93dc96a0a6a336419b2e80b221257a8a7532019a Mon Sep 17 00:00:00 2001
From: Aaron LI <aaronly.me@outlook.com>
Date: Mon, 16 Jan 2017 21:10:30 +0800
Subject: Add clusters/emission.py: Calculate synchrotron emissivity
Calculate the synchrotron emissivity from the given electron spectrum at
specified frequency.
Add the electron charge to units.Constants
---
fg21sim/extragalactic/clusters/emission.py | 117 +++++++++++++++++++++++++++++
fg21sim/extragalactic/clusters/units.py | 6 ++
2 files changed, 123 insertions(+)
create mode 100644 fg21sim/extragalactic/clusters/emission.py
(limited to 'fg21sim/extragalactic/clusters')
diff --git a/fg21sim/extragalactic/clusters/emission.py b/fg21sim/extragalactic/clusters/emission.py
new file mode 100644
index 0000000..e71632b
--- /dev/null
+++ b/fg21sim/extragalactic/clusters/emission.py
@@ -0,0 +1,117 @@
+# Copyright (c) 2017 Weitian LI <liweitianux@live.com>
+# MIT license
+
+"""
+Calculate the synchrotron emission for simulated radio halos.
+"""
+
+import logging
+
+import numpy as np
+import scipy.integrate
+import scipy.special
+
+from .units import (Units as AU, Constants as AC)
+
+
+logger = logging.getLogger(__name__)
+
+
+class SynchrotronEmission:
+ """
+ Calculate the synchrotron emission from a given population
+ of electrons.
+
+ References
+ ----------
+ [1] Cassano & Brunetti 2005, MNRAS, 357, 1313
+ http://adsabs.harvard.edu/abs/2005MNRAS.357.1313C
+ Appendix.C
+ """
+ def __init__(self, B):
+ # Uniform magnetic field strength
+ self.B = B # [uG]
+
+ @property
+ def frequency_larmor(self):
+ """
+ Electron Larmor frequency:
+ ν_L = e * B / (2*π * m0 * c) = e * B / (2*π * mec)
+
+ Unit: MHz
+ """
+ coef = AC.e / (2*np.pi * AU.mec) # [Hz/G]
+ coef *= 1e-12 # [MHz/uG]
+ nu = coef * self.B # [MHz]
+ return nu
+
+ def frequency_crit(self, p, theta=np.pi/2):
+ """
+ Synchrotron critical frequency.
+
+ Critical frequency:
+ ν_c = (3/2) * γ^2 * sin(θ) * ν_L
+
+ Parameters
+ ----------
+ p : float
+ Electron momentum (unit: mec), i.e., Lorentz factor γ
+ theta : float, optional
+ The angle between the electron velocity and the magnetic field.
+ (unit: radian)
+
+ Returns
+ -------
+ nu : float
+ Critical frequency, unit: MHz
+ """
+ nu_L = self.frequency_larmor
+ nu = (3/2) * p**2 * np.sin(theta) * nu_L
+ return nu
+
+ def F(self, x):
+ """
+ Synchrotron kernel function.
+ """
+ val = x * scipy.integrate.quad(lambda t: scipy.special.kv(5/3, t),
+ a=x, b=np.inf)[0]
+ return val
+
+ def emissivity(self, nu, p, n_e):
+ """
+ Calculate the synchrotron emissivity (power emitted per volume
+ and per frequency) at the specified frequency from the given
+ electron number & energy distribution
+
+ Parameters
+ ----------
+ nu : float
+ Frequency (unit: MHz) where to calculate the emissivity.
+ p : `~numpy.ndarray`
+ The momentum grid adopted when solving the Fokker-Planck equation.
+ Unit: [mec]
+ n_e : `~numpy.ndarray`
+ Electron spectrum by solving the Fokker-Planck equation.
+ Unit: [cm^-3 mec^-1]
+
+ Returns
+ -------
+ j_nu : float
+ Synchrotron emissivity at frequency ``nu``.
+ Unit: [erg/s/cm^3/Hz]
+ """
+ def func(theta, _p, _n_e):
+ nu_c = self.frequency_crit(_p, theta)
+ x = nu / nu_c
+ return (np.sin(theta)**2 * _n_e * self.F(x))
+
+ coef = np.sqrt(3) * AC.e**3 * self.B / AC.c # multiplied a [mec]
+ func_p = np.zeros(p.shape)
+ for i in range(len(p)):
+ # Integrate over ``theta``
+ func_p[i] = scipy.integrate.quad(
+ lambda t: func(t, p[i], n_e[i]),
+ a=0, b=np.pi/2)[0]
+ # Integrate over ``p``
+ j_nu = coef * scipy.integrate.trapz(func_p, p)
+ return j_nu
diff --git a/fg21sim/extragalactic/clusters/units.py b/fg21sim/extragalactic/clusters/units.py
index 01891fa..a7e4ef5 100644
--- a/fg21sim/extragalactic/clusters/units.py
+++ b/fg21sim/extragalactic/clusters/units.py
@@ -30,10 +30,13 @@ class UnitConversions:
Hold the conversion relations directly to avoid repeated/unnecessary
calculations.
"""
+ # Mass
Msun2g = au.solMass.to(au.g)
g2Msun = au.g.to(au.solMass)
+ # Time
Gyr2s = au.Gyr.to(au.s)
s2Gyr = au.s.to(au.Gyr)
+ # Length
kpc2cm = au.kpc.to(au.cm)
cm2kpc = au.cm.to(au.kpc)
Mpc2cm = au.Mpc.to(au.cm)
@@ -43,6 +46,7 @@ class UnitConversions:
kpc2km = au.kpc.to(au.km)
km2kpc = au.km.to(au.kpc)
km2cm = au.km.to(au.cm)
+ # Energy
keV2erg = au.keV.to(au.erg)
@@ -60,6 +64,8 @@ class Constants:
u = ac.u.cgs.value # [g]
# Gravitational constant
G = ac.G.cgs.value # [cm^3/g/s^2]
+ # Electron charge
+ e = ac.e.gauss.value # [Fr] = [esu]
# Mean molecular weight
# Ref.: Ettori et al, 2013, Space Science Review, 177, 119-154, Eq.(6)
--
cgit v1.2.2