4
$\begingroup$

I'm looking to add an absorption component to this Monte Carlo analysis, essentially I want there to be a 1% chance that a photon gets absorbed at each event as it travels through the medium, once the photon is absorbed the simulation should stop and a new photon brought into the system. The absorption option should come in to play before the scattering angle, I'm just not sure how to write the code to stop further iterations after the photon has been absorbed.

Here is the current code:

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
  ] 

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

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

getTrajectory[startV_, steps_] :=  
  Accumulate@NestList[stepVector, startV, steps];
Graphics3D[
 ParallelTable[
  Line@getTrajectory[{0, 0, -1}, 100],
  {20}
  ]
 ]
$\endgroup$
  • $\begingroup$ Something like Replace your NestList with NestWhileList with a test like RandomReal[]>0.01&. $\endgroup$ – Ymareth Mar 31 '15 at 10:44
3
$\begingroup$

I'll try it as an answer though it's short...

getTrajectory[startV_, steps_] := 
  Accumulate@NestWhileList[stepVector, startV, RandomReal[] >= 1/steps &];
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.