home /
academic notes /
The Monte Carlo simulation playbook
A practical introduction to Monte Carlo radiation/particle transport code in fortran
Basic overview
Step one: emit a particle from a location \((x, y, z\)) in a direction $$\begin{align} \textbf{n} &= \begin{bmatrix} \sin\theta \cos\phi \\ \sin\theta \sin\phi \\ \cos\theta \end{bmatrix} \end{align}$$ Step two: move particle to new position at a distance dependent on the particle energy and optical properties/composition of the medium.Step three: determine the nature of the particle's interaction with the medium at that point;
- Scattering event - change direction, and possibly energy of particle, then calculate new step.
- Absorption in to medium - absorb energy into medium, emit new particle
- Fission (only relevant for neutronics) - produces secondary neutrons and photons
Monte Carlo packets
The most common particle in MC simulation is the Monte Carlo luminosity packet. For a total luminosity of \(L\), each packet carries an energy \(E_i = L \Delta t/N\); where \(N\) is the total number of packets and \(\Delta t\) is the physical time over which the simulation is being performed. \(\Delta t\) is what allows us to compute absorbed power/energy, and physical properties of the medium do not change over this timescale.One MC packet represents \(N_\gamma\) real photons, where \(N_\gamma = E_i/h \nu_i\)
For neutronics simulation, one MC packet represents one neutron.
Sampling a random optical depth
If you don't know how to sample a random variable from an underlying distribution function, check out this quick primer on [ random sampling from PDFs ].The probability for a particle to travel a distance \(\tau\) before interaction is \(P(\tau) = e^{-\tau}\). So we sample a random \(\tau\) from this distribution using $$\xi = \int_0^\tau e^{-\tau '} d\tau ' = 1 - e^{-\tau} \Rightarrow \tau = -\log(1-\xi)$$ Seeing as \(1-\xi\) is also a random number distributed uniformly between zero and one, we can simply this to $$\tau = -log(\xi)$$ The physical distance that a packet travels, \(L\), can therefore be derived from the equation $$\tau = \int_0^L n\sigma ds$$ which is easy to solve in the case of a uniform density medium: $$L = \tau ( s_{max}/\tau_{max} )$$
Sampling a random direction
Given that directions in a sphere are determined by two angles \(\theta\) and \(\phi\), it's easy to imagine that you could sample a random solid angle by sampling those two angles from a uniform distribution. However, this causes a bunching up of point at higher latitudes (low \(\theta\)).\(\theta\) is normalised over \([0, \pi]\) and \(\phi\) is normalised over \([0, 2\pi]\), so their respetive PDFs are $$p(\theta) = 1/2 \sin\theta$$ $$p(\phi) = 1/2\pi$$ Using the sampling process as described earlier, we arrive at $$\theta = \cos^{-1}(2\xi - 1)$$ $$\phi = 2\pi \xi$$