<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 3D trajectories in a turbulent field
The following simulation tracks a single UHE nucleus and its secondary nucleons/nuclei through a turbulent magnetic field.
First we create a random realization of a turbulent field with a Kolmogorov power spectrum on 60-800 kpc lengthscales and an RMS field strength of 8 nG.
The field is stored on a $256^3$ grid with 30 kpc grid spacing, and thus has an extent of $(256 \cdot 30 \rm{kpc})^3$.
The field is by default periodically repeated in space to cover an arbitrary volume.
The chosen grid size consumes only very little memory. For practical purposes a larger grid is advised in order to represent more variations of turbulent modes, provide a larger turbulent range, or a higher resolution.
Step1:  Saving and loading fields
In addition to creating random turbulent fields, we can also load and save custom magnetic field grids.
As input and output we currently support binary files in single precision and ASCII files.
Step2:  Running the simulation
Now that we have our magnetic field ready we can fire up our simulation and hope that something visually interesting is going to happen.
Step3:  (Optional) Plotting
We plot the trajectory of our oxygen-16 nucleus. To distinguish between secondary nuclei the following colors are used
<ASSISTANT_TASK:>
Python Code:
from crpropa import *
randomSeed = 42
turbSpectrum = SimpleTurbulenceSpectrum(Brms=8*nG, lMin = 60*kpc, lMax=800*kpc, sIndex=5./3.)
gridprops = GridProperties(Vector3d(0), 256, 30*kpc)
BField = SimpleGridTurbulence(turbSpectrum, gridprops, randomSeed)
# print some properties of our field
print('Lc = {:.1f} kpc'.format(BField.getCorrelationLength() / kpc))  # correlation length
print('sqrt(<B^2>) = {:.1f} nG'.format(BField.getBrms() / nG))   # RMS
print('<|B|> = {:.1f} nG'.format(BField.getMeanFieldStrength() / nG))  # mean
print('B(10 Mpc, 0, 0) =', BField.getField(Vector3d(10,0,0) * Mpc) / nG, 'nG')
Explanation: 3D trajectories in a turbulent field
The following simulation tracks a single UHE nucleus and its secondary nucleons/nuclei through a turbulent magnetic field.
First we create a random realization of a turbulent field with a Kolmogorov power spectrum on 60-800 kpc lengthscales and an RMS field strength of 8 nG.
The field is stored on a $256^3$ grid with 30 kpc grid spacing, and thus has an extent of $(256 \cdot 30 \rm{kpc})^3$.
The field is by default periodically repeated in space to cover an arbitrary volume.
The chosen grid size consumes only very little memory. For practical purposes a larger grid is advised in order to represent more variations of turbulent modes, provide a larger turbulent range, or a higher resolution.
End of explanation
# save the field
# format: (Bx, By, Bz)(x, y, z) with z changing the quickest.
#dumpGrid(BField.getGrid(), 'myfield.dat')       # binary, single precision
#dumpGridToTxt(Bfield.getGrid(), 'myfield.txt')  # ASCII
# load your own field
#vgrid=Grid3f(gridprops)
#loadGrid(vgrid, 'myfield.dat')
#loadGridFromTxt(vgrid, 'myfield.txt')
Explanation: Saving and loading fields
In addition to creating random turbulent fields, we can also load and save custom magnetic field grids.
As input and output we currently support binary files in single precision and ASCII files.
End of explanation
sim = ModuleList()
sim.add(PropagationCK(BField))
sim.add(PhotoPionProduction(CMB()))
sim.add(PhotoPionProduction(IRB_Kneiske04()))
sim.add(PhotoDisintegration(CMB()))
sim.add(PhotoDisintegration(IRB_Kneiske04()))
sim.add(ElectronPairProduction(CMB()))
sim.add(ElectronPairProduction(IRB_Kneiske04()))
sim.add(NuclearDecay())
sim.add(MaximumTrajectoryLength(25 * Mpc))
output = TextOutput('trajectory.txt', Output.Trajectory3D) 
sim.add(output)
x = Vector3d(0,0,0)  # position
p = Vector3d(1,1,0)  # direction
c = Candidate(nucleusId(16, 8), 100 * EeV, x, p)
sim.run(c, True)
Explanation: Running the simulation
Now that we have our magnetic field ready we can fire up our simulation and hope that something visually interesting is going to happen.
End of explanation
%matplotlib inline
from pylab import *
from mpl_toolkits.mplot3d import axes3d
output.close()
data = genfromtxt('trajectory.txt', names=True)
# trajectory points
x, y, z = data['X'], data['Y'], data['Z']
# translate particle ID to charge number
Z = [chargeNumber(int(Id)) for Id in data['ID'].astype(int)]
# translate the charge number to color and size
# --> protons are blue, Helium is green, everthing else is red
colorDict = {0:'k', 1:'b', 2:'g', 3:'r', 4:'r', 5:'r', 6:'r', 7:'r', 8:'r'}
sizeDict = {0:4, 1:4, 2:8, 3:10, 4:10, 5:10, 6:10, 7:10, 8:10}
colors = [colorDict[z] for z in Z]
sizes  = [sizeDict[z] for z in Z]
fig = plt.figure(figsize=(12, 5))#plt.figaspect(0.5))
ax = fig.gca(projection='3d')# , aspect='equal'
ax.scatter(x,y,z+6, 'o', s=sizes, color=colors)
ax.set_xlabel('x / Mpc', fontsize=18)
ax.set_ylabel('y / Mpc', fontsize=18)
ax.set_zlabel('z / Mpc', fontsize=18)
ax.set_xlim((-1, 16))
ax.set_ylim((-1, 16))
ax.set_zlim((-1, 16))
ax.xaxis.set_ticks((0, 5, 10, 15))
ax.yaxis.set_ticks((0, 5, 10, 15))
ax.zaxis.set_ticks((0, 5, 10, 15))
show()
Explanation: (Optional) Plotting
We plot the trajectory of our oxygen-16 nucleus. To distinguish between secondary nuclei the following colors are used: protons are blue, alpha particles are green, everthing heavier is red.
End of explanation
<END_TASK>