# How to set up a simple Monte Carlo simulation?

there is a 1% chance that there will be a scattering event, if there is a scattering event the direction of scattering will be determined by a random sample of the custom[a_] function and if there is no scattering event i would like the photon to continue its previous trajectory.

"custom probability function of scattering angle"

custom[a_] :=
ProbabilityDistribution[
0.000569772*(350*Exp[-(x)^2/2*0.26^2] + 7.12*Exp[-0.105*x] + 0.0007),
{x, 0, 180}];


"chance of a scattering event occurring"

RandomChoice[{0.01, 0.99} -> {1, 0}


"current monte carlo walk simulation without if statement"

randomWalk[n_] := Accumulate[Prepend[RandomVariate[custom[a_], n], 0]]
ListLinePlot[Table[randomWalk, {5}]]

• (1) what is it exactly you want to calculate? Looking at your simulation results you wish to see a bunch of trajectories of photons moving though a medium. (2) what does "a" mean in custom[a]? – Dr. Wolfgang Hintze Mar 4 '15 at 11:03
• Is the problem 2 or 3 dimensional? Are your step lengths uniform? – djp Mar 4 '15 at 11:14
• If i want to modelling the scattering from neutrino using the considerations given above, what fundamental modifications i have to make. For instance, if i have a cross section for a given process, how i can implement her in the monte carlo calculations?? – user38192 Mar 2 '16 at 0:38
• If you have a new question, please ask it by clicking the Ask Question button. Include a link to this question if it helps provide context. - From Review – Michael E2 Mar 2 '16 at 1:01

My interpretation is: the photon will be scattered by an angle $\alpha$ (given by getScatterAngle), and the deviation will occur with equal probability in every direction. (For example, a photon initially going along the $z$ axis will be rotated by an angle $\alpha$ about a randomly chosen axis that lies in the $x,y$ plane.)

When I've written Monte Carlo simulations I always use this general scheme:

### Functions to support a "step" of the simulation

getScatterAngle :=
(\[Pi]/180) RandomVariate[
ProbabilityDistribution[
0.000569772*(350*Exp[-(x)^2/2*0.26^2] + 7.12*Exp[-0.105*x] +
0.0007), {x, 0, 180}]
]

getRotationAxis[v : {x_, y_, z_}] :=
Module[{perpVector},
If[x > 0,
perpVector = {y, -x, 0},
perpVector = {0, z, -y}];
RotationMatrix[RandomReal[2 \[Pi]], v].perpVector
] (* End Module *)

getScatteredVector[v : {x_, y_, z_}] :=
RotationMatrix[getScatterAngle, getRotationAxis[v]].v


### Function to perform a step of the simulation

stepVector[v_] :=
RandomChoice[{0.1, 0.9} -> {getScatteredVector[v], v}]


### Function to complete one simulation

getTrajectory[startV_, steps_] :=
Accumulate@NestList[stepVector, startV, steps];


### Function to complete multiple simulation chains

To run it all with 20 photons and 100 steps each:

Graphics3D[
ParallelTable[
Line@getTrajectory[{0, 0, 1}, 100],
{20}
]
] If I wrote the algorithm I wouldn't "scatter or advance", but rather I would normally advance by an exponentially-distributed distance and then scatter. It makes little difference here.

However, the reason the code is so slow is that it's very expensive to calculate scatter angles like this (38 ms per call on my computer). If you were sampling from a simple distribution (like a Normal distribution), things would be a whole lot quicker.

I am very suspicious of your

Exp[-(x)^2/2*0.26^2]

Is the 0.26^2 meant to be a numerator or denominator? It's a numerator at present.

• Thanks so much for your help, given that all the data is compressed into one variable how would i go about getting the dimensions of the plot and putting axis labels and whatever else on top? Sorry for being a pain, i've been trying to figure it out for a while now with no success – AW1991 Mar 5 '15 at 2:11
• I think this is a research problem now - you need to figure out "exactly what is my problem?" The diagram is schematic as the dimensions are meaningless (the step length is "1" -- 1 what? metre? millimetre?). But have a read of the help for Graphics3D. – djp Mar 5 '15 at 2:29
• @djp This is very interesting.. how will it be the code if using the GPU for paralization? – Jose Enrique Calderon Jul 12 '19 at 4:18

Let me give a model solution which can easily be adapted.

2 dimensions

Consider a photon which moves in the x-y-plane, starting at time t = 0 in the origin {0,0} and moving towards the positive x-axis. At each tick of the clock, corresponding to a constant distance 1 travelled by the photon, the photon will experience a scattering event which leads to a deviation of the original direction by an angle "a". The angle "a" is a random variable with a certain given distibution dist.

We ask for the random trajectory of the photon after a given number n of scattering events. We would like to see several of those trajectories in one picture.

For simplicity, I assume that the angle distribution is a symmetric NormaDistribution[0,s] with a given width s.

dist[s_] = NormalDistribution[0, s];


The scattering angle is then

a[s_] := RandomVariate[dist[s]]


A sample of angles is

Array[a &, 10];


We need to cumulate these angeles because the new one gives the deviation with respect to the current one:

Accumulate[%];


Now the x-y-coordinates of each subtrajectory of one tick of the clock are

{Cos[#], Sin[#]} & /@ %;


This again needs to be cumulated in order to get the coordinates of the trajectory:

Accumulate[%];


Now we have one trajectory an plot it

ListLinePlot[%]
(* 150304 plot_1_trajectory.jpg *) Putting everything together gives us a trajectory by

tr[s_, n_] :=
Accumulate[{Cos[#], Sin[#]} & /@ (Accumulate[Array[a[s] &, n]])]


Plotting a bunch of k = 11 trajectories, each of length n = 100 with a width s of the normal Distribution

With[{s = 0.15, n = 100, k = 11}, ListLinePlot[Table[tr[s, n], {k}]]]
(* 150304 plot_k_trajectories.jpg *) Varying s from small values (s = .01, prefers forward directions) to larger ones (s = 1, gives strong scattering) shows interesting outcomes.

3 dimensions

dist[s_] =
NormalDistribution[0, s]; (* distribution of scattering angle *)

a[s_] := RandomVariate[dist[s]]; (* scattering angle *)

f := 2 \[Pi] Random []; (* rotation angle *)

k = 30; (* length of a trajectory *)
m = 8; (* number of trajectories *)
s = 0.3; (* width of distribution of scattering angle *)

(* run starts here *)
tbtrk = Table[af = Accumulate[Array[{a[s], f} &, k]];
{Hue[j/m],
trk = Line @@ {Accumulate[
Prepend[Table[{Sin[#1] Cos[#2], Sin[#1] Sin[#2],
Cos[#1]} &@(Sequence @@ af[[i]]), {i, 1,
Length[af]}], {0, 0, 0}]]}}, {j, 0, m}];

d = 5; (* size of the display box *)
Graphics3D[{{PointSize[0.02], Point[{0, 0, 0}]}, tbtrk},
PlotRange -> {{-d, d}, {-d, d}, {-5, k}}, Boxed -> True]
(* 150304_plot _ 3D_k _trajectories.jpg *) • I just saw that 6 minutes before me, the user djp has given a (seemingly) complete 3D simulation result. – Dr. Wolfgang Hintze Mar 4 '15 at 12:16
• I am getting no output in the 2D code. Is there a typo in rv[s]? – Jose Enrique Calderon Jul 12 '19 at 12:30
• @ Jose Enrique Calderon You are right, there's a typo in the 2D code: rv[s] should read a[s]. – Dr. Wolfgang Hintze Jul 12 '19 at 16:59