subroutine aptscar (kph, ph1, ph2, np, au, av, aw, pm,
& bu, bv, bw, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTSCAR
c
c call aptscar (kph, ph1, ph2, np, au, av, aw, pm,
c & bu, bv, bw, nerr)
c
c Version: aptscar Updated 1992 June 3 18:00.
c aptscar Originated 1990 August 15 11: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 b = (bu, bv, bw)
c at angles between ph1 and ph2 from the direction of the unit
c vector a = (au, av, aw), in either a uniform distribution
c (pm = 0.0) or a cosine**pm distribution centered on "a".
c The units of ph1 and ph2 are determined by the value of kph.
c Flag nerr indicates any input error.
c
c Note: For no restrictions ph1 and ph2, use subroutine aptscat or
c aptscap. For ph1 = ph2, use subroutine aptscan.
c
c History: 1992 June 3 18:00. Moved rotm calculation out of loop 170.
c
c Input: kph, ph1, ph2, np, au, av, aw, pm.
c
c Output: bu, bv, bw, nerr.
c
c Glossary:
c
c au,av,aw Input The u, v and w components of the unit vector in the
c direction of the center of a cosine**pm distribution.
c
c bu,bv,bw Output The u, v, w components of a unit vector representing a
c direction chosen randomly from a cosine**pm
c distribution in 3-D space, centered in the direction
c of unit vector "a". Coordinates u, v and w may
c be any 3 orthogonal coordinates. Size np.
c
c kph Input Specifies the meaning of ph1, ph2:
c -1 to ignore input, use ph1 = 0, ph2 = 180 degrees.
c 0 if ph1 and ph2 are the limiting angles in degrees.
c 1 if ph1 and ph2 are the limiting angles in radians.
c 2 if ph1 and ph2 are cosines of the limiting angles.
c
c nerr Output Indicates an input error, if not 0.
c 1 if np is not positive.
c 2 if kph is not in the range from -1 to 2.
c 3 if pm is -1.0 or less.
c 4 if pm is not zero, and ph1 and/or ph2 exceeds
c 90 degrees or pi / 2 radians, or has a negative
c cosine.
c
c np Input Size of arrays bu, bv, bw.
c
c ph1 Input Specifies the first limiting 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 If pm is not zero, must be between 0 and 90 degrees.
c
c ph2 Input Specifies the second limiting 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 If pm is not zero, must be between 0 and 90 degrees.
c
c pm Input Power used for the cosine**pm spatial distribution
c from which unit vector "b" is to be chosen.
c Zero for a uniform spatial distribution.
c Must be greater than -1.0.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c---- Component u of random unit vector.
dimension bu (1)
c---- Component v of random unit vector.
dimension bv (1)
c---- Component w of random unit vector.
dimension bw (1)
c.... Local variables.
c---- Random value between cosp1 and cosp2.
common /laptscar/ cosp
c---- Function cosph1**(1.0 + pm).
common /laptscar/ cosp1
c---- Function cosph2**(1.0 + pm).
common /laptscar/ cosp2
c---- Cosine of angle phi.
common /laptscar/ cosph
c---- Cosine of first limiting angle.
common /laptscar/ cosph1
c---- Cosine of second limiting angle.
common /laptscar/ cosph2
c---- Cosine of angle theta.
common /laptscar/ costh
c---- Component u of unit vector "c".
common /laptscar/ cu (64)
c---- Component v of unit vector "c".
common /laptscar/ cv (64)
c---- Component w of unit vector "c".
common /laptscar/ cw (64)
c---- Function 1.0 / (1.0 + aw + fuz).
common /laptscar/ faw
c---- A very small number.
common /laptscar/ fuz
c---- Index in unit vector array.
common /laptscar/ n
c---- First index of subset of data.
common /laptscar/ n1
c---- Last index of subset of data.
common /laptscar/ n2
c---- Index in external array.
common /laptscar/ nn
c---- Size of current subset of data.
common /laptscar/ ns
c---- Numerical constant pi.
common /laptscar/ pi
c---- Random number (0.0 to 1.0).
common /laptscar/ ranfp1 (64)
c---- Random number (0.0 to 1.0).
common /laptscar/ ranfp2 (64)
c---- Rotation matrix, (0,0,1) to "a".
common /laptscar/ rotm (3,3)
c---- Sine of angle phi.
common /laptscar/ sinph
c---- Sine of angle of "a" from w axis.
common /laptscar/ sinr
c---- Function sinr + fuz.
common /laptscar/ sinr1
c---- Sine of angle theta.
common /laptscar/ sinth
c---- Angle from v axis in vw plane.
common /laptscar/ theta
cbugc***DEBUG begins.
cbug common /laptscar/ avgu, avgv, avgw, devu, devv, devw
cbug common /laptscar/ avgu2, avgv2, avgw2, sumsqs
cbug 9901 format (/ 'aptscar finding cos**pm limited random directions.' /
cbug & ' np=',i3,' pm=',1pe22.14 /
cbug & ' au,av,aw=',1p3e22.14 /
cbug & ' kph=',i2,' ph1,ph2=',1p2e22.14)
cbug write (3, 9901) np, pm, au, av, aw, kph, ph1, ph2
cbugc***DEBUG ends.
c.... initialize.
c---- A very small number.
fuz = 1.e-99
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
if ((kph .lt. -1) .or. (kph .gt. 2)) then
nerr = 2
go to 210
endif
if (pm .le. -1.0) then
nerr = 3
go to 210
endif
c.... Find the cosines of the limiting angles.
c---- No limiting angles.
if (kph .eq. -1) then
if (pm .eq. 0.0) then
cosph1 = 1.0
cosph2 = -1.0
else
cosph1 = 1.0
cosph2 = 0.0
endif
c---- Units of ph1, ph2 are degrees.
elseif (kph .eq. 0) then
cosph1 = cos (ph1 * pi / 180.0)
cosph2 = cos (ph2 * pi / 180.0)
c---- Units of ph1, ph2 are radians.
elseif (kph .eq. 1) then
cosph1 = cos (ph1)
cosph2 = cos (ph2)
c---- Units of ph1, ph2 are cosines.
elseif (kph .eq. 2) then
cosph1 = ph1
cosph2 = ph2
endif
c.... Truncate small values of cosph1, cosph2 to zero.
if (abs (cosph1) .lt. 1.e-11) then
cosph1 = 0.0
endif
if (abs (cosph2) .lt. 1.e-11) then
cosph2 = 0.0
endif
cbugc***DEBUG begins.
cbug 9902 format (/ ' cosph1,cosph2= ',1p2e22.14)
cbug write ( 3, 9902) cosph1, cosph2
cbugc***DEBUG ends.
c---- Non-uniform spatial distribution.
if (pm .ne. 0.0) then
if ((cosph1 .lt. 0.0) .or. (cosph2 .lt. 0.0)) then
nerr = 4
go to 210
endif
cosp1 = cosph1**(1.0 + pm)
cosp2 = cosph2**(1.0 + pm)
c---- Tested pm.
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.... 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---- Distribution is uniform.
if (pm .eq. 0.0) then
c.... Randomly sample the unit vectors "c", with the w component sampled
c.... from a uniform spatial distribution, between the angles ph1 and ph2.
c---- Loop over vectors.
do 140 n = 1, ns
cosph = cosph1 + ranfp1(n) * (cosph2 - cosph1)
sinph = sqrt (1.0 - cosph**2)
theta = 2.0 * pi * ranfp2(n)
costh = cos (theta)
sinth = sin (theta)
cu(n) = costh * sinph
cv(n) = sinth * sinph
cw(n) = cosph
c---- End of loop over vectors.
140 continue
c---- Distribution is directed toward "a".
else
c.... Randomly sample the unit vectors "c", with the w component sampled
c.... from a cos**pm distribution, between the angles ph1 and ph2.
c---- Loop over vectors.
do 150 n = 1, ns
cosp = cosp1 + ranfp1(n) * (cosp2 - cosp1)
cosph = cosp**(1.0 / (1.0 + pm))
sinph = sqrt (1.0 - cosph**2)
theta = 2.0 * pi * ranfp2(n)
costh = cos (theta)
sinth = sin (theta)
cu(n) = costh * sinph
cv(n) = sinth * sinph
cw(n) = cosph
c---- End of loop over vectors.
150 continue
c---- Tested pm.
endif
c.... Rotate all the unit vectors "c" by the angle required to align the
c.... vector (0., 0., 1.) with the unit vector a = (au, av, aw).
c---- No rotation needed.
if (aw .eq. 1.0) then
c---- Loop over vectors.
do 160 n = 1, ns
nn = n + n1 - 1
bu(nn) = cu(n)
bv(nn) = cv(n)
bw(nn) = cw(n)
c---- End of loop over vectors.
160 continue
c---- Must rotate.
else
c.... Find rotation matrix.
sinr = sqrt (au**2 + av**2)
sinr1 = sinr + fuz
faw = (1.0 - aw) / sinr1**2
rotm(1,1) = av**2 * faw + aw
rotm(1,2) = -au * av * faw
rotm(1,3) = au * sinr / sinr1
rotm(2,1) = -au * av * faw
rotm(2,2) = au**2 * faw + aw
rotm(2,3) = av * sinr / sinr1
rotm(3,1) = -au * sinr / sinr1
rotm(3,2) = -av * sinr / sinr1
rotm(3,3) = aw
c---- Loop over vectors.
do 170 n = 1, ns
nn = n + n1 - 1
bu(nn) = cu(n) * rotm(1,1) + cv(n) * rotm(1,2) +
& cw(n) * rotm(1,3)
bv(nn) = cu(n) * rotm(2,1) + cv(n) * rotm(2,2) +
& cw(n) * rotm(2,3)
bw(nn) = cu(n) * rotm(3,1) + cv(n) * rotm(3,2) +
& cw(n) * rotm(3,3)
c---- End of loop over vectors.
170 continue
c---- Tested aw.
endif
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 9903 format (/ 'aptscar results:')
cbug 9904 format (i4,' bu,bv,bw=',1p3e22.14 /
cbug & ' sumsq= ',1pe22.14)
cbug 9905 format (/ 'aptscar mean and deviation:' /
cbug & ' avg(u,v,w)= ',1p3e22.14 /
cbug & ' dev(u,v,w)= ',1p3e22.14)
cbug
cbug write ( 3, 9903)
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 + bu(n)
cbug avgv = avgv + bv(n)
cbug avgw = avgw + bw(n)
cbug avgu2 = avgu2 + bu(n)**2
cbug avgv2 = avgv2 + bv(n)**2
cbug avgw2 = avgw2 + bw(n)**2
cbug sumsqs = bu(n)**2 + bv(n)**2 + bw(n)**2
cbug write ( 3, 9904) n, bu(n), bv(n), bw(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 (avgu2 - avgu**2)
cbug devv = sqrt (avgv2 - avgv**2)
cbug devw = sqrt (avgw2 - avgw**2)
cbug
cbug write ( 3, 9905) avgu, avgv, avgw, devu, devv, devw
cbugc***DEBUG ends.
210 continue
cbugc***DEBUG begins.
cbug 9906 format (/ ' nerr=',i3)
cbug write ( 3, 9906) nerr
cbugc***DEBUG ends.
return
c.... End of subroutine aptscar. (+1 line.)
end
UCRL-WEB-209832