Created
September 27, 2014 20:48
-
-
Save cdeil/1b9e58f2e9e5c4096f03 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
"""Simple analytical PWN evolution and spectral models. | |
- Plot total energy released by a pulsar via spindown | |
as a function of age. | |
See Equation (9) and Figure (2) from Mattana et al. | |
http://adsabs.harvard.edu/abs/2009ApJ...694...12M | |
- Compute exponential cutoff broken power-law electron | |
spectrum resulting from a continuous injection for time T | |
with luminosity L and a power law with spectral index gamma | |
and cutoff E_cut into magnetic field B and photon field U_ph. | |
- Compute resulting synchrotron and inverse Compton spectra. | |
References: | |
Hinton: http://arjournals.annualreviews.org/doi/pdf/10.1146/annurev-astro-082708-101816 | |
Reynolds: http://adsabs.harvard.edu/abs/2009ApJ...703..662R | |
Beckmann: http://asd.gsfc.nasa.gov/Volker.Beckmann/school/download/Longair_Radiation2.pdf | |
Hoppe: Stefan Hoppe's thesis (@todo: remove this, only use formulae from refereed papers!) | |
Mattana: http://adsabs.harvard.edu/abs/2009ApJ...694...12M | |
""" | |
from numpy import exp, where | |
from scipy.integrate import quad | |
from mpmath import gammainc | |
NO_EMIN = 1e-10 | |
NO_EMAX = 1e20 | |
electron_rest_energy = 510.999e3 # eV | |
year_to_sec = 365. * 24 * 3600 | |
sec_to_year = 1 / year_to_sec | |
eV_to_erg = 1.6022e-12 | |
erg_to_eV = 1 / eV_to_erg | |
def L_to_F(L, d, eV_to_erg_power=1): | |
""" | |
@param L: Luminosity | |
@param d: distance in kpc | |
@return: flux at distance d in units of L per cm^2""" | |
d_cm = 3.08568025e21 * d | |
return eV_to_erg ** eV_to_erg_power * L / (4 * pi * d_cm ** 2) | |
def B_energy_density(B): | |
"""@param B: magnetic field strength in uG | |
@return: energy density in eV cm^-3 | |
@note The basic formula B ** 2 / (8 * pi) needs Gauss | |
as input unit and has erg m^3 as output""" | |
return erg_to_eV * (1e-6 * B) ** 2 / (8 * pi) | |
''' | |
Doesn't give correct results and CombinedSyIC().cooling_time() | |
is simpler. | |
@todo Remove at some point. | |
class Cooling: | |
def __init__(self, B, U_ph=0.26): | |
"""Combined cooling energy from synchtrotron and | |
IC losses. Mattana Eq. (12). | |
@param B: magnetic field strength (uG) | |
@param U_ph: photon field energy density (eV cm^-3)""" | |
self.B = float(B) | |
self.U_B = B_energy_density(self.B) | |
self.U_ph = U_ph | |
self.U = self.U_B + self.U_ph | |
self.zeta = self.U_ph / self.U_B | |
self.factor = (1e7 * electron_rest_energy * | |
24.5 * 1e3 / (1 + self.zeta) / | |
(self.B / 10) ** 2) | |
def energy(self, T): | |
"""@param T: cooling time (yr) | |
@return: energy with cooling time T (eV)""" | |
return self.factor / T | |
def time(self, E): | |
"""@param E: energy (eV) | |
@return cooling time (yr)""" | |
return self.factor * E | |
''' | |
class ECPL: | |
"""Representation of exponential cutoff power law | |
f(E) = norm * E ** -gamma * exp(-E / Ecut) | |
Used to compute the integral""" | |
def __init__(self, norm=1, gamma=2, Ecut=1, Eref=1): | |
self.norm = float(norm) | |
self.gamma = float(gamma) | |
self.Ecut = float(Ecut) | |
self.Eref = float(Eref) | |
def __call__(self, E, power=0): | |
return(E ** power * self.norm * (E / self.Eref) ** -self.gamma * | |
exp(-E / self.Ecut)) | |
def indef_integral(self, E, power=0): | |
"""Computed with www.wolframalpha.com | |
Input (needs simple variable names): | |
int t^p n (t/r)^-g * e^(-t/a) dt | |
Output: | |
- a n t^p (t/r)^-g (t/a)^(g-p) gamma(p-g+1, t/a) | |
""" | |
term1 = -self.Ecut * self.norm * E ** power | |
term2 = (E / self.Eref) ** -self.gamma | |
term3 = (E / self.Ecut) ** (self.gamma - power) | |
arg1 = power - self.gamma + 1 | |
arg2 = E / self.Ecut | |
#print('E = %g, Ecut = %g' % (E, self.Ecut)) | |
#print('arg1 = %g, arg2 = %g' % (arg1, arg2)) | |
# @note scipy.special.gammainc is the regularized version | |
# and will not work for negative first arguments. | |
# We need mpmath.gammainc here! | |
term4 = gammainc(arg1, arg2) | |
#print('term1 = %g, term2 = %g, term3 = %g' % | |
# (term1, term2, term3)) | |
return term1 * term2 * term3 * term4 | |
def analytic_integral(self, Emin, Emax, power=0): | |
Imax = self.indef_integral(Emax, power) | |
Imin = self.indef_integral(Emin, power) | |
#print 'Imax = %g, Imin = %g' % (Imax, Imin) | |
return Imax - Imin | |
def numeric_integral(self, Emin, Emax, power): | |
"""This is numerically not very stable!""" | |
return quad(self.__call__, Emin, Emax, (power))[0] | |
def integral(self, Emin, Emax, power=0): | |
return self.analytic_integral(Emin, Emax, power) | |
def test_integral(self, Emin=0.1, Emax=10, power=0): | |
analytic = self.analytic_integral(Emin, Emax, power) | |
numeric = self.numeric_integral(Emin, Emax, power) | |
print('numeric = %g, analytic = %g, ratio = %g' % | |
(numeric, analytic, numeric / analytic)) | |
class Model: | |
"""Exponential cutoff broken power law | |
@param norm: differential flux at Ebreak | |
http://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/source_models.html#BPLExpCutoff""" | |
def __init__(self, norm=1, gamma=2, gamma2=2, | |
Ebreak=1e12, Ecut=1e15): | |
self.norm = float(norm) | |
self.gamma = float(gamma) | |
self.gamma2 = float(gamma2) | |
self.Ebreak = float(Ebreak) | |
self.Ecut = float(Ecut) | |
def __call__(self, E, power=0): | |
term1 = E ** power * self.norm * exp(-E / self.Ecut) | |
term2 = where(E < self.Ebreak, | |
(E / self.Ebreak) ** (-self.gamma), | |
(E / self.Ebreak) ** (-self.gamma2)) | |
return term1 * term2 | |
def numeric_integral(self, Emin=NO_EMIN, Emax=NO_EMAX, power=0): | |
"""This is numerically not very stable!""" | |
return quad(self.__call__, Emin, Emax, (power))[0] | |
def analytic_integral(self, Emin=NO_EMIN, Emax=NO_EMAX, power=0): | |
"""ECPL integral computed analytically""" | |
# Make sure the ranges make sense | |
Esplit = self.Ebreak if Emin <= self.Ebreak <= Emax else Emin | |
#Eref = self.Ebreak | |
ecpl = ECPL(self.norm, self.gamma, self.Ecut, self.Ebreak) | |
term1 = ecpl.integral(Emin, Esplit, power) | |
ecpl.gamma = self.gamma2 | |
term2 = ecpl.integral(Esplit, Emax, power) | |
#print('Emin = %g, Esplit = %g, Emax = %g' % | |
# (Emin, Esplit, Emax)) | |
#print('term1 = %g, term2 = %g' % (term1, term2)) | |
return term1 + term2 | |
def integral(self, Emin=NO_EMIN, Emax=NO_EMAX, power=0): | |
return self.analytic_integral(Emin, Emax, power) | |
def set_total_energy(self, total_energy, Emin=NO_EMIN, Emax=NO_EMAX): | |
"""Set norm so that the total energy in the spectrum matches the given value.""" | |
self.norm *= total_energy / self.total_energy(Emin, Emax) | |
def total_energy(self, Emin=NO_EMIN, Emax=NO_EMAX): | |
return self.integral(Emin, Emax, 1) | |
def __str__(self): | |
return ('norm = %g, gamma = %g, gamma2 = %g, Ebreak = %g, Ecut = %g' % | |
(self.norm, self.gamma, self.gamma2, self.Ebreak, self.Ecut)) | |
def test_integral(self, Emin=0.1, Emax=10, power=0): | |
analytic = self.analytic_integral(Emin, Emax, power) | |
numeric = self.numeric_integral(Emin, Emax, power) | |
print('Emin = %g, Emax = %g, power = %g, numeric = %g, analytic = %g, ratio = %g' % | |
(Emin, Emax, power, numeric, analytic, numeric / analytic)) | |
class ContinuousInjectionElectrons: | |
def __init__(self, L=1e37, T=1e4, U_ph=0.26, B=100, | |
gamma=2, Ecut=1e15, Emin=NO_EMIN, Emax=NO_EMAX): | |
""" | |
@param L: source power in erg s^-1 | |
@param T: age in years | |
@param U_ph: photon field energy density (eV cm^-3) | |
@param B: magnetic field strength (1e-6 Gauss) | |
""" | |
self.L = float(L) | |
self.T = float(T) | |
self.U_ph = float(U_ph) | |
self.B = float(B) | |
self.gamma = float(gamma) | |
self.Ecut = float(Ecut) | |
self.Emin, self.Emax = float(Emin), float(Emax) | |
# The spectral index steepens by 1 through cooling | |
self.gamma2 = gamma + 1 | |
self.U_B = B_energy_density(B) | |
#self.Ebreak = Cooling(self.B, self.U_ph).energy(self.T) | |
self.Ebreak = CombinedSyIC(self.B, self.U_ph).cooling_energy(self.T) | |
# Compute the total energy in the electrons | |
self.injected_energy = erg_to_eV * self.L * self.T * year_to_sec # in eV | |
self.set_norm() | |
def set_norm(self): | |
# Here is the procedure we use to compute the norm of the cooled spectrum. | |
# - an uncooled spectrum is created with arbitrary norm = 1. | |
# - then its norm is scaled such that its total energy matches the total injected energy | |
# - the same norm is taken for the cooled spectrum, since the norm is the flux | |
# at the break energy and there the cooled and uncooled spectra match | |
uncooled_spectrum = Model(1, self.gamma, self.gamma, self.Ebreak, self.Ecut) | |
uncooled_spectrum.set_total_energy(self.injected_energy, self.Emin, self.Emax) | |
norm = uncooled_spectrum.norm | |
uncooled_total_energy = uncooled_spectrum.integral(self.Emin, self.Emax, 1) | |
self.spectrum = Model(uncooled_spectrum.norm, self.gamma, self.gamma2, self.Ebreak, self.Ecut) | |
self.total_energy = self.spectrum.integral(self.Emin, self.Emax, 1) | |
self.remaining_fraction = self.total_energy / self.injected_energy | |
print('norm = %g, injected_energy = %g, uncooled_total_energy = %g, total_energy = %g, remaining_fraction = %g' % | |
(norm, self.injected_energy, uncooled_total_energy, self.total_energy, self.remaining_fraction)) | |
def __call__(self, E, power=0): | |
return self.spectrum(E, power) | |
Electrons = ContinuousInjectionElectrons | |
class RadiationMechanism: | |
"""Common stuff for Synchrotron and InverseCompton | |
C1 and C2 are defined by the following relations and are | |
computed in the Synchrotron and InverseCompton constructors. | |
The electron spectrum F_e(E_e) = dN_e / dE_e is related | |
to the photon spectrum F(E) = dN / (dE dt) as follows: | |
(1) E = C1 E_e ^ 2 (this is the correspondence of electron and photon energies) | |
(2) dE_e / dt = C2 E_e ^ 2 (this is the energy loss rate of one electron) | |
C1 is in eV^-1 | |
C2 is in eV^-1 s^-1 | |
""" | |
def compute_parameters(self, electron_spec): | |
gamma = self.index(electron_spec.gamma) | |
gamma2 = self.index(electron_spec.gamma2) | |
Ebreak = self.energy(electron_spec.Ebreak) | |
Ecut = self.energy(electron_spec.Ecut) | |
# The norm is computed by first setting to 1 and then | |
# scaling to the desired flux | |
# This is needed because of the cutoff, so norm != flux | |
E = electron_spec.Ebreak | |
if 0: | |
"""Correct method, currently being debugged""" | |
self.spectrum = Model(1, gamma, gamma2, Ebreak, Ecut) | |
expected_flux = self.flux(E, electron_spec(E)) | |
current_flux = self.spectrum(Ebreak) | |
self.spectrum.norm *= expected_flux / current_flux | |
else: | |
norm = self.flux(E, electron_spec(E)) | |
self.spectrum = Model(norm, gamma, gamma2, Ebreak, Ecut) | |
def energy(self, E_e): | |
"""Photon energy corresponding to given electron energy""" | |
return self.C1 * E_e ** 2 | |
def flux(self, E_e, F_e): | |
"""Photon spectrum flux at energy(E_e) for a given F_e(E_e). | |
Here is the derivation of this formula: | |
Take the derivative of (1): dE = 2 C1 E_e dE_e | |
Set the energy lost by electrons equal to the energy | |
gained by photons in a time bin dt: | |
dN E = dt dN_e (dE_e / dt) = dt dN_e (C2 E_e ^ 2) | |
Putting it all together we find: | |
F = dN / (dE dt) = (dt dN_e (C2 E_e ^ 2)) / (E dE dt) | |
= (dN_e C2 E_e ^ 2) / (C1 E_e ^ 2 2 C1 E_e dE_e) | |
= C2 / (2 C1 ^ 2 E_e) * (dN_e / dE_e) | |
= C2 / (2 C1 ^ 2 E_e) * F_e""" | |
return self.C2 / (2 * self.C1 ** 2 * E_e) * F_e | |
def cooling_time(self, E_e): | |
"""@param E_e: (eV) | |
@return: Cooling time (yr)""" | |
return sec_to_year / (E_e * self.C2) | |
def cooling_energy(self, t): | |
"""@param time: Cooling time (yr) | |
@return: Electron energy with that cooling time (eV)""" | |
return sec_to_year / (t * self.C2) | |
def index(self, gamma_e): | |
return (gamma_e + 1) / 2. | |
def __call__(self, E, power=0): | |
return self.spectrum(E, power) | |
def info(self): | |
print('U = %g, C1 = %g, C2 = %g' % | |
(self.U, self.C1, self.C2)) | |
class Synchrotron(RadiationMechanism): | |
def __init__(self, B, electron_spec=None): | |
""" | |
@param B: magnetic field strength (uG)""" | |
self.B = B # if B else electron_spec.B | |
self.U = B_energy_density(B) | |
# Constant relating photon and electron energy | |
self.C1 = 0.087 * 1e-24 * self.B # Hinton Eq. (11) | |
# Constant for electron energy loss rate | |
#self.C2 = year_to_sec * 1.3e7 * 1e-12 * self.B ** 2 # Hoppe Eq. (TODO) | |
#self.C2 = 1 * (2 / 3. * 2.37e-3) * self.B ** 2 # Reynolds Eq. (4) | |
# Note that the same formula for C2 applies for Synchrotron | |
# and InverseCompton, the only difference is that U is different! | |
self.C2 = 1e-12 * self.U / (year_to_sec * 3.1e5) # Hinton Eq. (5) | |
if electron_spec: self.compute_parameters(electron_spec) | |
class InverseCompton(RadiationMechanism): | |
def __init__(self, U=0.26, T=2.7, E_T=None, electron_spec=None): | |
""" | |
@param U: target photon field energy density (eV cm^-3) | |
@param T: temperature of blackbody target photons (K)""" | |
self.U = U | |
if E_T: | |
self.T = E_T / 8.6173e-5 | |
self.E_T = E_T | |
else: | |
self.T = T | |
# Target photon energy (eV) | |
"""There is some ambiguity in how to assign one photon energy | |
to a given temperature, since the maximum or median energy depends | |
if you look at the blackbody spectrum per wavelength or per energy | |
and if you do it on a linear or log scale. | |
The maximum per energy on a linear scale is given here | |
http://en.wikipedia.org/wiki/Black_body#Wien.27s_displacement_law | |
as f_max = (59 GHz K^-1) T, i.e. | |
E_max = (2.44e-4 eV K^-1) T""" | |
#self.E_T = 2.44e-4 * T # from Wien's law | |
self.E_T = 8.6173e-5 * T # using E_T = k * T | |
# Constant relating photon and electron energy | |
# @todo: The following relation is completely wrong in the KN regime. | |
# Maybe we can simply use (10) from Hinton & Hofmann, TeV Astronomy? | |
# This formula is from Hoppe's thesis (factor 1.4 higher than the formula by Hinton et al.) | |
#self.C1 = 5e9 * 1e4 * 1e-24 * self.E_T | |
self.C1 = 33 * 1e+12 * 1e-24 * self.E_T # Hinton Eq. (8) | |
# Constant for electron energy loss rate | |
self.C2 = 1e-12 * self.U / (year_to_sec * 3.1e5) # Hinton Eq. (5) | |
if electron_spec: self.compute_parameters(electron_spec) | |
class InverseComptonIC(InverseCompton): | |
def __init__(self, U=0.26, T=2.7, E_T=None, electron_spec=None): | |
super().__init__(U, T, E_T, electron_spec) | |
# @todo implement KN | |
def b(self, E_e): | |
return 4 * E_e * self.E_T / electron_rest_energy ** 2 # Hinton Section 2.2.1 | |
def f_KN(self, E_e): | |
return (1 + self.b(E_e)) ** -1.5 # Hinton Eq. (6) | |
def energy(self, E_e): | |
"""Photon energy corresponding to given electron energy""" | |
#return self.C1 * E_e ** 2 | |
raise NotImplementedError | |
def flux(self, E_e, F_e): | |
"""Photon spectrum flux at energy(E_e) for a given F_e(E_e). | |
Here is the derivation of this formula: | |
Take the derivative of (1): dE = 2 C1 E_e dE_e | |
Set the energy lost by electrons equal to the energy | |
gained by photons in a time bin dt: | |
dN E = dt dN_e (dE_e / dt) = dt dN_e (C2 E_e ^ 2) | |
Putting it all together we find: | |
F = dN / (dE dt) = (dt dN_e (C2 E_e ^ 2)) / (E dE dt) | |
= (dN_e C2 E_e ^ 2) / (C1 E_e ^ 2 2 C1 E_e dE_e) | |
= C2 / (2 C1 ^ 2 E_e) * (dN_e / dE_e) | |
= C2 / (2 C1 ^ 2 E_e) * F_e""" | |
raise NotImplemtedError | |
#return self.C2 / (2 * self.C1 ** 2 * E_e) * F_e | |
class CombinedSyIC(RadiationMechanism): | |
"""Useful for computing e.g. combined cooling time""" | |
def __init__(self, B, U=0.26, T=2.7, E_T=None, electron_spec=None): | |
"""@param B: magnetic field strength (uG) | |
@param U: target photon field energy density (eV cm^-3) | |
@param T: temperature of blackbody target photons (K)""" | |
sy = Synchrotron(B) | |
ic = InverseCompton(U, T, E_T) | |
self.C1 = None # Shall not be used | |
# Inverses of cooling times can be added: | |
# 1 / T = 1 / T_sy + 1 / T_ic | |
# and C ~ 1 / T | |
self.C2 = sy.C2 + ic.C2 | |
if electron_spec: self.compute_parameters(electron_spec) | |
class MattanaPWN: | |
def __init__(self, Edot0=1e39, tdec=1e2, B=10): | |
"""Reasonable values: | |
@param Edot0 1e38 -- 1e40 | |
@param tdec 1e2 -- 1e3""" | |
self.Edot0 = float(Edot0) | |
self.tdec = float(tdec) | |
self.B = float(B) | |
self.espec = CombinedSyIC(B) | |
def t_c(self, E): | |
"""Synchrotron plus IC on CMB cooling time""" | |
return self.espec.cooling_time(E) | |
def Edot(self, t): | |
"""Energy loss rate at time t""" | |
return self.Edot0 / (1 + t / self.tdec) ** 2 | |
def E(self, t2, t1=0): | |
"""Accumulated energy injected in time interval t1 to t2 | |
Computed with www.wolframalpha.com | |
Input (needs simple variable names): | |
int A / (1 + t/B)^2 dt from C to D | |
Output: | |
A B^2 (1 / (B+C) - 1 / (B + D)) | |
""" | |
#return self.Edot0 * self.tdec * (t / (t + self.tdec)) | |
term1 = self.Edot0 * self.tdec ** 2 | |
term2 = 1 / (self.tdec + t1) | |
term3 = 1 / (self.tdec + t2) | |
return term1 * (term2 - term3) | |
def n_u(self, t): | |
"""Injected energy since birth""" | |
return self.E(t) | |
def n_c_analytic(self, t, E): | |
t_c = self.t_c(E) | |
term1 = self.Edot0 * self.tdec ** 2 * t_c | |
term2 = t - t_c + self.tdec | |
term3 = t + self.tdec | |
return np.abs(term1 / term2 / term3) | |
def n_c_limit(self, t, E, match_at_tc=True, formula=0): | |
"""Limit for t >> max(tdec, t_c). | |
See Mattana after Eq. (13).""" | |
t_c = self.t_c(E) | |
if formula == 0: | |
value = self.Edot0 * self.tdec ** 2 * t_c * t ** -2 | |
else: | |
value = self.Edot(t) * t_c | |
if match_at_tc: | |
return value * self.n_u(t_c) / self.n_c_limit(t_c, E, False, formula) | |
else: | |
return value | |
def n_c(self, t, E): | |
"""Injected energy since one cooling time ago""" | |
t_c = self.t_c(E) | |
t_1 = np.where(t > t_c, t - t_c, 0) | |
#print 't = %g, t_c = %g, t_1 = %g' % (t, t_c, t_1) | |
return self.E(t, t_1) | |
def n(self, t, E): | |
return np.where(t < self.t_c(E), | |
self.n_u(t), | |
self.n_c_limit(t, E)) | |
def efficiency(self, t, E): | |
return self.n(t, E) / self.Edot(t) | |
def test(self): | |
"""Run a few consistency checks""" | |
t = 0 | |
dt = 1e-6 | |
print self.Edot(t) | |
print self.E(t+dt, t) / dt | |
print self.n_c(t + dt, 1e12) / dt | |
class SpectrumTester: | |
def __init__(self, electron_spec, sy_spec, ic_spec): | |
self.electron = electron_spec | |
self.sy = sy_spec | |
self.ic = ic_spec | |
def test_flux_ratio(self, E_e=None): | |
"""See Hoppe thesis (1.19)""" | |
E_e = E_e if E_e else self.electron.Ebreak | |
E_ic = self.ic.energy(E_e) | |
E_sy = self.sy.energy(E_e) | |
f_sy = E_sy ** 2 * self.sy(E_sy) | |
f_ic = E_ic ** 2 * self.ic(E_ic) | |
print ('E_e = %g, E_sy = %g, f_sy = %g, E_ic = %g, f_ic = %g' % | |
(E_e, E_sy, f_sy, E_ic, f_ic)) | |
B = self.sy.B | |
f_ratio = f_ic / f_sy | |
f_expected_ratio = 10 * B ** -2 | |
print('B = %g, f_ratio = %g, f_expected_ratio = %g, ratio = %g' % | |
(B, f_ratio, f_expected_ratio, f_ratio / f_expected_ratio)) | |
def print_basic_info(self): | |
self.sy.info() | |
print self.sy.spectrum | |
self.ic.info() | |
print self.ic.spectrum | |
class RelationsTests: | |
def __init__(self): | |
#self.test_ic_target_photon_energy() | |
#self.test_magnetic_field_energy_density() | |
#self.test_sy_energy() | |
#self.test_ic_energy() | |
self.test_sy_cooling_time() | |
self.test_ic_cooling_time() | |
self.test_combined_cooling_time() | |
#self.test_flux_ratio() | |
def test_ic_target_photon_energy(self): | |
"""Check against value given in caption of Figure 2 in Hinton et al.""" | |
print '*** test_CMB_photon_energy ***' | |
ic = InverseCompton() | |
E_T = ic.E_T | |
E_T_expected = 2.35e-4 | |
ratio = E_T / E_T_expected | |
print('E_T = %g, E_T_expected = %g, ratio = %g' % | |
(E_T, E_T_expected, ratio)) | |
def test_magnetic_field_energy_density(self, B=10): | |
"""Hoppe states right after Eq. (1.14) that | |
U = B ** 2 / (8 * pi) = 0.024 eV cm^-3 * (B / 1 uG)""" | |
print '*** test_magnetic_field_energy_density ***' | |
sy = Synchrotron(B=B) | |
U = sy.U | |
U_expected = 0.024 * B ** 2 | |
ratio = U / U_expected | |
print('B = %g, U = %g, U_expected = %g, ratio = %g' % | |
(sy.B, U, U_expected, ratio)) | |
def test_ic_energy(self): | |
"""Check equation (9) of | |
http://arjournals.annualreviews.org/doi/abs/10.1146/annurev-astro-082708-101816""" | |
print '*** test_ic_energy ***' | |
ic = InverseCompton() | |
E_e = 1 # in TeV | |
E = 1e-12 * ic.energy(1e12 * E_e) # in TeV | |
E_expected = (E_e / 11.) ** 2 | |
ratio = E / E_expected | |
print('E_e = %g, E = %g, E_expected = %g, ratio = %g' % | |
(E_e, E, E_expected, ratio)) | |
def test_sy_energy(self): | |
"""Check roughly against peak positions read off Fig. 2a in Hinton et al. | |
Results should be good to within a factor of 2.""" | |
print '*** test_sy_energy ***' | |
sy = Synchrotron(B=3) | |
for E_e, E_expected in zip([1e12, 100e12], [2e-1, 2e3]): # in eV | |
E = sy.energy(E_e) # in eV | |
ratio = E / E_expected | |
print('E_e = %g, E = %g, E_expected = %g, ratio = %g' % | |
(E_e, E, E_expected, ratio)) | |
def test_flux_ratio(self): | |
"""Check that f_IC / f_Sy is as given in Hoppe thesis (1.19). | |
Note that the precision in (1.19) is only one digit!""" | |
print '*** test_flux_ratio ***' | |
for B in [1, sqrt(10), 10]: | |
electron_spec = Model(norm=1, gamma=2, Ecut=1e15) | |
ic_spec = InverseCompton(electron_spec=electron_spec) | |
sy_spec = Synchrotron(B=B, electron_spec=electron_spec) | |
st = SpectrumTester(electron_spec, sy_spec, ic_spec) | |
st.test_flux_ratio(electron_spec.Ecut) | |
def test_sy_cooling_time(self): | |
"""Test synchrotron cooling time""" | |
print '*** test_sy_cooling_time ***' | |
spec = Synchrotron(B=1) | |
spec.info() | |
tau = spec.cooling_time(1e12) | |
tau_expected = 1.3e7 | |
ratio = tau / tau_expected | |
print('tau = %g, tau_expected = %g, ratio = %g' % | |
(tau, tau_expected, ratio)) | |
def test_ic_cooling_time(self): | |
"""Test IC cooling time""" | |
print '*** test_ic_cooling_time ***' | |
spec = InverseCompton(U=1) | |
spec.info() | |
tau = spec.cooling_time(E_e=1e12) | |
tau_expected = 3.1e5 | |
ratio = tau / tau_expected | |
print('tau = %g, tau_expected = %g, ratio = %g' % | |
(tau, tau_expected, ratio)) | |
def test_combined_cooling_time(self, B=10, U_ph=0.26, | |
E=1e7*electron_rest_energy): | |
"""Test Sy + IC combined cooling time.""" | |
print '*** test_combined_cooling_time ***' | |
#tau = Cooling(B, U_ph).time(E) | |
tau = CombinedSyIC(B, U_ph).cooling_time(E) | |
tau_expected = 24.5e3 | |
ratio = tau / tau_expected | |
print('tau = %g, tau_expected = %g, ratio = %g' % | |
(tau, tau_expected, ratio)) | |
def plot_aharonian_electrons(L=1e37, T=1e4, Emin=1e8, Emax=1e16): | |
"""Plot Figure 6.19 in Aharonian's book""" | |
E = logspace(log10(Emin), log10(Emax), 100) | |
#def __init__(self, L=1e37, T=1e4, U_ph=0.26, B=100, | |
# gamma=2, Ecut=1e15, Emin=1, Emax=1e20): | |
for B in [0.1, 1, 3, 10, 30, 100]: | |
electron_spec = Electrons(L=L, T=T, U_ph=0, B=B, Emin=Emin, Emax=Emax) | |
print 'B = %g, Ebreak = %g' % (B, electron_spec.Ebreak) | |
N = electron_spec(E, 2) | |
plot(E, N, label='B = %s' % B) | |
xlabel('E_e (eV)') | |
ylabel('E^2 N(E) (erg)') | |
loglog() | |
#ylim(1e42, 1e48) | |
legend(loc='best') | |
def plot_aharonian_photons(T=1e4, B=5, Ecut=1e15, | |
L=1e37, d=1, Emin=1e-7, Emax=1e16): | |
"""Plot (something similar to) Figure 6.21 in Aharonian's book""" | |
E = logspace(log10(Emin), log10(Emax) + 3, 200) | |
#for B, color in zip([3, 10, 100], ['black', 'red', 'blue']): | |
#for Ecut, color in zip([1e15, 1e14, 1e13], ['black', 'red', 'blue']): | |
for T, color in zip([1e6, 1e4, 1e3], ['black', 'red', 'blue']): | |
electron_spec = Electrons(T=T, B=B, Ecut=Ecut) | |
ic_spec = InverseCompton(electron_spec=electron_spec.spectrum) | |
sy_spec = Synchrotron(B=B, electron_spec=electron_spec.spectrum) | |
F_ic = L_to_F(ic_spec(E, 2), d) | |
F_sy = L_to_F(sy_spec(E, 2), d) | |
plot(E, F_ic, E, F_sy, color=color, label='%g' % Ecut) | |
# Debugging | |
print('B = %s, F_ic = %s, F_sy = %s' % | |
(B, F_ic.max(), F_sy.max())) | |
st = SpectrumTester(electron_spec, sy_spec, ic_spec) | |
st.test_flux_ratio(electron_spec.Ecut) | |
st.print_basic_info() | |
xlabel('E (eV)') | |
ylabel('E F(E) (erg cm^-2 s^-1)') | |
loglog() | |
ylim(1e-14, 1e-7) | |
def plot_rainbow(B=10, Ecut=1e15, | |
energy=1e12, d=1, Emin=1e3, Emax=1e17): | |
"""Plot PWN evolution as a function of time | |
(either spectrum for each time or functions of time. | |
@todo Efficiencies are ridiculous. Needs debugging!!! | |
- Start by checking that for continuous injection the | |
radiated energy equals the injected energy above the | |
breaking energy. | |
- Then look at the electron spectra!""" | |
from matplotlib.cm import rainbow | |
E = logspace(log10(Emin), log10(Emax), 100) | |
pwn = MattanaPWN() | |
times = logspace(0, 6, 1000) | |
fluxes = np.empty_like(times) | |
W = pwn.n_c(times, energy) | |
L = pwn.Edot(times) | |
for ii, T in enumerate(times): | |
#print pwn.t_c(energy) | |
#print pwn.E(T) | |
#print pwn.n_c(T, energy) | |
#print pwn.n_c(T, energy) / pwn.E(T) | |
# @todo Make PWNElectrons class that follows | |
# time-dependent injection energy | |
electron_spec = Electrons(L=W[ii]/T, T=T, B=B, Ecut=Ecut) | |
ic_spec = InverseCompton(electron_spec=electron_spec.spectrum) | |
F_ic = L_to_F(ic_spec(E, 2), d) | |
# Must scale to 0 .. 1 ourselves | |
color_value = log10(T) / 6 | |
color = rainbow(color_value) | |
#color = colorConverter.to_rgb(str(color_value)) | |
#plot(E, F_ic, color=color) | |
#fluxes[ii] = ic_spec.spectrum.total_energy() | |
#fluxes[ii] = ic_spec.spectrum.total_energy(1e12, 1e13) | |
fluxes[ii] = electron_spec.spectrum.total_energy(1e12, 1e13) | |
#fluxes[ii] = ic_spec(1e12, 2) | |
#xlabel('E (eV)') | |
#ylabel('E F(E) (erg cm^-2 s^-1)') | |
#loglog() | |
#ylim(1e-14, 1e-7) | |
y = W | |
figure(); title('W'); plot(times, y); loglog(); ylim(y.max()/1e5, y.max()) | |
y = fluxes | |
figure(); title('fluxes'); plot(times, y); loglog(); ylim(y.max()/1e5, y.max()) | |
y = L | |
figure(); title('L'); plot(times, y); loglog(); ylim(y.max()/1e5, y.max()) | |
y = fluxes / L | |
figure(); title('efficiency'); plot(times, y); loglog(); ylim(y.max()/1e5, y.max()) | |
loglog() | |
#colorbar() | |
def plot_hinton_spectra(total_energy=1e48, gamma=2, B=3, d=1, Emin=1e-7, Emax=1e16): | |
"""Plot Fig 2c in Hinton et al. | |
@todo: I don't understand why the spectrum with the lower cutoff | |
doesn't have a much larger flux at a given energy. | |
It is very slightly larger, but I don't know why this leads | |
to the same total energy. | |
How does Emin affect the normalization!??? | |
""" | |
E = logspace(log10(Emin), log10(Emax), 200) | |
for Ecut, color in zip([1e12, 100e12], ['red', 'blue']): | |
electron_spec = Model(norm=1, gamma=gamma, gamma2=gamma, Ecut=Ecut) | |
#electron_spec.norm = 1e34 | |
electron_spec.set_total_energy(erg_to_eV * total_energy) | |
ic_spec = InverseCompton(electron_spec=electron_spec) | |
ic_spec_dust = InverseCompton(U=0.26, E_T=0.02, electron_spec=electron_spec) | |
sy_spec = Synchrotron(B=B, electron_spec=electron_spec) | |
F_el = L_to_F(electron_spec(E, 2), d) | |
F_ic = L_to_F(ic_spec(E, 2), d) | |
F_ic_dust = L_to_F(ic_spec_dust(E, 2), d) | |
F_sy = L_to_F(sy_spec(E, 2), d) | |
figure(1) | |
plot(E, F_ic, E, F_sy, E, F_ic_dust, label='Ecut = %g' % Ecut, color=color) | |
figure(2) | |
plot(E, F_el, label='Ecut = %g' % Ecut, color=color) | |
# Debugging | |
#print('B = %s, L_ic = %s, F_ic = %s, L_sy = %s, F_sy = %s' % | |
# (B, L_ic.max(), F_ic.max(), L_sy.max(), F_sy.max())) | |
#print sy_spec(1e-3) | |
#st = SpectrumTester(electron_spec, sy_spec, ic_spec) | |
#st.test_flux_ratio() | |
#print electron_spec, electron_spec.total_energy() | |
print ic_spec.spectrum | |
print ic_spec_dust.spectrum | |
figure(1) | |
xlabel('E (eV)') | |
ylabel('E F(E) (erg cm^-2)') | |
loglog() | |
ylim(1e-14, 1e-10) | |
figure(2) | |
xlabel('E (eV)') | |
ylabel('E F(E) (erg cm^-2)') | |
loglog() | |
ylim(1e2) | |
def plot_mattana_evolution(E = 1e7 * electron_rest_energy): | |
pwn = MattanaPWN() | |
t_c = pwn.t_c(E) | |
print 1e-3 * t_c | |
# @todo Mattana have 25 kyr, I have 22.1 kyr. Why??? | |
t = logspace(1, 6, 100) | |
n = pwn.n(t, E) | |
Edot = pwn.Edot(t) | |
n_u = pwn.n_u(t) | |
n_c = pwn.n_c(t, E) | |
n_c_analytic = pwn.n_c_analytic(t, E) | |
n_c_limit = pwn.n_c_limit(t, E) | |
Edot /= Edot.max() | |
n_u /= n.max() | |
n_c /= n.max() | |
n_c_analytic /= n.max() | |
n_c_limit /= n.max() | |
n /= n.max() | |
efficiency = n / Edot | |
efficiency /= efficiency.max() | |
#plot(t, n_u, t, n_c, t, n_c_analytic, t, n) | |
#plot(t, n_c_limit, t, n) | |
#plot(t, n, t, Edot, t, efficiency) | |
print n_c | |
plot(t, n_c) | |
loglog() | |
xlabel('Time (years)') | |
ylabel('Number of particles') | |
#ylim(1e-6, 1e3) | |
#ylim(1e-10, 1e10) | |
vlines(t_c, ylim()[0], ylim()[1]) | |
if __name__ == '__main__': | |
from pylab import * | |
#plot_aharonian_electrons() | |
#plot_aharonian_photons() | |
#plot_rainbow() | |
plot_hinton_spectra() | |
#plot_mattana_evolution() | |
# Tests | |
#ECPL(gamma=3, Eref=2).test_integral(Emin=1e-2, Emax=1e2, power=1) | |
#Model(gamma2=3, Ebreak=2, Ecut=10).test_integral(1e-2, 1e2) | |
#Model(gamma2=3, Ebreak=2, Ecut=10).test_integral(1e-2, 1e2, 1) | |
#RelationsTests() | |
#MattanaPWN().test() | |
#show() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment