1 -- 2 -- Patricia Fasel 3 -- Los Alamos National Laboratory 4 -- 1990 August 5 -- 6 module TransPort (transPort) where 7 8 import GamtebType 9 import Consts 10 import Utils 11 import Distance 12 import RoulSplit 13 import PhotoElec 14 import Compton 15 import Pair 16 (=:) a b = (a,b) 17 18 -- transport a particle 19 20 transPort :: Particle -> Probability -> ([Result], [Stat]) 21 transPort p prob = 22 if (dColl < dSurf) 23 then (let 24 pos' = transPos pos dir dColl 25 p' = Part pos' dir w e eIndx cell seed 26 doCompton = (r1 < (pComp / (pTot-pPhot))) 27 (res, stat) = collision p' prob doCompton 28 in 29 (res, [nc=:1]++stat)) 30 else (let 31 pos' = transPos pos dir dSurf 32 p' = Part pos' dir w e eIndx cell seed 33 (res, stat) = noCollision p' prob surf 34 in 35 (res, [nnc=:1]++stat)) 36 where 37 (Part pos dir w e eIndx cell seed) = p 38 (pComp,pPair,pPhot,pTot) = prob 39 (r, r1) = genRand seed 40 (dSurf, surf) = distSurf pos dir 41 dColl = -((log r)/ pTot) 42 43 44 -- no collision inside cylinder 45 46 noCollision :: Particle -> Probability -> Int -> ([Result], [Stat]) 47 noCollision p@(Part pos dir w e eIndx cell seed) prob surf = 48 case surf of 49 1 -> ([(scatter, eIndx)=:w], [ns=:1]) 50 2 -> ([(escape, eIndx)=:w], [ne=:1]) 51 4 -> ([(transit, eIndx)=:w], [nt=:1]) 52 3 -> -- cross internal surface 53 -- particle will split, die in russian roulette, or continue 54 -- cell = [1..] causes roulet or split to alternate 55 if (cell == 1) 56 then (let 57 (p1, p2) = split p 58 (r1, s1) = transPort p1 prob 59 (r2, s2) = transPort p2 prob 60 in 61 (r1++r2, [nsp=:1]++s1++s2)) 62 else 63 (let 64 (p', stat', roulKill) = roulet p 65 in 66 if (roulKill) 67 then ([], stat') 68 else (let 69 (res, stat) = transPort p' prob 70 in 71 (res, stat++stat'))) 72 73 74 -- collision is in cylinder, do collision physics 75 76 collision :: Particle -> Probability -> Bool -> ([Result], [Stat]) 77 collision p prob doCompton = 78 if (wgtKill) 79 then ([], [nwk=:1]) 80 else 81 if (doCompton) 82 then -- compton scattering 83 (let 84 (p'', prob', comptonCut) = compton p' 85 in 86 if (comptonCut) 87 then ([], [nek=:1]) 88 else transPort p'' prob') 89 else -- pair production 90 (let 91 (p'', prob', pairCut) = pair p' 92 in 93 if (pairCut) 94 then ([], [nek=:1]) 95 else transPort p'' prob') 96 where 97 (Part pos dir w e eIndx cell seed) = p 98 (pComp,pPair,pPhot,pTot) = prob 99 (p', absorb, wgtKill) = photoElec p prob 100 101 102 -- translate a particle position 103 104 transPos :: Point -> Point -> Value -> Point 105 transPos (x,y,z) (u,v,w) dist = 106 (x',y',z') 107 where 108 x' = x + u*dist 109 y' = y + v*dist 110 z' = z + w*dist