Spectrum from continuous sources

In this example the spectrum from continuous sources is calculated in a model with zero magnetic field. The script can be found in the PARSEC code under ./examples/spectrumContinuousSources.py .

#26.08.11, Tobias Winchen
#PARSEC example - Calculates the spectrum for continuous sources

import parsec
from numpy import linspace, zeros
from matplotlib import pyplot as mpl

#Define Model
attenuationLengthData = parsec.AttenuationLengthData_Berezinsky2006()
energyLossModel = parsec.ContinuousEnergyLoss(attenuationLengthData)
propagationModel = parsec.RandomWalkPropagation(energyLossModel)
propagationModel.setBField(0,1)
spectralIndex = -2.7
maximumEnergy = 1000
sourceModel = parsec.PowerLawSourceModel(energyLossModel, propagationModel,
     spectralIndex, maximumEnergy)

#Borders of the energy bins used in th calculation
energyBins = 10**linspace(18.0, 21,100)
flux = zeros(len(energyBins)-1)
distances = linspace(0.01, 1000, 1000)

for i in range(len(energyBins)-1):
  for d in distances:
    contribution = sourceModel.getRelativeSourceContribution(energyBins[i]/1E18,
        energyBins[i+1]/1E18, d) * d**2
    #sourceModel returns neg. value if no contribution
    if contribution > 0:
      flux[i]+= contribution
    else:
      break

# calculate differential flux and energy bin centers
energies = zeros(len(flux))
for i in range(len(flux)):
  dE = energyBins[i+1]-energyBins[i]
  flux[i]/=dE
  energies[i] = energyBins[i] + dE/2.

#plot data
f = mpl.figure()
mpl.loglog()
mpl.plot(energies, flux*energies**-spectralIndex,'r')
mpl.xlabel('Energy [eV]')
mpl.ylabel('Flux  $\cdot$ E$^{%.2f}$ [a.u.]' % (-spectralIndex))
mpl.show()