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)