This notebook explains how to use cubic results of 3D MHD models on a uniform grid in CRPropa.
The fields need to be supplied in a raw binary file that contains only single floats, arranged as follows: Starting with the cell values (Bx,By,Bz for magnetic field or rho for density) at the origin of the box, the code continues to read along z, then y and finally x.
On https://crpropa.desy.de/ under "Additional resources" you can find a number of MHD models used with CRPropa in the literature.
In [1]:
from crpropa import *
## settings for MHD model (must be set according to model)
filename_bfield = "clues_primordial.dat" ## filename of the magnetic field
gridOrigin = Vector3d(0,0,0) ## origin of the 3D data, preferably at boxOrigin
gridSize = 1024 ## size of uniform grid in data points
size = 249.827*Mpc ## physical edgelength of volume in Mpc
b_factor = 1. ## global renormalizatino factor for the field
## settings of simulation
boxOrigin = Vector3d( 0, 0, 0,) ## origin of the full box of the simulation
boxSize = Vector3d( size, size, size ) ## end of the full box of the simulation
## settings for computation
minStep = 10.*kpc ## minimum length of single step of calculation
maxStep = 4.*Mpc ## maximum length of single step of calculation
tolerance = 1e-2 ## tolerance for error in iterative calculation of propagation step
spacing = size/(gridSize) ## resolution, physical size of single cell
m = ModuleList()
## instead of computing propagation without Lorentz deflection via
# m.add(SimplePropagation(minStep,maxStep))
## initiate grid to hold field values
vgrid = VectorGrid( gridOrigin, gridSize, spacing )
## load values to the grid
loadGrid( vgrid, filename_bfield, b_factor )
## use grid as magnetic field
bField = MagneticFieldGrid( vgrid )
## add propagation module to the simulation to activate deflection in supplied field
m.add(PropagationCK( bField, tolerance, minStep, maxStep))
#m.add(DeflectionCK( bField, tolerance, minStep, maxStep)) ## this was used in older versions of CRPropa
to make use of periodicity of the provided data grid, use
In [2]:
m.add( PeriodicBox( boxOrigin, boxSize ) )
to not follow particles forever, use
In [3]:
m.add( MaximumTrajectoryLength( 400*Mpc ) )
In [4]:
source = Source()
source.add( SourceUniformBox( boxOrigin, boxSize ))
In [65]:
filename_density = "mass-density_clues.dat" ## filename of the density field
source = Source()
## initialize grid to hold field values
mgrid = ScalarGrid( gridOrigin, gridSize, spacing )
## load values to grid
loadGrid( mgrid, filename_density )
## add source module to simulation
source.add( SourceDensityGrid( mgrid ) )
In [67]:
filename_halos = 'clues_halos.dat'
# read data from file
data = np.loadtxt(filename_halos, unpack=True, skiprows=39)
sX = data[0]
sY = data[1]
sZ = data[2]
mass_halo = data[5]
## find only those mass halos inside the provided volume (see Hackstein et al. 2018 for more details)
Xdown= sX >= 0.25
Xup= sX <= 0.75
Ydown= sY >= 0.25
Yup= sY <= 0.75
Zdown= sZ >= 0.25
Zup= sZ <= 0.75
insider= Xdown*Xup*Ydown*Yup*Zdown*Zup
## transform relative positions to physical positions within given grid
sX = (sX[insider]-0.25)*2*size
sY = (sY[insider]-0.25)*2*size
sZ = (sZ[insider]-0.25)*2*size
## collect all sources in the multiple sources container
smp = SourceMultiplePositions()
for i in range(0,len(sX)):
pos = Vector3d( sX[i], sY[i], sZ[i] )
smp.add( pos, 1. )
## add collected sources
source = Source()
source.add( smp )
additional source properties
In [5]:
## use isotropic emission from all sources
source.add( SourceIsotropicEmission() )
## set particle type to be injected
A, Z = 1, 1 # proton
source.add( SourceParticleType( nucleusId(A,Z) ) )
## set injected energy spectrum
Emin, Emax = 1*EeV, 1000*EeV
specIndex = -1
source.add( SourcePowerLawSpectrum( Emin, Emax, specIndex ) )
In [6]:
filename_output = 'output.txt'
filename_output = 'data/output_MW.txt'
obsPosition = Vector3d(0.5*size,0.5*size,0.5*size) # position of observer, MW is in center of constrained simulations
obsSize = 800*kpc ## physical size of observer sphere
## initialize observer that registers particles that enter into sphere of given size around its position
obs = Observer()
obs.add( ObserverSmallSphere( obsPosition, obsSize ) )
## write registered particles to output file
obs.onDetection( TextOutput( filename_output ) )
## choose to not further follow particles paths once detected
obs.setDeactivateOnDetection(True)
## add observer to module list
m.add(obs)
finally run the simulation by
In [7]:
N = 1000
m.showModules() ## optional, see summary of loaded modules
m.setShowProgress(True) ## optional, see progress during runtime
m.run(source, N, True) ## perform simulation with N particles injected from source