subroutine aptscma (kph, ph, np, au, av, aw, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTSCMA
c
c call aptscma (kph, ph, np, au, av, aw, nerr)
c
c Version: aptscma Updated 1991 August 5 14:00.
c aptscma Originated 1991 August 5 14:00.
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 cone at an angle phi from the u axis in 3-D space.
c The angle phi is specified by "ph", based on the option "kph".
c Flag nerr indicates any input error (np not positive).
c
c Input: kph, ph, np.
c
c Output: au, av, aw, nerr.
c
c Glossary:
c
c au Output The u component of a unit vector sampled randomly from
c all directions in 3-D space that are a specified
c angle (see kph, ph) from the u axis. Size np.
c
c av Output The v component of a unit vector sampled randomly from
c all directions in 3-D space that are a specified
c angle (see kph, ph) from the u axis. Size np.
c
c aw Output The w component of a unit vector sampled randomly from
c all directions in 3-D space that are a specified
c angle (see kph, ph) from the u axis. Size np.
c
c kph Input Specifies the meaning of ph:
c 0 if ph is the scattering angle in degrees.
c 1 if ph is the scattering angle in radians.
c 2 if ph is the cosine of the scattering angle.
c
c nerr Output Indicates an input error, if not 0.
c 1 if np is not positive.
c 2 if kph is not 0, 1 or 2.
c 3 if kph = 2, and ph is not between -1.0 and 1.0.
c
c np Input Size of arrays au, av and aw.
c
c ph Input Specifies the scattering angle in degrees, from
c 0.0 to 180.0 (kph = 0), in radians, from 0.0 to pi
c (kph = 1), or its cosine, from -1.0 to 1.0 (kph = 2).
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 scattering angle.
common /laptscma/ cosph
c---- Cosine of azimuth angle.
common /laptscma/ costh
c---- Index in local array.
common /laptscma/ n
c---- First index of subset of data.
common /laptscma/ n1
c---- Last index of subset of data.
common /laptscma/ n2
c---- Index in external array.
common /laptscma/ nn
c---- Size of current subset of data.
common /laptscma/ ns
c---- Scattering angle in radians.
common /laptscma/ phi
c---- Mathematical constant, pi.
common /laptscma/ pi
c---- Random number (0.0 to 1.0).
common /laptscma/ ranfp (64)
c---- Sine of scattering angle.
common /laptscma/ sinph
c---- Sine of azimuth angle.
common /laptscma/ sinth
c---- Azimuth angle around the u axis.
common /laptscma/ theta
cbugc***DEBUG begins.
cbug common /laptscma/ avgu, avgv, avgw, devu, devv, devw
cbug common /laptscma/ avgu2, avgv2, avgw2, sumsqs
cbug 9901 format (/ 'aptscma scattering by PH around u axis:' /
cbug & ' kph=',i2,' ph=',1pe22.14)
cbug write (3, 9901) kph, ph
cbugc***DEBUG ends.
c.... initialize.
c---- Mathematical constant, pi.
pi = 3.14159265358979323
nerr = 0
c.... Test for input errors.
if (np .le. 0) then
nerr = 1
go to 210
endif
if ((kph .lt. 0) .or. (kph .gt. 2)) then
nerr = 2
go to 210
endif
c.... Set up the indices of the first subset of data.
n1 = 1
n2 = min (np, 64)
c.... Loop over the data subsets.
110 ns = n2 - n1 + 1
c.... Find the scattering angle cosine and sine.
c---- Unit of ph is degrees.
if (kph .eq. 0) then
phi = ph * pi / 180.0
cosph = cos (phi)
sinph = sin (phi)
c---- Unit of ph is radians.
elseif (kph .eq. 1) then
phi = ph
cosph = cos (phi)
sinph = sin (phi)
c---- Unit of ph is cosine of angle.
elseif (kph .eq. 2) then
cosph = ph
if (abs (cosph) .gt. 1.0) then
nerr = 3
go to 210
endif
sinph = sqrt (1.0 - cosph**2)
endif
c.... Generate the needed random numbers.
c---- Loop over subset of data.
do 150 n = 1, ns
ranfp(n) = ranf( )
c---- End of loop over subset of data.
150 continue
c.... Sample randomly on a cone at angle phi from u axis.
c---- Loop over vectors.
do 160 n = 1, ns
theta = 2.0 * pi * ranfp(n)
costh = cos (theta)
sinth = sin (theta)
nn = n + n1 - 1
au(nn) = cosph
av(nn) = costh * sinph
aw(nn) = sinth * sinph
c---- End of loop over vectors.
160 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 (/ 'aptscma results:')
cbug 9903 format (i3,' au,av,aw=',1p3e22.14 /
cbug & ' sumsq= ',1pe22.14)
cbug 9904 format (/ 'aptscma 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 180 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 180 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 (amax1 (0.0, avgu2 - avgu**2))
cbug devv = sqrt (amax1 (0.0, avgv2 - avgv**2))
cbug devw = sqrt (amax1 (0.0, avgw2 - avgw**2))
cbug
cbug write ( 3, 9904) avgu, avgv, avgw, devu, devv, devw
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptscma. (+1 line.)
end
UCRL-WEB-209832