forked from lynch829/Monte-Carlo-Radiation-Transport-code
-
Notifications
You must be signed in to change notification settings - Fork 0
/
stokes.f
153 lines (124 loc) · 3.27 KB
/
stokes.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
subroutine stokes(nxp,nyp,nzp,sint,cost,sinp,cosp,phi,
+ hgg,g2,pi,twopi,iseed)
implicit none
integer iseed
real*8 nxp,nyp,nzp,sint,cost,sinp,cosp,phi
real*8 hgg,g2,pi,twopi
real*8 costp,sintp,phip
real*8 bmu,b,ri1,ri3,cosi3,sini3,cosb2,sinbt,sini2,bott,cosdph
real*8 cosi2,sini1,cosi1
real ran2
c***** isotropic scattering if g = 0.0 ******************************
if(hgg.eq.0.0) then
cost=2.*ran2(iseed)-1.
sint=(1.-cost*cost)
if(sint.le.0.)then
sint=0.
else
sint=sqrt(sint)
endif
phi=twopi*ran2(iseed)
sinp=sin(phi)
cosp=cos(phi)
nxp=sint*cosp
nyp=sint*sinp
nzp=cost
else
c***** heyney greenstein scattering ********************************
costp=cost
sintp=sint
phip=phi
bmu=((1.+g2)-((1.-g2)/(1.-hgg+2.*hgg*ran2(iseed)))**2)/(2.*hgg)
cosb2=bmu**2
b=cosb2-1.
if(abs(bmu).gt.1.) then
if(bmu.gt.1.) then
bmu=1.
cosb2=1.
b=0.
else
bmu=-1.
cosb2=1.
b=0.
end if
end if
sinbt=sqrt(1.-cosb2)
ri1=twopi*ran2(iseed)
if(ri1.gt.pi) then
ri3=twopi-ri1
cosi3=cos(ri3)
sini3=sin(ri3)
if(bmu.eq.1.) then
goto 100
else
if(bmu.eq.-1.) then
goto 100
end if
end if
cost=costp*bmu+sintp*sinbt*cosi3
if(abs(cost).lt.1.) then
sint=abs(sqrt(1.-cost*cost))
sini2=sini3*sintp/sint
bott=sint*sinbt
cosi2=costp/bott-cost*bmu/bott
else
sint=0.
sini2=0.
if(cost.ge.1.) cosi2=-1.
if(cost.le.-1.) cosi2=1.
end if
cosdph=-cosi2*cosi3+sini2*sini3*bmu
if(abs(cosdph).gt.1.) then
if(cosdph.gt.1.) then
cosdph=1.
else
cosdph=-1.
end if
end if
phi=phip+acos(cosdph)
if(phi.gt.twopi) phi=phi-twopi
if(phi.lt.0.) phi=phi+twopi
c elseif(ri1.le.pi) then
else
cosi1=cos(ri1)
sini1=sin(ri1)
if(bmu.eq.1.) then
goto 100
else
if(bmu.eq.-1.) then
goto 100
end if
end if
cost=costp*bmu+sintp*sinbt*cosi1
if(abs(cost).lt.1.) then
sint=abs(sqrt(1.-cost*cost))
sini2=sini1*sintp/sint
bott=sint*sinbt
cosi2=costp/bott-cost*bmu/bott
else
sint=0.
sini2=0.
if(cost.ge.1.) cosi2=-1.
if(cost.le.-1.) cosi2=1.
end if
cosdph=-cosi1*cosi2+sini1*sini2*bmu
if(abs(cosdph).gt.1.) then
if(cosdph.gt.1.) then
cosdph=1.
else
cosdph=-1.
end if
end if
phi=phip-acos(cosdph)
if(phi.gt.twopi) phi=phi-twopi
if(phi.lt.0.) phi=phi+twopi
end if
cosp=cos(phi)
sinp=sin(phi)
nxp=sint*cosp
nyp=sint*sinp
nzp=cost
end if
100 continue
return
end