subroutine aptbrkp (ksys, iunit, u, du, v, dv, w, dw, fu, fv, fw,
& np, tol, pu, pv, pw, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTBRKP
c
c call aptbrkp (ksys, iunit, u, du, v, dv, w, dw, fu, fv, fw,
c & np, tol, pu, pv, pw, nerr)
c
c Version: aptbrkp Updated 1994 March 3 12:00.
c aptbrkp Originated 1994 March 3 12:00.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To find np points p = (pu, pv, pw) in the inside of the
c volume element defined in the coordinate system ksys
c (0 = rectangular, 1 = cylindrical, 2 = spherical) with angle
c units iunit (0 = degrees, 1 = radians), bounded by the
c coordinate surfaces at u, u + du, v, v + dv, w and w + dw.
c The np points will be at the fractional distances across the
c volume element, fu, fv and fw.
c Coordinates of points "p" less than the estimated error in their
c calculation, based on tol, will be truncated to zero.
c Flag nerr indicates any input error.
c
c Input: ksys, iunit, u, du, v, dv, w, dw, fu, fv, fw, np, tol.
c
c Output: pu, pv, pw, nerr.
c
c Calls: None
c
c Glossary:
c
c fu,fv,fw Input The fractional distance (on a volume basis) of a
c point "p" across the volume element. Size np.
c
c iunit Input Indicates the unit used for angles:
c 0 if angles theta and phi are in degrees.
c 1 if angles theta and phi are in radians.
c Not used if ksys = 0.
c
c ksys Input Indicates the coordinate system type:
c 0 for Cartesian coordinates. u = x, v = y, w = z.
c 1 for cylindrical coordinates. u = radius from z
c axis, v = angle theta in xy plane, counterclockwise
c from x axis, w = z.
c 2 for spherical coordinates. u = radius from origin,
c v = angle theta in xy plane, counterclockwise from
c x axis, w = angle phi from z axis.
c
c nerr Output Indicates an input error, if not 0.
c 1 if ksys is not 0, 1 or 2.
c 2 if iunit is not 0 or 1.
c 3 if np is not positive.
c
c np Input Size of arrays pu, pv, pw, fu, fv, fw.
c Number of points "p" to find. Must be positive.
c
c pu,pv,pw Output The coordinates of a spatial point "p" in the volume
c element defined in coordinate system ksys, with angle
c units iunit, bounded by the coordinate surfaces at
c u = u, u = u + du, v = v, v = bvr + dv, w = w and
c w = dw, at volume fractions fu, fv and fw across the
c volume element from (u, v, w). Values less than the
c estimated error in their calculation, based on tol,
c will be truncated to zero. Size np.
c
c tol Input Numerical tolerance limit.
c On Cray computers, recommend 1.e-5 to 1.e-11.
c
c u, du Input The coordinate and increment of the bounding coordinate
c surfaces of the volume element in the u direction.
c
c v, dv Input The coordinate and increment of the bounding coordinate
c surfaces of the volume element in the v direction.
c
c w, dw Input The coordinate and increment of the bounding coordinate
c surfaces of the volume element in the w direction.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c____ Mathematical constant pi. Dimensionless.
parameter (pi = 3.14159265358979323)
c.... Dimensioned arguments.
c---- Coordinate u of point "p".
dimension pu (1)
c---- Coordinate v of point "p".
dimension pv (1)
c---- Coordinate w of point "p".
dimension pw (1)
c---- Fractional u distance of point "p".
dimension fu (1)
c---- Fractional v distance of point "p".
dimension fv (1)
c---- Fractional w distance of point "p".
dimension fw (1)
c.... Local variables.
c---- Cosine of angle phi.
common /laptbrkp/ cosw
c---- Cosine of minimum angle phi.
common /laptbrkp/ cosw1
c---- Cosine of maximum angle phi.
common /laptbrkp/ cosw2
c---- Angle conversion factor.
common /laptbrkp/ fact
c---- Index in output array.
common /laptbrkp/ n
c---- Estimated error in pu.
common /laptbrkp/ puerr
c---- Estimated error in pv.
common /laptbrkp/ pverr
c---- Estimated error in pw.
common /laptbrkp/ pwerr
c---- Square of minimum radius.
common /laptbrkp/ r1sq
c---- Square of maximum radius.
common /laptbrkp/ r2sq
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptbrkp finding points in a brick. np=',i3,
cbug & ' ksys=',i3,' iunit=',i2 /
cbug & ' tol=',1pe22.14 /
cbug & ' u, du= ',1p2e22.14 /
cbug & ' v, dv= ',1p2e22.14 /
cbug & ' w, dw= ',1p2e22.14 )
cbug
cbug write ( 3, 9901) np, ksys, iunit, tol, u, du, v, dv, w, dw
cbug
cbug 9902 format (i3,' fu,fv,fw=',1p3e22.14 )
cbug
cbug do 110 n = 1, np
cbug write ( 3, 9902) n, fu(n), fv(n), fw(n)
cbug 110 continue
cbugc***DEBUG ends.
c.... initialize.
nerr = 0
c.... Test for input errors.
if ((ksys .lt. 0) .or. (ksys .gt. 2)) then
nerr = 1
go to 210
endif
if ((iunit .ne. 0) .and. (iunit .ne. 1)) then
nerr = 2
go to 210
endif
if (np .le. 0) then
nerr = 3
go to 210
endif
c.... Find the angle conversion factor (degrees to radians).
if (iunit .eq. 0) then
fact = pi / 180.0
else
fact = 1.0
endif
c.... Find points "p" in the brick.
if (ksys .eq. 0) then
c---- Loop over data.
do 150 n = 1, np
pu(n) = u + fu(n) * du
pv(n) = v + fv(n) * dv
pw(n) = w + fw(n) * dw
c---- End of loop over data.
150 continue
elseif (ksys .eq. 1) then
c---- Loop over data.
do 160 n = 1, np
r1sq = u**2
r2sq = (u + du)**2
pu(n) = sqrt (r1sq + fu(n) * (r2sq - r1sq))
pv(n) = v + fv(n) * dv
pw(n) = w + fw(n) * dw
c---- End of loop over data.
160 continue
elseif (ksys .eq. 2) then
c---- Loop over data.
do 170 n = 1, np
r1cb = u**3
r2cb = (u + du)**3
pu(n) = (r1cb + fu(n) * (r2cb - r1cb))**(1.0 / 3.0)
pv(n) = v + fv(n) * dv
cosw1 = cos (fact * w)
cosw2 = cos (fact * (w + dw))
cosw = cosw1 + fw(n) * (cosw2 - cosw1)
pw(n) = acos (cosw) / fact
c---- End of loop over data.
170 continue
c---- Tested ksys.
endif
c.... See if truncation to zero is allowed.
c---- Truncation to zero allowed.
if (tol .gt. 0.0) then
if (ksys .eq. 0) then
c---- Loop over data.
do 180 n = 1, np
puerr = tol * (abs (u) + fu(n) * abs (du))
pverr = tol * (abs (v) + fv(n) * abs (dv))
pwerr = tol * (abs (w) + fw(n) * abs (dw))
if (abs (pu(n)) .lt. puerr) then
pu(n) = 0.0
endif
if (abs (pv(n)) .lt. pverr) then
pv(n) = 0.0
endif
if (abs (pw(n)) .lt. pwerr) then
pw(n) = 0.0
endif
c---- End of loop over data.
180 continue
elseif (ksys .eq. 1) then
c---- Loop over data.
do 190 n = 1, np
pverr = tol * (abs (v) + fv(n) * abs (dv))
pwerr = tol * (abs (w) + fw(n) * abs (dw))
if (abs (pv(n)) .lt. pverr) then
pv(n) = 0.0
endif
if (abs (pw(n)) .lt. pwerr) then
pw(n) = 0.0
endif
c---- End of loop over data.
190 continue
elseif (ksys .eq. 2) then
c---- Loop over data.
do 200 n = 1, np
cosw1 = cos (fact * w)
cosw2 = cos (fact * (w + dw))
cosw = cosw1 + fw(n) * (cosw2 - cosw1)
pverr = tol * (abs (v) + fv(n) * abs (dv))
cwerr = tol * (abs (cosw1) +
& fw(n) * (abs (cosw2) + abs (cosw1)))
if (abs (pv(n)) .lt. pverr) then
pv(n) = 0.0
endif
if (abs (cosw - 1.0) .lt. cwerr) then
pw(n) = 0.0
endif
if (abs (cosw) .lt. cwerr) then
pw(n) = 0.5 * pi / fact
endif
if (abs (cosw + 1.0) .lt. cwerr) then
pw(n) = pi / fact
endif
c---- End of loop over data.
200 continue
c---- Tested ksys.
endif
c---- Tested tol.
endif
cbugc***DEBUG begins.
cbug
cbug 9903 format (/ 'aptbrkp results:')
cbug 9904 format (i3,' pu,pv,pw=',1p3e22.14 )
cbug
cbug write ( 3, 9903)
cbug do 205 n = 1, np
cbug write ( 3, 9904) n, pu(n), pv(n), pw(n)
cbug 205 continue
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptbrkp. (+1 line.)
end
UCRL-WEB-209832