[project @ 1996-07-25 21:02:03 by partain]
[nofib.git] / real / gamteb / TransPort.hs
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