pwn_spectra.py

"""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()