subroutine aptbrkv (ksys, iunit, u, du, v, dv, w, dw, tol,
& umin, vmin, wmin, ucen, vcen, wcen,
& umax, vmax, wmax,
& du11, dv11, dw11, du12, dv12, dw12,
& ducn, dvcn, dwcn, du21, dv21, dw21,
& du22, dv22, dw22,
& aumin, avmin, awmin, aucen, avcen, awcen,
& aumax, avmax, awmax, vol, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTBRKV
c
c call aptbrkv (ksys, iunit, u, du, v, dv, w, dw, tol,
c & umin, vmin, wmin, ucen, vcen, wcen,
c & umax, vmax, wmax,
c & du11, dv11, dw11, du12, dv12, dw12,
c & ducn, dvcn, dwcn, du21, dv21, dw21,
c & du22, dv22, dw22,
c & aumin, avmin, awmin, aucen, avcen, awcen,
c & aumax, avmax, awmax, vol, nerr)
c
c Version: aptbrkv Updated 1994 June 6 10:50.
c aptbrkv Originated 1994 March 3 12:00
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: For the 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 to find the minimum, maximum and average coordinates, the
c 12 edge lengths, and the 3 distances through the centroid,
c the 6 face areas, and the 3 areas through the centroid,
c and the volume. Flag nerr indicates any input error.
c
c Input: ksys, iunit, u, du, v, dv, w, dw
c
c Output:
c
c Calls: None
c
c Glossary:
c
c aumin,... Output The face areas at umin, ucen and umax.
c
c avmin,... Output The face areas at vmin, vcen and vmax.
c
c awmin,... Output The face areas at wmin, wcen and wmax.
c
c du11,... Output The edge lengths of the brick in the u direction.
c 11 indicates at minimum v, minimum w.
c 12 indicates at minimum v, maximum w.
c cn indicates at centroid v, centroid w.
c 21 indicates at maximum v, minimum w.
c 22 indicates at maximum v, maximum w.
c
c dv11,... Output The edge lengths of the brick in the v direction.
c 11 indicates at minimum u, minimum w, etc.
c
c dw11,... Output The edge lengths of the brick in the w direction.
c 11 indicates at minimum u, minimum v, etc.
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 ksys is 1 or 2 and u or u + du is negative.
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 u,v,wcen Output The centroid values of the coordinates u, v, w.
c
c u,v,wmax Output The maximum values of the coordinates u, v, w.
c
c u,v,wmin Output The minimum values of the coordinates u, v, w.
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 vol Output The volume of the brick.
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.... Local variables.
common /laptbrkv/ cosw
common /laptbrkv/ cosw1
common /laptbrkv/ cosw2
common /laptbrkv/ dcos
common /laptbrkv/ errd
common /laptbrkv/ pdu
common /laptbrkv/ pdv
common /laptbrkv/ pdw
common /laptbrkv/ r1cb
common /laptbrkv/ r2cb
common /laptbrkv/ r1sq
common /laptbrkv/ r2sq
common /laptbrkv/ swcen
common /laptbrkv/ swmax
common /laptbrkv/ swmin
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptbrkv finding data for a brick. ksys=',i3,
cbug & ' 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) ksys, iunit, tol, u, du, v, dv, w, dw
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
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.... Find the angle conversion factor.
if (iunit .eq. 0) then
fact = pi / 180.0
else
fact = 1.0
endif
c.... Find the limiting coordinate system values.
if (du .ge. 0.0) then
umin = u
umax = u + du
if (abs (umax) .le. tol * (abs (u) + du)) umax = 0.0
else
umin = u + du
umax = u
if (abs (umin) .le. tol * (abs (u) - du)) umin = 0.0
endif
if (dv .ge. 0.0) then
vmin = v
vmax = v + dv
if (abs (vmax) .le. tol * (abs (v) + dv)) vmax = 0.0
else
vmin = v + dv
vmax = v
if (abs (vmin) .le. tol * (abs (v) - dv)) vmin = 0.0
endif
if (dw .ge. 0.0) then
wmin = w
wmax = w + dw
if (abs (wmax) .le. tol * (abs (w) + dw)) wmax = 0.0
else
wmin = w + dw
wmax = w
if (abs (wmin) .le. tol * (abs (w) - dw)) wmin = 0.0
endif
c.... Find the absolute coordinate increments.
pdu = abs (du)
pdv = abs (dv)
pdw = abs (dw)
c.... Find the centroid in the brick.
if (ksys .eq. 0) then
ucen = u + 0.5 * du
if (abs (ucen) .le. tol * (abs (u) + 0.5 * abs (du))) ucen = 0.0
vcen = v + 0.5 * dv
if (abs (vcen) .le. tol * (abs (v) + 0.5 * abs (dv))) vcen = 0.0
wcen = w + 0.5 * dw
if (abs (wcen) .le. tol * (abs (w) + 0.5 * abs (dw))) wcen = 0.0
elseif (ksys .eq. 1) then
if (umin .lt. 0.0) then
nerr = 3
go to 210
endif
r1sq = u**2
r2sq = (u + du)**2
if (abs (u + du) .le. tol * (abs (u) + abs (du))) r2sq = 0.0
ucen = sqrt (0.5 * (r1sq + r2sq))
vcen = v + 0.5 * dv
if (abs (vcen) .le. tol * (abs (v) + 0.5 * abs (dv))) vcen = 0.0
wcen = w + 0.5 * dw
if (abs (wcen) .le. tol * (abs (w) + 0.5 * abs (dw))) wcen = 0.0
elseif (ksys .eq. 2) then
if (umin .lt. 0.0) then
nerr = 3
go to 210
endif
r1cb = u**3
r2cb = (u + du)**3
if (abs (u + du) .le. tol * (abs (u) + abs (du))) r2cb = 0.0
cosw1 = cos (fact * w)
arg = w + dw
if (abs (arg) .le. tol * (abs (w) + abs (dw))) arg = 0.0
cosw2 = cos (fact * arg)
cosw = 0.5 * (cosw1 + cosw2)
errc = 0.5 * tol * (abs (cosw1) + abs (cosw2))
if (abs (cosw) .le. errc) cosw = 0.0
ucen = (0.5 * (r1cb + r2cb))**(1.0 / 3.0)
vcen = v + 0.5 * dv
if (abs (vcen) .le. tol * (abs (v) + 0.5 * abs (dv))) vcen = 0.0
wcen = acos (cosw) / fact
endif
c.... Find the edge lengths.
if (ksys .eq. 0) then
du11 = pdu ! dx.
du12 = pdu ! dx.
du21 = pdu ! dx.
du22 = pdu ! dx.
ducn = pdu ! dx.
dv11 = pdv ! dy.
dv12 = pdv ! dy.
dv21 = pdv ! dy.
dv22 = pdv ! dy.
dvcn = pdv ! dy.
dw11 = pdw ! dz.
dw12 = pdw ! dz.
dw21 = pdw ! dz.
dw22 = pdw ! dz.
dwcn = pdw ! dz.
elseif (ksys .eq. 1) then
du11 = pdu ! d(rcyl).
du12 = pdu ! d(rcyl).
du21 = pdu ! d(rcyl).
du22 = pdu ! d(rcyl).
ducn = pdu ! d(rcyl).
dv11 = umin * fact * pdv ! rcyl * d(theta).
dv12 = umin * fact * pdv ! rcyl * d(theta).
dv21 = umax * fact * pdv ! rcyl * d(theta).
dv22 = umax * fact * pdv ! rcyl * d(theta).
dvcn = ucen * fact * pdv ! rcyl * d(theta).
dw11 = pdw ! dz.
dw12 = pdw ! dz.
dw21 = pdw ! dz.
dw22 = pdw ! dz.
dwcn = pdw ! dz.
elseif (ksys .eq. 2) then
swmin = sin (fact * wmin)
if (abs (swmin) .le. tol) swmin = 0.0
swmax = sin (fact * wmax)
if (abs (swmax) .le. tol) swmax = 0.0
swcen = sin (fact * wcen)
if (abs (swcen) .le. tol) swcen = 0.0
du11 = pdu ! d(rsph).
du12 = pdu ! d(rsph).
du21 = pdu ! d(rsph).
du22 = pdu ! d(rsph).
ducn = pdu ! d(rsph).
dv11 = umin * fact * pdv * swmin ! rsph1 * d(theta) * sin (phi1).
dv12 = umin * fact * pdv * swmax ! rsph1 * d(theta) * sin (phi2).
dv21 = umax * fact * pdv * swmin ! rsph2 * d(theta) * sin (phi1).
dv22 = umax * fact * pdv * swmax ! rsph2 * d(theta) * sin (phi2).
dvcn = ucen * fact * pdv * swcen ! rsphcn * d(theta) * sin (phicn).
dw11 = umin * fact * pdw ! rsph1 * d(phi).
dw12 = umin * fact * pdw ! rsph1 * d(phi).
dw21 = umax * fact * pdw ! rsph2 * d(phi).
dw22 = umax * fact * pdw ! rsph2 * d(phi).
dwcn = ucen * fact * pdw ! rsphcn * d(phi).
endif
c.... Find the face areas, volume.
if (ksys .eq. 0) then
aumin = pdv * pdw
aumax = pdv * pdw
aucen = pdv * pdw
avmin = pdw * pdu
avmax = pdw * pdu
avcen = pdw * pdu
awmin = pdu * pdv
awmax = pdu * pdv
awcen = pdu * pdv
vol = pdu * pdv * pdw
elseif (ksys .eq. 1) then
aumin = umin * fact * pdv * pdw
aumax = umax * fact * pdv * pdw
aucen = ucen * fact * pdv * pdw
avmin = pdw * pdu
avmax = pdw * pdu
avcen = pdw * pdu
arg = umax**2 - umin**2
if (abs (arg) .le. tol * (umax**2 + umin**2)) arg = 0.0
awmin = arg * fact * pdv / 2.0
awmax = arg * fact * pdv / 2.0
awcen = arg * fact * pdv / 2.0
vol = arg * fact * pdv * pdw / 2.0
elseif (ksys .eq. 2) then
cwmin = cos (fact * wmin)
cwmax = cos (fact * wmax)
dcos = cwmin - cwmax
errd = tol * (abs (cwmin) + abs (cwmax))
if (abs (dcos) .le. errd) dcos = 0.0
swmin = sin (fact * wmin)
if (abs (swmin) .le. tol) swmin = 0.0
swmax = sin (fact * wmax)
if (abs (swmax) .le. tol) swmax = 0.0
swcen = sin (fact * wcen)
if (abs (swcen) .le. tol) swcen = 0.0
aumin = umin**2 * fact * pdv * dcos
aumax = umax**2 * fact * pdv * dcos
aucen = ucen**2 * fact * pdv * dcos
arg = umax**2 - umin**2
if (abs (arg) .le. tol * (umax**2 + umin**2)) arg = 0.0
avmin = (umax**2 - umin**2) * fact * pdw / 2.0
avmax = (umax**2 - umin**2) * fact * pdw / 2.0
avcen = (umax**2 - umin**2) * fact * pdw / 2.0
awmin = arg * swmin * fact * pdv / 2.0
awmax = arg * swmax * fact * pdv / 2.0
awcen = arg * swcen * fact * pdv / 2.0
arg = umax**3 - umin**3
if (abs (arg) .le. tol * (umax**3 + umin**3)) arg = 0.0
vol = arg * fact * pdv * dcos / 3.0
endif
210 continue
cbugc***DEBUG begins.
cbug 9903 format (/ 'aptbrkv results: nerr=',i3)
cbug 9904 format (' u,v,w min',1p3e22.14)
cbug 9905 format (' u,v,w cen',1p3e22.14)
cbug 9906 format (' u,v,w max',1p3e22.14)
cbug 9907 format (' du,dv,dw 11',1p3e22.14)
cbug 9908 format (' du,dv,dw 12',1p3e22.14)
cbug 9909 format (' du,dv,dw cn',1p3e22.14)
cbug 9910 format (' du,dv,dw 21',1p3e22.14)
cbug 9911 format (' du,dv,dw 22',1p3e22.14)
cbug 9912 format (' au,av,aw min',1p3e22.14)
cbug 9913 format (' au,av,aw cen',1p3e22.14)
cbug 9914 format (' au,av,aw max',1p3e22.14)
cbug 9915 format (' volume ',1pe22.14)
cbug write ( 3, 9903) nerr
cbug write ( 3, 9904) umin, vmin, wmin
cbug write ( 3, 9905) ucen, vcen, wcen
cbug write ( 3, 9906) umax, vmax, wmax
cbug write ( 3, 9907) du11, dv11, dw11
cbug write ( 3, 9908) du12, dv12, dw12
cbug write ( 3, 9909) ducn, dvcn, dwcn
cbug write ( 3, 9910) du21, dv21, dw21
cbug write ( 3, 9911) du22, dv22, dw22
cbug write ( 3, 9912) aumin, avmin, awmin
cbug write ( 3, 9913) aucen, avcen, awcen
cbug write ( 3, 9914) aumax, avmax, awmax
cbug write ( 3, 9915) vol
cbugc***DEBUG ends.
return
c.... End of subroutine aptbrkv. (+1 line.)
end
UCRL-WEB-209832