subroutine aptscat (np, au, av, aw, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTSCAT
c
c call aptscat (np, au, av, aw, nerr)
c
c Version: aptscat Updated 1990 May 7 10:40.
c aptscat Originated 1990 January 10 10:40.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To randomly sample np unit direction vectors a = (au, av, aw)
c from a uniform distribution in 3-D space.
c Flag nerr indicates any input error (np not positive).
c
c History: 1990 April 30 16:00. Switched from rejection to trig method.
c 1990 May 7 10::40. Removed arguments tolx, toly, tolz.
c This option may be exercized by calling aptvlim or aptvtol.
c
c Input: np.
c
c Output: au, av, aw, nerr.
c
c Glossary:
c
c au,av,aw Output The u, v, w components of a unit vector representing a
c random direction in 3-d space. Coordinates u, v and
c w may be any 3 orthogonal coordinates. Size np.
c
c nerr Output Indicates an input error, if not 0.
c 1 if np is not positive.
c
c np Input Size of arrays.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c---- Component u of random unit vector.
dimension au (1)
c---- Component v of random unit vector.
dimension av (1)
c---- Component w of random unit vector.
dimension aw (1)
c.... Local variables.
c---- Cosine of angle phi.
common /laptscat/ cosph
c---- Cosine of angle theta.
common /laptscat/ costh
c---- Index in unit vector array.
common /laptscat/ n
c---- First index of subset of data.
common /laptscat/ n1
c---- Last index of subset of data.
common /laptscat/ n2
c---- Index in external array.
common /laptscat/ nn
c---- Size of current subset of data.
common /laptscat/ ns
c---- Numerical constant pi.
common /laptscat/ pi
c---- Random number (0.0 to 1.0).
common /laptscat/ ranfp1 (64)
c---- Random number (0.0 to 1.0).
common /laptscat/ ranfp2 (64)
c---- Sine of angle phi.
common /laptscat/ sinph
c---- Sine of angle theta.
common /laptscat/ sinth
c---- Angle from v axis in vw plane.
common /laptscat/ theta
cbugc***DEBUG begins.
cbug common /laptscat/ avgu, avgv, avgw, devu, devv, devw
cbug common /laptscat/ avgu2, avgv2, avgw2, sumsqs
cbug 9901 format (/ 'aptscat finding isotropic random directions.',
cbug & ' np=',i3)
cbug write (3, 9901) np
cbugc***DEBUG ends.
c.... initialize.
c++++ Dimensionless. 18 digits.
pi = 3.14159265358979323
nerr = 0
c.... Test for input errors.
if (np .le. 0) then
nerr = 1
go to 210
endif
c.... Set up the indices of the first subset of data.
n1 = 1
n2 = min (np, 64)
110 ns = n2 - n1 + 1
c.... Generate the needed random numbers.
c---- Loop over subset of data.
do 120 n = 1, ns
ranfp1(n) = ranf( )
c---- End of loop over subset of data.
120 continue
c---- Loop over subset of data.
do 130 n = 1, ns
ranfp2(n) = ranf( )
c---- End of loop over subset of data.
130 continue
c.... Find the randomly oriented unit vectors.
c---- Loop over vectors.
do 140 n = 1, ns
nn = n + n1 - 1
cosph = 1.0 - 2.0 * ranfp1(n)
sinph = sqrt (1.0 - cosph**2)
theta = 2.0 * pi * ranfp2(n)
costh = cos (theta)
sinth = sin (theta)
au(nn) = cosph
av(nn) = costh * sinph
aw(nn) = sinth * sinph
c---- End of loop over vectors.
140 continue
c.... See if all data subsets are done.
c---- Do another subset of data.
if (n2 .lt. np) then
n1 = n2 + 1
n2 = min (np, n1 + 63)
go to 110
endif
cbugc***DEBUG begins.
cbug
cbug 9902 format (/ 'aptscat results:')
cbug 9903 format (i3,' au,av,aw=',1p3e22.14 /
cbug & ' sumsq= ',1pe22.14)
cbug 9904 format (/ 'aptscat mean and deviation:' /
cbug & ' avg(u,v,w)=',1p3e22.14 /
cbug & ' dev(u,v,w)=',1p3e22.14)
cbug
cbug write ( 3, 9902)
cbug
cbug avgu = 0.0
cbug avgv = 0.0
cbug avgw = 0.0
cbug avgu2 = 0.0
cbug avgv2 = 0.0
cbug avgw2 = 0.0
cbug
cbug do 160 n = 1, np
cbug avgu = avgu + au(n)
cbug avgv = avgv + av(n)
cbug avgw = avgw + aw(n)
cbug avgu2 = avgu2 + au(n)**2
cbug avgv2 = avgv2 + av(n)**2
cbug avgw2 = avgw2 + aw(n)**2
cbug sumsqs = au(n)**2 + av(n)**2 + aw(n)**2
cbug write ( 3, 9903) n, au(n), av(n), aw(n), sumsqs
cbug 160 continue
cbug
cbug avgu = avgu / np
cbug avgv = avgv / np
cbug avgw = avgw / np
cbug avgu2 = avgu2 / np
cbug avgv2 = avgv2 / np
cbug avgw2 = avgw2 / np
cbug devu = sqrt (avgu2 - avgu**2)
cbug devv = sqrt (avgv2 - avgv**2)
cbug devw = sqrt (avgw2 - avgw**2)
cbug
cbug write ( 3, 9904) avgu, avgv, avgw, devu, devv, devw
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptscat. (+1 line.)
end
UCRL-WEB-209832