Bielajew A.F. Fundamentals of the Monte Carlo method for neutral and charged particle transport
Подождите немного. Документ загружается.
18.3. SUBROUTINES FOR PARTICLE TRANSPORT AND DEFLECTION 321
18.3 Subroutines for particle transport and deflection
azimuthal0.f
Uniform sampling of azimuthal direction vectors in 2-D.
The direction cosines are determined using the following method:
u =cos(2πr)andv =sin(2πr)
where r is a random number uniform on the range 0 ≤ r<1, or 0 <r≤ 1, or 0 ≤ r ≤ 1.
C23456789|123456789|123456789|123456789|123456789|123456789|123456789|12
subroutine azimuthal(u,v)
implicit none
C Computes unit normalized random transverse direction cosines u,v
C normalized so that u**2 + v**2 = 1
real
* u, ! Output: final x-axis direction cosine
* v ! Output: final y-axis direction cosine
real
* phi, ! Internal: azimuthal angle phi
* twopi, ! Internal: value of 2*pi
* rng ! Internal: a random number uniform on the range [0,1) or
! (0,1] or [0,1].
! (,) => exclusive
! [,] => inclusive
data twopi/6.283185/
save twopi ! Saved for repeated use.
phi = twopi*rng()
u = cos(phi)
v = sin(phi)
return
end
322 CHAPTER 18. CODE LIBRARY
18.3. SUBROUTINES FOR PARTICLE TRANSPORT AND DEFLECTION 323
rotate.f
Rotates ~u −→ ~u
0
according to cos(Θ).
C23456789|123456789|123456789|123456789|123456789|123456789|123456789|12
subroutine rotate(u,v,w,ws)
implicit none
C Rotates the direction vectors (u,v,w) after a scatter by polar cosine ws
C Some notation conventions:
C Theta, Phi refer to the scattering angles in the reference frame where the
C z-axis is aligned with the initial particle direction
C theta, phi refer to the original particle directions with respect to the
C laboratory frame of reference
real
* u, ! Input: initial x-axis direction cosine
* ! Output: final x-axis direction cosine
* v, ! Input: initial y-axis direction cosine
* ! Output: final y-axis direction cosine
* w, ! Input: initial z-axis direction cosine
* ! Output: final z-axis direction cosine
* ws ! Input: cos(Theta)
real
* sinThetas, ! Internal: sin(Theta), scattering angle
* sintheta, ! Internal: sin(theta), initial angle theta
* sinthetai, ! Internal: 1/sin(theta)
* cosphi, ! Internal: cos(phi), initial angle phi
* sinphi, ! Internal: sin(phi), initial angle phi
* us, ! Internal: scattering x-direction cosine,
! sin(Theta)*cos(Phi)
* vs, ! Internal: scattering y-direction cosine
! sin(Theta)*sin(Phi)
* u2v2 ! Internal: u**2 + v**2 = sin(Theta)**2
C Determine the azimuthal direction cosines
324 CHAPTER 18. CODE LIBRARY
call azimuthal(us,vs) ! Unit normalized 2-D transverse direction cosines
! i.e. cos(Phi), sin(Phi)
C Convert 2-D transverse to 3-D transverse direction cosines
sinThetas = sqrt(1 - ws**2)
us = sinThetas*us ! i.e. us = sin(Theta)*cos(Phi)
vs = sinThetas*vs ! i.e. vs = sin(Theta)*sin(Phi)
u2v2 = u**2 + v**2 ! i.e. u2v2 = sin(theta)**2
if (u2v2 .gt. 0e0) then ! General scattering case
sintheta = sqrt(u2v2) ! i.e. sin(theta)
sinthetai = 1/sintheta
cosphi = u*sinthetai ! i.e. cos(phi)
sinphi = v*sinthetai ! i.e. sin(phi)
u = u*ws + us*w*cosphi - vs*sinphi
v = v*ws + us*w*sinphi + vs*cosphi
w = w*ws - us*sintheta
return
end if
C This is the scattering from initial Z-direction case
u=us
v=vs
w=ws
return
end
18.4. SUBROUTINES FOR MODELING INTERACTIONS 325
18.4 Subroutines for modeling interactions
isotropic.f
Samples the z-axis direction cosine in isotropic scattering from the probability distribution:
p(θ)dθ =
sin θ dθ
2
0 ≤ θ ≤ π
where θ is the polar angle. The direct method is employed giving
cos θ =2r − 1
where r is a random number uniform on the range 0 ≤ r<1, or 0 <r≤ 1, or 0 ≤ r ≤ 1.
C23456789|123456789|123456789|123456789|123456789|123456789|123456789|12
subroutine isotropic(costhe)
implicit none
C Computes direction cosine for isotropic scattering
real
* costhe ! Output: cos(Theta)
real
* rng ! Internal: random number
costhe = 2e0*rng() - 1e0
return
end
326 CHAPTER 18. CODE LIBRARY
rutherford.f
Samples the z-axis direction cosine in “Rutherfordian” scattering from the probability dis-
tribution:
p(µ)=
a(2 + a)
2
1
(1 − µ + a)
2
; − 1 ≤ µ ≤ 1,
where µ = cos(θ). The direct method is employed giving
µ =1− 2a
1 − r
a +2r
where r is a random number uniform on the range 0 ≤ r<1, or 0 <r≤ 1, or 0 ≤ r ≤ 1.
Restrictions: a ≥ 0 otherwise the routine will stop execution. a = 0 is treated as a special
case (no scattering).
C23456789|123456789|123456789|123456789|123456789|123456789|123456789|12
subroutine rutherford(costhe,eta)
implicit none
C Rutherford-type scattering
real
* costhe, ! Output: cos(Theta)
* eta ! Input: screening factor
real
* rng ! Internal: random number
if (eta .gt. 0e0) then
costhe = rng()
costhe = 1e0 - 2e0*eta*(1e0 - costhe)/(eta + 2e0*costhe)
return
elseif (eta .eq. 0e0) then
costhe = 1e0
return
else
write(6,*) ’Stopped in subroutine rutherford.’
write(6,*) ’Eta was less than zero.’
stop 1
endif
18.4. SUBROUTINES FOR MODELING INTERACTIONS 327
end
328 CHAPTER 18. CODE LIBRARY
18.5 Subroutines for modeling geometry
xplane.f
Geometry routine for a plane normal to the x-axis.
C23456789|123456789|123456789|123456789|123456789|123456789|123456789|12
subroutine xplane(p0,x0,u,s,inside,hit)
implicit none
C Calculates the distance of a particle to a planar surface normal to the
C x-axis
logical
* inside ! Input: inside = .true. => particle thinks it is inside
* , ! Input: inside = .false. => particle thinks it is outside
* hit ! Output: hit = .true. => particle would hit the surface
! Output: hit = .false. => particle would miss the surface
real
* p0, ! Input: Point at which the plane cuts the x-axis
* x0, ! Input: x-coordinate of the particle
* u , ! Input: x-axis direction cosine of the particle
* s ! Output: Distance to the surface (if hit)
if (
* (inside.and.u.gt.0e0) !headed towards the surface
* .or.
* (.not.inside.and.u.lt.0e0) !headed towards the surface
*)
*then
hit = .true.
s = max(0e0,(p0-x0)/u)
return
else
hit = .false.
return
endif
end
18.5. SUBROUTINES FOR MODELING GEOMETRY 329
yplane.f
Geometry routine for a plane normal to the y-axis.
C23456789|123456789|123456789|123456789|123456789|123456789|123456789|12
subroutine yplane(p0,y0,v,s,inside,hit)
implicit none
C Calculates the distance of a particle to a planar surface normal to the
C y-axis
logical
* inside ! Input: inside = .true. => particle thinks it is inside
* , ! Input: inside = .false. => particle thinks it is outside
* hit ! Output: hit = .true. => particle would hit the surface
! Output: hit = .false. => particle would miss the surface
real
* p0, ! Input: Point at which the plane cuts the y-axis
* y0, ! Input: y-coordinate of the particle
* v , ! Input: y-axis direction cosine of the particle
* s ! Output: Distance to the surface (if hit)
if (
* (inside.and.v.gt.0e0) !headed towards the surface
* .or.
* (.not.inside.and.v.lt.0e0) !headed towards the surface
*)
*then
hit = .true.
s = max(0e0,(p0-y0)/v)
return
else
hit = .false.
return
endif
end
330 CHAPTER 18. CODE LIBRARY
zplane.f
Geometry routine for a plane normal to the z-axis.
C23456789|123456789|123456789|123456789|123456789|123456789|123456789|12
subroutine zplane(p0,z0,w,s,inside,hit)
implicit none
C Calculates the distance of a particle to a planar surface normal to the
C z-axis
logical
* inside ! Input: inside = .true. => particle thinks it is inside
* , ! Input: inside = .false. => particle thinks it is outside
* hit ! Output: hit = .true. => particle would hit the surface
! Output: hit = .false. => particle would miss the surface
real
* p0, ! Input: Point at which the plane cuts the z-axis
* z0, ! Input: z-coordinate of the particle
* w , ! Input: z-axis direction cosine of the particle
* s ! Output: Distance to the surface (if hit)
if (
* (inside.and.w.gt.0e0) !headed towards the surface
* .or.
* (.not.inside.and.w.lt.0e0) !headed towards the surface
*)
*then
hit = .true.
s = max(0e0,(p0-z0)/w)
return
else
hit = .false.
return
endif
end