subroutine aptscad (np, au, av, aw, bu, bv, bw, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTSCAD
c
c call aptscad (np, au, av, aw, bu, bv, bw, nerr)
c
c Version: aptscad Updated 1990 May 7 16:00.
c aptscad Originated 1990 February 27 10:20.
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 from a uniform distribution in a plane perpendicular to the
c vector a = (au, av, aw) in 3-D space.
c Flag nerr indicates any input error (np not positive).
c
c History: 1990 May 7. Deleted arguments tolu, tolv, tolw from call
c to aptscat.
c
c Calls: aptscat, aptvxun
c
c
c Input: np, au, av, aw.
c
c Output: bu, bv, bw, nerr.
c
c Glossary:
c
c au,av,aw Input The u, v and w components of the vector normal to
c the plane in which the vectors "b" are to be.
c
c bu,bv,bw Output The u, v, w components of a unit vector representing a
c direction chosen randomly from a uniform
c distribution in 3-D space, in the plane with normal
c vector "a". Coordinates u, v and w may
c 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 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---- Component u of random unit vector.
common /laptscad/ cu (64)
c---- Component v of random unit vector.
common /laptscad/ cv (64)
c---- Component w of random unit vector.
common /laptscad/ cw (64)
c---- Component u of vector "a".
common /laptscad/ du (64)
c---- Component v of vector "a".
common /laptscad/ dv (64)
c---- Component w of vector "a".
common /laptscad/ dw (64)
c---- Index in unit vector array.
common /laptscad/ n
c---- First index of subset of data.
common /laptscad/ n1
c---- Last index of subset of data.
common /laptscad/ n2
c---- Size of current subset of data.
common /laptscad/ ns
c---- Magnitude of a vector.
common /laptscad/ vlen (64)
cbugc***DEBUG begins.
cbug common /laptscad/ avgu, avgv, avgw, devu, devv, devw
cbug common /laptscad/ avgu2, avgv2, avgw2, sumsqs
cbug 9901 format (/ 'aptscad finding random directions in a plane.' /
cbug & i3,' au,av,aw=',1p3e22.14)
cbug write (3, 9901) np, au, av, aw
cbugc***DEBUG ends.
c.... initialize.
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)
c.... Loop over the data subsets.
110 ns = n2 - n1 + 1
c.... Find directions sampled uniformly in 3-D space.
call aptscat (ns, cu, cv, cw, nerr)
c.... Shift the unit vectors to the plane perpendicular to vector "a".
c---- Loop over subset of data.
do 120 n = 1, ns
du(n) = au
dv(n) = av
dw(n) = aw
c---- End of loop over subset of data.
120 continue
call aptvxun (cu, cv, cw, du, dv, dw, ns, 0.0,
& bu(n1), bv(n1), bw(n1), vlen, nerr)
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 (/ 'aptscad results:')
cbug 9903 format (i3,' bu,bv,bw=',1p3e22.14 /
cbug & ' sumsq= ',1pe22.14)
cbug 9904 format (/ 'aptscad 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 + 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, 9903) n, bu(n), bv(n), bw(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 aptscad. (+1 line.)
end
UCRL-WEB-209832