From 13229ee2a69c37cc3c5a1f002efce1ff0d66de98 Mon Sep 17 00:00:00 2001 From: Aaron LI Date: Wed, 20 Sep 2017 22:13:25 +0800 Subject: clusters/emission: Rewrite synchrotron kernel function Add asymptotic functions to calculate the values beyond the interpolation bounds (e.g., <1e-3 and >10), otherwise, the calculated synchrotron emissivity is overestimated at the higher frequencies. By rewrite this synchrotron kernel function, the calculated results is consistent with the theoretical/analytical results, e.g., the synchrotron radiation of a population of electrons of power-law index n is also a power-law with index (n-1)/2. --- fg21sim/extragalactic/clusters/emission.py | 47 ++++++++++++++++++++++++------ 1 file changed, 38 insertions(+), 9 deletions(-) diff --git a/fg21sim/extragalactic/clusters/emission.py b/fg21sim/extragalactic/clusters/emission.py index baab1b9..d620f80 100644 --- a/fg21sim/extragalactic/clusters/emission.py +++ b/fg21sim/extragalactic/clusters/emission.py @@ -17,6 +17,11 @@ References Essential Radio Astronomy https://science.nrao.edu/opportunities/courses/era/ Chapter.5 + +.. [you1998] + You 1998 + The Radiation Mechanisms in Astrophysics, 2nd Edition, Beijing + Sec.4.2.3, p.187 """ import logging @@ -32,20 +37,27 @@ from ...utils.units import (Units as AU, Constants as AC) logger = logging.getLogger(__name__) -def _interp_sync_kernel(xmin=1e-5, xmax=20.0, xsample=128): +def _interp_sync_kernel(xmin=1e-3, xmax=10.0, xsample=256): """ Sample the synchrotron kernel function at the specified X positions and make an interpolation, to optimize the speed when invoked to calculate the synchrotron emissivity. + WARNING + ------- + Do NOT simply bound the synchrotron kernel within the specified + [xmin, xmax] range, since it decreases as a power law of index + 1/3 at the left end, and decreases exponentially at the right end. + Bounding it with interpolation will cause the synchrotron emissivity + been *overestimated* on the higher frequencies. + Parameters ---------- xmin, xmax : float, optional The lower and upper cuts for the kernel function. - Default: [1e-5, 20.0] + Default: [1e-3, 10.0] xsample : int, optional Number of samples within [xmin, xmax] used to do interpolation. - NOTE: The kernel function is quiet smooth and slow-varying. Returns ------- @@ -56,9 +68,8 @@ def _interp_sync_kernel(xmin=1e-5, xmax=20.0, xsample=128): Fxx = [xp * integrate.quad(lambda t: scipy.special.kv(5/3, t), a=xp, b=np.inf)[0] for xp in xx] - F_interp = interpolate.interp1d( - xx, Fxx, kind="quadratic", bounds_error=False, - fill_value=(Fxx[0], Fxx[-1]), assume_sorted=True) + F_interp = interpolate.interp1d(xx, Fxx, kind="quadratic", + bounds_error=True, assume_sorted=True) return F_interp @@ -78,8 +89,13 @@ class SynchrotronEmission: The assumed uniform magnetic field within the cluster ICM. Unit: [uG] """ - # The interpolated synchrotron kernel function ``F(x)``. - F_interp = _interp_sync_kernel() + # The interpolated synchrotron kernel function ``F(x)`` within + # the specified range. + # NOTE: See the *WARNING* above. + F_xmin = 1e-3 + F_xmax = 10.0 + F_xsample = 256 + F_interp = _interp_sync_kernel(F_xmin, F_xmax, F_xsample) def __init__(self, gamma, n_e, B): self.gamma = np.asarray(gamma) @@ -152,8 +168,21 @@ class SynchrotronEmission: ------- y : `~numpy.ndarray` Calculated kernel function values. + + References: Ref.[you1998] """ - return cls.F_interp(x) + x = np.array(x, ndmin=1) + y = np.zeros(x.shape) + idx = (x >= cls.F_xmin) & (x <= cls.F_xmax) + y[idx] = cls.F_interp(x[idx]) + # Left end: power law of index 1/3 + idx = (x < cls.F_xmin) + A = cls.F_interp(cls.F_xmin) + y[idx] = A * (x[idx] / cls.F_xmin)**(1/3) + # Right end: exponentially decrease + idx = (x > cls.F_xmax) + y[idx] = (0.5*np.pi * x[idx])**0.5 * np.exp(-x[idx]) + return y def emissivity(self, frequencies): """ -- cgit v1.2.2