subroutine aptbrkn (ksys, iunit, u, du, v, dv, w, dw, pu, pv, pw,
& np, tol, nside, fu, fv, fw, dsurf, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTBRKN
c
c call aptbrkn (ksys, iunit, u, du, v, dv, w, dw, pu, pv, pw,
c & np, tol, nside, fu, fv, fw, dsurf, nerr)
c
c Version: aptbrkn Updated 1995 February 23 11:10.
c aptbrkn 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 if np points p = (pu, pv, pw) are 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 fractional distances across the volume element fu, fv and
c fw, and the distances dsurf of each point from the six bounding
c surfaces of the brick will also be returned. Points on the
c surface of the volume element, based on tol, will be considered
c to be inside.
c Flag nerr indicates any input error.
c
c Input: ksys, iunit, u, du, v, dv, w, dw, pu, pv, pw, np, tol.
c
c Output: nside, fu, fv, fw, dsurf, nerr.
c
c Calls: None
c
c Glossary:
c
c dsurf Output The distances of a point "p" from the bounding
c surfaces at u, u + du, v, v + dv, w and w + dw,
c respectively. Size 6 by np.
c
c fu,fv,fw Output 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, nside, fu, fv, fw.
c Number of points "p".
c
c nside Output Indicates the point is in the volume element, if 1.
c Size np.
c
c pu,pv,pw Input The coordinates of a spatial point "p". 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
c Changes: 1995 February 23 11:10. Added argument dsurf, to store
c distances of the point from the bounding surfaces.
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---- Distances from 6 bounding surfaces.
dimension dsurf (6,1)
c---- Location flag, 1 for inside.
dimension nside (1)
c.... Local variables.
c---- Angular equivalent of 360 degrees.
common /laptbrkn/ circle
c---- Cosine of angle phi at point.
common /laptbrkn/ cosph
c---- Cosine of minimum angle phi.
common /laptbrkn/ cosph1
c---- Cosine of maximum angle phi.
common /laptbrkn/ cosph2
c---- Distance of point from max u.
common /laptbrkn/ dumax
c---- Distance of point from min u.
common /laptbrkn/ dumin
c---- Distance of point from max v.
common /laptbrkn/ dvmax
c---- Distance of point from min v.
common /laptbrkn/ dvmin
c---- Distance of point from max w.
common /laptbrkn/ dwmax
c---- Distance of point from min w.
common /laptbrkn/ dwmin
c---- Estimated error in location of pu.
common /laptbrkn/ erru
c---- Estimated error in location of pv.
common /laptbrkn/ errv
c---- Estimated error in location of pw.
common /laptbrkn/ errw
c---- Angle conversion factor.
common /laptbrkn/ fact
c---- Index in output array.
common /laptbrkn/ n
c---- Square of minimum radius.
common /laptbrkn/ r1sq
c---- Square of maximum radius.
common /laptbrkn/ r2sq
c---- Square of minimum radius.
common /laptbrkn/ r1cb
c---- Square of maximum radius.
common /laptbrkn/ r2cb
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptbrkn finding if points are 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 9902 format (i3,' pu,pv,pw=',1p3e22.14 )
cbug
cbug write ( 3, 9901) np, ksys, iunit, tol, u, du, v, dv, w, dw
cbug do 110 n = 1, np
cbug write ( 3, 9902) n, pu(n), pv(n), pw(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
circle = 360.0
else
fact = 1.0
circle = 2.0 * pi
endif
c.... See if the points are in the brick, and their fractional distances.
do 180 n = 1, np ! Loop over points.
nside(n) = 0
c.... Find the distances of the point from the ends of the brick.
if (du .ge. 0.0) then
dumin = pu(n) - u
dumax = u + du - pu(n)
else
dumin = pu(n) - u - du
dumax = u - pu(n)
endif
if (dv .ge. 0.0) then
dvmin = pv(n) - v
dvmax = v + dv - pv(n)
else
dvmin = pv(n) - v - dv
dvmax = v - pv(n)
endif
c.... Rotate the point 360 degrees, if necessary.
if (ksys .ne. 0) then
150 if (dvmin .lt. 0.0) then
dvmin = dvmin + circle
dvmax = dvmax - circle
go to 150
endif
160 if (dvmax .lt. 0.0) then
dvmin = dvmin - circle
dvmax = dvmax + circle
go to 160
endif
endif
if (dw .ge. 0.0) then
dwmin = pw(n) - w
dwmax = w + dw - pw(n)
else
dwmin = pw(n) - w - dw
dwmax = w - pw(n)
endif
c.... See if the point is in the brick.
erru = tol * (abs (pu(n)) + abs (u) + abs (du))
errv = tol * (abs (pv(n)) + abs (v) + abs (dv))
errw = tol * (abs (pw(n)) + abs (w) + abs (dw))
if (abs (dumin) .lt. erru) dumin = 0.0
if (abs (dumax) .lt. erru) dumax = 0.0
if (abs (dvmin) .lt. errv) dvmin = 0.0
if (abs (dvmax) .lt. errv) dvmax = 0.0
if (abs (dwmin) .lt. errw) dwmin = 0.0
if (abs (dwmax) .lt. errw) dwmax = 0.0
if ((dumin .ge. 0.0) .and.
& (dumax .ge. 0.0) .and.
& (dvmin .ge. 0.0) .and.
& (dvmax .ge. 0.0) .and.
& (dwmin .ge. 0.0) .and.
& (dwmax .ge. 0.0)) then
nside(n) = 1
endif
c.... Find the fractional distances across the brick, on a volume basis.
fu(n) = 0.0
fv(n) = 0.0
fw(n) = 0.0
if (ksys .eq. 0) then ! Rectangular.
if (du .ne. 0.0) fu(n) = dumin / du
if (dv .ne. 0.0) fv(n) = dvmin / dv
if (dw .ne. 0.0) fw(n) = dwmin / dw
elseif (ksys .eq. 1) then ! Cylindrical.
if (du .ne. 0.0) then
r1sq = u**2
r2sq = (u + du)**2
fu(n) = (pu(n)**2 - r1sq) / (r2sq - r1sq)
endif
if (dv .ne. 0.0) fv(n) = dvmin / dv
if (dw .ne. 0.0) fw(n) = dwmin / dw
elseif (ksys .eq. 2) then ! Spherical
if (du .ne. 0.0) then
r1cb = u**3
r2cb = (u + du)**3
fu(n) = (pu(n)**3 - r1cb) / (r2cb - r1cb)
endif
if (dv .ne. 0.0) fv(n) = dvmin / dv
if (dw .ne. 0.0) then
cosph = cos (fact * pw(n))
cosph1 = cos (fact * w)
cosph2 = cos (fact * (w + dw))
fw(n) = (cosph - cosph1) / (cosph2 - cosph1)
endif
endif
c.... Find the distances of the point from the six bounding surfaces.
do 170 ns = 1, 6
dsurf(ns,n) = -1.e99
170 continue
if (ksys .eq. 0) then ! Rectangular.
dsurf(1,n) = u - pu(n)
dsurf(2,n) = u - pu(n) + du
dsurf(3,n) = v - pv(n)
dsurf(4,n) = v - pv(n) + dv
dsurf(5,n) = w - pw(n)
dsurf(6,n) = w - pw(n) + dw
err = tol * (abs (u) + abs (pu(n)))
if (abs (dsurf(1,n)) .le. err) dsurf(1,n) = 0.0
err = tol * (abs (u) + abs (pu(n)) + abs (du))
if (abs (dsurf(2,n)) .le. err) dsurf(2,n) = 0.0
err = tol * (abs (v) + abs (pv(n)))
if (abs (dsurf(3,n)) .le. err) dsurf(3,n) = 0.0
err = tol * (abs (v) + abs (pv(n)) + abs (dv))
if (abs (dsurf(4,n)) .le. err) dsurf(4,n) = 0.0
err = tol * (abs (w) + abs (pw(n)))
if (abs (dsurf(5,n)) .le. err) dsurf(5,n) = 0.0
err = tol * (abs (w) + abs (pw(n)) + abs (dw))
if (abs (dsurf(6,n)) .le. err) dsurf(6,n) = 0.0
elseif (ksys .eq. 1) then ! Cylindrical.
costhmin = cos (fact * v)
sinthmin = sin (fact * v)
costhmax = cos (fact * (v + dv))
sinthmax = sin (fact * (v + dv))
costhp = cos (fact * pv(n))
sinthp = sin (fact * pv(n))
dsurf(1,n) = u - pu(n)
dsurf(2,n) = u - pu(n) + du
dsurf(3,n) = pu(n) * (sinthmin * sinthp - costhmin * costhp)
dsurf(4,n) = pu(n) * (sinthmax * sinthp - costhmax * costhp)
dsurf(5,n) = w - pw(n)
dsurf(6,n) = w - pw(n) + dw
err = tol * (abs (u) + abs (pu(n)))
if (abs (dsurf(1,n)) .le. err) dsurf(1,n) = 0.0
err = tol * (abs (u) + abs (pu(n)) + abs (du))
if (abs (dsurf(2,n)) .le. err) dsurf(2,n) = 0.0
err = tol * abs (pu(n)) *
& (abs (sinthmin * sinthp) + abs (costhmin * costhp))
if (abs (dsurf(3,n)) .le. err) dsurf(3,n) = 0.0
err = tol * abs (pu(n)) *
& (abs (sinthmax * sinthp) + abs (costhmax * costhp))
if (abs (dsurf(4,n)) .le. err) dsurf(4,n) = 0.0
err = tol * (abs (w) + abs (pw(n)))
if (abs (dsurf(5,n)) .le. err) dsurf(5,n) = 0.0
err = tol * (abs (w) + abs (pw(n)) + abs (dw))
if (abs (dsurf(6,n)) .le. err) dsurf(6,n) = 0.0
elseif (ksys .eq. 2) then ! Spherical.
costhmin = cos (fact * v)
sinthmin = sin (fact * v)
costhmax = cos (fact * (v + dv))
sinthmax = sin (fact * (v + dv))
costhp = cos (fact * pv(n))
sinthp = sin (fact * pv(n))
cosphmin = cos (fact * w)
sinphmin = sin (fact * w)
cosphmax = cos (fact * (w + dw))
sinphmax = sin (fact * (w + dw))
cosphp = cos (fact * pw(n))
sinphp = sin (fact * pw(n))
dsurf(1,n) = u - pu(n)
dsurf(2,n) = u - pu(n) + du
dsurf(3,n) = pu(n) * sinphp *
& (sinthmin * costhp - costhmin * sinthp)
dsurf(4,n) = pu(n) * sinphp *
& (sinthmax * costhp - costhmax * sinthp)
dsurf(5,n) = pu(n) * (sinphmin * cosphp - cosphmin * sinphp)
dsurf(6,n) = pu(n) * (sinphmax * cosphp - cosphmax * sinphp)
err = tol * (abs (u) + abs (pu(n)))
if (abs (dsurf(1,n)) .le. err) dsurf(1,n) = 0.0
err = tol * (abs (u) + abs (pu(n)) + abs (du))
if (abs (dsurf(2,n)) .le. err) dsurf(2,n) = 0.0
err = tol * abs (pu(n) * sinphp) *
& (abs (sinthmin * costhp) + abs (costhmin * sinthp))
if (abs (dsurf(3,n)) .le. err) dsurf(3,n) = 0.0
err = tol * abs (pu(n) * sinphp) *
& (abs (sinthmax * costhp) + abs (costhmax * sinthp))
if (abs (dsurf(4,n)) .le. err) dsurf(4,n) = 0.0
err = tol * abs (pu(n)) *
& (abs (sinphmin * cosphp) + abs (cosphmin * sinphp))
if (abs (dsurf(5,n)) .le. err) dsurf(5,n) = 0.0
err = tol * abs (pu(n)) *
& (abs (sinphmax * cosphp) + abs (cosphmax * sinphp))
if (abs (dsurf(6,n)) .le. err) dsurf(6,n) = 0.0
endif ! Tested coordinate system.
180 continue ! End of loop over points.
cbugc***DEBUG begins.
cbug 9903 format (/ 'aptbrkn results:')
cbug 9904 format (i3,' fu,fv,fw=',1p3e22.14 )
cbug 9905 format (i3,' nside=',i2)
cbug 9906 format (i3,' distances of point from surfaces:' /
cbug & ' d umin,umax=',1p2e22.14 /
cbug & ' d vmin,vmax=',1p2e22.14 /
cbug & ' d wmin,wmax=',1p2e22.14 )
cbug
cbug write ( 3, 9903)
cbug do 205 n = 1, np
cbug write ( 3, 9904) n, fu(n), fv(n), fw(n)
cbug write ( 3, 9905) n, nside(n)
cbug write ( 3, 9906) n, (dsurf(nn,n), nn = 1, 6)
cbug 205 continue
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptbrkn. (+1 line.)
end
UCRL-WEB-209832