Bielajew A.F. Fundamentals of the Monte Carlo method for neutral and charged particle transport
Подождите немного. Документ загружается.
18.5. SUBROUTINES FOR MODELING GEOMETRY 331
csphere.f
Geometry routine for a sphere, centered at the origin.
C23456789|123456789|123456789|123456789|123456789|123456789|123456789|12
subroutine csphere(R,x0,y0,z0,u,v,w,s,inside,hit)
implicit none
C Calculates the distance to a sphere centered at (0,0,0)
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
* R, ! Input: Radius of the cylinder
* x0, ! Input: x-coordinate of the particle
* y0, ! Input: y-coordinate of the particle
* z0, ! Input: z-coordinate of the particle
* u , ! Input: x-axis direction cosine of the particle
* v , ! Input: y-axis direction cosine of the particle
* w , ! Input: w-axis direction cosine of the particle
* s ! Output: Distance to the surface (if hit)
real
* A, ! Internal: Quadratic coefficient A
* B, ! Internal: Quadratic coefficient B
* C ! Internal: Quadratic coefficient C
A = 1e0 ! i.e. u**2 + v**2 + w**2 = 1
B = u*x0 + v*y0 + w*z0
C = x0**2 + y0**2 + z0**2 - R**2
call quadric(A,B,C,s,inside,hit) ! Get the generic quadric solution
return
end
332 CHAPTER 18. CODE LIBRARY
ccylz.f
Geometry routine for a circular cylinder, aligned and concentric with the z-axis
C23456789|123456789|123456789|123456789|123456789|123456789|123456789|12
subroutine ccylz(R,x0,y0,u,v,s,inside,hit)
implicit none
C Calculates the distance to a circular cylinder centered and aligned along
C the 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
* R, ! Input: Radius of the cylinder
* x0, ! Input: x-coordinate of the particle
* y0, ! Input: y-coordinate of the particle
* u , ! Input: x-axis direction cosine of the particle
* v , ! Input: y-axis direction cosine of the particle
* s ! Output: Distance to the surface (if hit)
real
* A, ! Internal: Quadratic coefficient A
* B, ! Internal: Quadratic coefficient B
* C ! Internal: Quadratic coefficient C
A = u**2 + v**2
B = u*x0 + v*y0
C = x0**2 + y0**2 - R**2
call quadric(A,B,C,s,inside,hit) ! Get the generic quadric solution
return
end
18.5. SUBROUTINES FOR MODELING GEOMETRY 333
quadric.f
Geometry routine for the general solution of closest point-of-flight positive intersection with
a general quadric surface
C23456789|123456789|123456789|123456789|123456789|123456789|123456789|12
subroutine quadric(A,B,C,s,inside,hit)
implicit none
C Calculates the first positive distance to an arbitrary quadric surface
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
* A, ! Input: Quadratic coefficient A
* B, ! Input: Quadratic coefficient B
* C, ! Input: Quadratic coefficient C
* s ! Output: Distance to the surface (if hit)
real
* Q ! Internal: quadratic coefficient
Q = B**2 - A*C
if (Q .lt. 0e0) then
hit = .false.
return
end if
if (inside) then !inside the surface
if (B .le. 0e0) then !headed away from surface
if (A .gt. 0e0) then !but, surface curving up
hit = .true. !always a hit in this case
s = (sqrt(Q) - B)/A
return
else !surface curving away and headed in
334 CHAPTER 18. CODE LIBRARY
hit = .false. !never a hit in this case
return
end if
else !headed toward the surface
hit = .true. !always a hit in this case
s = max(0e0, -C/(sqrt(Q) + B))
return
end if
end if
!Must be outside the surface
if (B .ge. 0e0) then !headed away
if (A .ge. 0e0) then !surface curves away
hit = .false. !never a hit in this case
return
else !surface curves up
hit = .true. !always a hit in this case
s = -(sqrt(Q) + B)/A
return
end if
else !headed toward the surface
hit = .true. !always a hit in this case
s = max(0e0, C/(sqrt(Q) - B))
return
end if
end
18.6. TEST ROUTINES 335
18.6 Test routines
testiso.f
Tests isotropic scattering
C23456789|123456789|123456789|123456789|123456789|123456789|123456789|12
implicit none
C test isotropic scattering
real
* costhe, ! Output: cos(Theta)
* sinthe, ! Output: sin(Theta)
* cosphi, ! Output: cos(Phi)
* sinphi ! Output: sin(Phi)
integer i
do i = 1,1000
call isotropic(costhe)
sinthe = sqrt(1 - costhe**2)
call azimuthal(cosphi,sinphi)
write(1,*) 0e0, 0e0
write(1,*) costhe, sinthe*cosphi
write(1,*) 0e0,0e0
write(2,*) 0e0, 0e0
write(2,*) costhe, sinthe*sinphi
write(2,*) 0e0,0e0
write(3,*) 0e0, 0e0
write(3,*) sinthe*cosphi, sinthe*sinphi
write(3,*) 0e0,0e0
end do
end
336 CHAPTER 18. CODE LIBRARY
testrng.f
Tests the random number generators
implicit none
integer i
real rng, rn0, rn
i=0
rn0 = rng()
write(6,*) ’ Initial RN = ’,rn0
1 continue
i=i+1
rn = rng()
if (i .le. 50 .and. i .gt. 0) write(6,*) rn
if (rn .eq. rn0) then
write(6,*) ’RNG looped, stopping, i = ’, i
write(6,*) ’Final RN = ’, rn
stop
endif
goto 1
end
Bibliography
[Jam94] F. James. RANLUX: A FORTRAN implementation of the high-quality pseudoran-
dom number generator of L¨uscher. Computer Physics Communications, 79:111 –
114, 1994.
[L
¨
94] M. L¨uscher. A portable high-quality random number generator for lattice field
theory simulations. Computer Physics Communications, 79:100 – 110, 1994.
337
338 BIBLIOGRAPHY