The following is a simple 4D simulation where cosmic rays are emitted by a source at a specified spatial position at a specified time-point. A cosmic rays is detected if it arrives at the observer position within a specified time window.
Note: In CRPropa, time is always expressed in terms of redshift $z$, whereas positions are always expressed in terms of comoving coordinates as Cartesian 3-vectors.
The simulation setup is that of a 3D simulation with a few addition:
SourceRedshift
, SourceUniformRedshift
or SourceRedshiftEvolution
.FutureRedshift
implements adiabatic energy loss and updates the redshift. In contrast to Redshift
it allows particles to be propagated into the future $z < 0$ which enables faster convergence for finite observation windows.ObserverRedshiftWindow
specifies a time window $z_\rm{min} < z < z_\rm{max}$ in which particles are detected if they hit the observer. Note that this can also be done after the simulation by cutting on the redshifts at observation. For this we also output the current redshift at observation.Due to the additional time dimension, particles are detected much less often. In order to increase the otherwhise horrible simulation efficiency, a PeriodicBox
is defined: Particles that leave this simulation volume, enter again from the opposite side and their source position is moved accordingly.
As a result the periodic boundaries keep the particles close to the observer and therefore increase the chance of detection. A careful setup is required however:
initTurbulence
as long as the simulation volume coincides with (multiples of) the magnetic field grid.In the example below, a single source is defined. For specifying multiple identical discrete sources SourceMultiplePositions
can be used. Multiple non-identical sources can be added to a SourceList
. For continous source distributions SourceUniformSphere
, SourceUniformSphere
, SourceUniformBox
and SourceUniformCylinder
can be used. SourceDensityGrid
allows to specify a source distribution via a 3D grid.
In [ ]:
from crpropa import *
# set up random turbulent field
vgrid = VectorGrid(Vector3d(0), 512, 30 * kpc)
initTurbulence(vgrid, 8 * nG, 60 * kpc, 800 * kpc, -11. / 3, 42)
Bfield = MagneticFieldGrid(vgrid)
# simulation setup
sim = ModuleList()
sim.add(PropagationCK(Bfield))
sim.add(FutureRedshift())
sim.add(PhotoPionProduction(CMB))
sim.add(PhotoPionProduction(IRB))
sim.add(PhotoDisintegration(CMB))
sim.add(PhotoDisintegration(IRB))
sim.add(ElectronPairProduction(CMB))
sim.add(ElectronPairProduction(IRB))
sim.add(NuclearDecay())
sim.add(MinimumEnergy(1 * EeV))
sim.add(MinimumRedshift(-0.05))
# periodic boundaries
extent = 512 * 30 * kpc # size of the magnetic field grid
sim.add(PeriodicBox(Vector3d(-extent), Vector3d(2 * extent)))
# define the observer
obs = Observer()
obs.add(ObserverSmallSphere(Vector3d(0, 0, 0), 0.5 * Mpc))
obs.add(ObserverRedshiftWindow(-0.05, 0.05))
output = TextOutput('output.txt', Output.Event3D)
output.enable(output.RedshiftColumn)
obs.onDetection(output)
sim.add(obs)
# define the source(s)
source = Source()
source.add(SourcePosition(Vector3d(10, 0, 0) * Mpc))
source.add(SourceIsotropicEmission())
source.add(SourceParticleType(nucleusId(1, 1)))
source.add(SourcePowerLawSpectrum(1 * EeV, 200 * EeV, -1))
source.add(SourceRedshiftEvolution(1.5, 0.001, 3))
# run simulation
sim.setShowProgress(True)
sim.run(source, 10000)