1 -- 2 -- Patricia Fasel 3 -- Los Alamos National Laboratory 4 -- 1990 August 5 -- 6 module Compton (compton) where 7 8 import GamtebType 9 import Consts 10 import Utils 11 12 -- compton scattering 13 14 compton :: Particle -> (Particle, Probability, Bool) 15 compton (Part pos dir w e eIndx cell seed) = 16 if (e' <= ergCut) 17 then (Part pos dir w e' eIndx' cell seed', prob', True) 18 else (Part pos dir' w e' eIndx' cell seed', prob', False) 19 where 20 (seed', r2) = genRand seed 21 (r3, r4) = genRand r2 22 eIn = 1.956917 * e 23 eOut = klein eIn r3 24 angle = 1 + 1/eIn - 1/eOut 25 e' = 0.511008 * eOut 26 (eIndx', prob') = xsectInterp e' 27 dir' = rotas dir angle r4 28 29 30 -- rotate a point through a polar angle whose cosine is c 31 -- and through an azimuthal angle sampled uniformly 32 33 rotas :: Point -> Angle -> Random -> Point 34 rotas (u,v,w) a rn = 35 if (r > 1) 36 then rotas (u,v,w) a rn' 37 else 38 (let 39 r' = sqrt ((1 - a*a) / r) 40 t1' = t1*r' 41 t2' = t2*r' 42 wsq = 1 - w*w 43 s = sqrt wsq 44 u' = u*a + (t1'*u*w - t2'*v) / s 45 v' = v*a + (t1'*v*w - t2'*u) / s 46 w' = w*a - t1'*s 47 in 48 if (wsq < small) 49 then (t1',t2',(w*a)) 50 else (u',v',w') 51 ) 52 where 53 (r1, r2) = genRand rn 54 (rn', r3) = genRand r2 55 t1 = 2*r1 - 1 56 t2 = 2*r3 - 1 57 r = t1*t1 + t2*t2 58 59 60 -- sample from klein-nishina using inverse fit 61 -- e = energy in, units of the rest mass of an electron 62 63 klein :: Energy -> Random -> Energy 64 klein e r = 65 if (e > 1.16666667) 66 then 67 (let 68 a' = 1.65898 + a*(0.62537*a - 1.00796) 69 b' = a'/f 70 in 71 if (r > b') 72 then (let 73 c' = (d-1.20397) / (1-b') 74 x' = 0.3 * exp (c'*(b'- r)) 75 in 76 x'*e) 77 else (let 78 c' = a'/(3.63333 + a*(5.44444*a - 4.66667)) 79 x' = klein1 (r/b') 2.1 c' 1.4 (0.5*a') 80 in 81 x'*e)) 82 else (let 83 x' = klein1 r (3*c') a' (2*c') b' 84 a' = f/(b+c) 85 b' = 0.5*f 86 c' = 1-c 87 in 88 x'*e) 89 where 90 a = 1/e 91 b = 2*e + 1 92 c = 1/b 93 d = log b 94 f = 2*e*(1+e)*c*c + 4*a + (1-2*a*(1+a))*d 95 klein1 x2 x3 x4 x5 x7 = 96 1 + x2*(x2*(2*x7 + x4 - x3 + x2*(x5 - x7 - x4)) - x7)