subroutine aptelis (cenw, aw, ar, bw, br, np, tol, qc, qu, qv, qw,
& quv, qvw, qwu, quu, qvv, qww, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTELIS
c
c call aptelis (cenw, aw, ar, bw, br, np, tol, qc, qu, qv, qw,
c & quv, qvw, qwu, quu, qvv, qww, nerr)
c
c Version: aptelis Updated 1990 December 12 14:30.
c aptelis Originated 1990 December 10 17:30.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: For each of np sets of input data, to find the equation of the
c quadric surface symmetric around the w axis, with a symmetry
c center at cenw on the w axis, and passing through the two
c distinct points a = (aw, ar) and b = (bw, br), where
c ar = sqrt (bu**2 + bv**2) and br = sqrt (cu**2 + cv**2)
c are radial distances from the w axis.
c The coordinates (u,v,w) are along the three major orthogonal
c axes, and may represent (x,y,z), (y,z,x) or (z,x,y).
c
c The equation for the quadric surface is as follows:
c
c f(u,v,w) = qc + qu * u + qv * v + qw * w +
c quv * u * v + qvw * v * w + qwu * w * u +
c quu * u**2 + qvv * v**2 + qww * w**2 = 0
c
c where qu = qv = quv = qvw = qwu = 0, and quu = qvv.
c
c The surface crosses the w axis (r = 0) at the coordinates
c w = cenw +/- sqrt (abm / (ar**2 - br**2)),
c where abm = ar**2 * (bw - cenw)**2 - br**2 * (aw - cenw)**2.
c
c General cases:
c
c If the radial distances ar and br decrease with axial distance
c from cenw, the surface will be an ellipsoid which is circular in
c a w plane.
c
c If the radial distances ar and br increase with axial distance
c from cenw, the surface will be a hyperboloid of two sheets, also
c circular in a w plane.
c
c Indeterminate cases:
c
c If ar**2 = br**2, and (aw - cenw)**2 = (bw - cens)**2, within
c the limits determined by tol, the surface is indeterminate, and
c all coefficients will be zero.
c
c Other singular cases:
c
c If (aw - cenw)**2 = (bw - cens)**2, within the limits determined
c by tol, the surface will be two parallel w planes.
c
c If ar**2 = br**2, within the limits determined by tol, the
c surface will be a right circular cylinder on the w axis.
c
c If ar**2 * (bw - cenw)**2 = br**2 * (aw - cenw)**2, within the
c limits determined by tol, the surface will be a right circular
c cone on the w axis with the vertex at cenw.
c
c
c Note: Use apttris and aptrois to translate and rotate the surface
c into any desired position and orientation.
c Use aptesis to find the equation of an ellipsoid given the
c three semiaxes on the major axes.
c Use aptcnis or aptcois to find the equation of a cone
c symmetric around a major axis.
c Use aptcyis to find the equation of a cylinder symmetric
c around a major axis.
c
c Input: cenw, aw, ar, bw, br, np, tol.
c
c Output: qc, qu, qv, qw, quv, qvw, qwu, quu, qvv, qww, nerr.
c
c Calls: aptvdic, aptvaxc
c
c
c Glossary:
c
c ar, aw Input Coordinates r and w of point "a" on the surface.
c Size np. Must be distinct from point "b".
c
c br, bw Input Coordinates r and w of point "b" on the surface.
c Size np. Must be distinct from point "a".
c
c cenw Input Coordinate w of the symmetry center of the surface.
c Size np. On the w axis.
c
c nerr Output Indicates an input error, if not zero:
c 1 if np is not positive.
c
c np Input Size of arrays cenw, aw, ar, bw, br, qc, qu, qv, qw,
c quv, qvw, qwu, quu, qvv, qww.
c
c qc Output Constant term in the equation of the surface. Size np.
c Zero only if the surface is centered on the origin.
c qc = (aw - cenw)**2 * br**2 - (bw - bcen)**2 * ar**2
c + (ar**2 - br**2) * cenw**2.
c
c qu,qv,qw Output Coefficients of u, v, w in the equation of the surface.
c Size np.
c qu = qv = 0,
c qw = 2.0 * (br**2 - ar**2) * cenw.
c
c q.. Output Coefficients of second-order terms in the equation of
c the surface: quv, qvw, qwu, quu, qvv, qww.
c Coefficients of u*v, v*w, w*u, u**2, v**2, w**2,
c respectively. Size np.
c quv = qvw = qwu = 0,
c quu = qvv = (cz - cenz)**2 - (bz - cenz)**2,
c qww = ar**2 - br**2.
c
c tol Input Numerical tolerance limit.
c On Cray computers, recommend 1.e-5 to 1.e-11.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c---- Coordinate r of point "a".
dimension ar (1)
c---- Coordinate w of point "a".
dimension aw (1)
c---- Coordinate r of point "b".
dimension br (1)
c---- Coordinate w of point "b".
dimension bw (1)
c---- Coordinate w of center of surface.
dimension cenw (1)
c---- Coefficient of term 1.0.
dimension qc (1)
c---- Coefficient of term u.
dimension qu (1)
c---- Coefficient of term u**2.
dimension quu (1)
c---- Coefficient of term u * v.
dimension quv (1)
c---- Coefficient of term v.
dimension qv (1)
c---- Coefficient of term v**2.
dimension qvv (1)
c---- Coefficient of term v * w.
dimension qvw (1)
c---- Coefficient of term w.
dimension qw (1)
c---- Coefficient of term w * u.
dimension qwu (1)
c---- Coefficient of term w**2.
dimension qww (1)
c.... Local variables.
c---- Component r of "aa", ar**2.
common /laptelis/ aar (64)
c---- Component w of "aa", (aw - cenw)**2.
common /laptelis/ aaw (64)
c---- Estimated error in aaw, based on tol.
common /laptelis/ aawerr
c---- Magnitude of "bb" - "aa".
common /laptelis/ abd (64)
c---- Cross product of "aa" and "bb".
common /laptelis/ abm (64)
c---- Difference br**2 - ar**2.
common /laptelis/ abr (64)
c---- Difference (bw-cenw)**2-(aw-cenw)**2.
common /laptelis/ abw (64)
c---- Component r of "bb", br**2.
common /laptelis/ bbr (64)
c---- Component w of "bb", (bw - cenw)**2.
common /laptelis/ bbw (64)
c---- Estimated error in bbw, based on tol.
common /laptelis/ bbwerr
c---- Index in local array.
common /laptelis/ n
c---- First index of subset of data.
common /laptelis/ n1
c---- Last index of subset of data.
common /laptelis/ n2
c---- Index in external array.
common /laptelis/ nn
c---- Size of current subset of data.
common /laptelis/ ns
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptelis finding the equation of an quadric surface.',
cbug & ' tol=',1pe13.5 /
cbug & ' Surface has center, two points:')
cbug 9902 format (
cbug & i3,' cenw=',1pe22.14 /
cbug & ' aw,ar=',1p2e22.14 /
cbug & ' bw,br=',1p2e22.14)
cbug write ( 3, 9901) tol
cbug write ( 3, 9902) (n, cenw(n), aw(n), ar(n), bw(n), br(n),
cbug & n = 1, np)
cbugc***DEBUG ends.
c.... Initialize.
nerr = 0
c.... Test for input errors.
if (np .le. 0) then
nerr = 1
cbugc***DEBUG begins.
cbug 9903 format (/ "aptelis fatal. Non-positive np =",i5)
cbug write ( 3, 9903) np
cbugc***DEBUG ends.
go to 210
endif
c.... Set up the indices of the first subset of data.
n1 = 1
n2 = min (np, 64)
110 ns = n2 - n1 + 1
c.... Find vectors "aa" and "bb".
c---- Loop over subset of data.
do 120 n = 1, ns
nn = n + n1 - 1
aaw(n) = (aw(nn) - cenw(nn))**2
aar(n) = ar(nn)**2
bbw(n) = (bw(nn) - cenw(nn))**2
bbr(n) = br(nn)**2
c---- End of loop over subset of data.
120 continue
c.... See if truncation tests are to be done.
c---- Test for truncation to zero.
if (tol .gt. 0.0) then
c---- Loop over subset of data.
do 130 n = 1, ns
nn = n + n1 - 1
aawerr = 2.0 * tol * (abs (aw(nn)) + abs (cenw(nn)))**2
bbwerr = 2.0 * tol * (abs (bw(nn)) + abs (cenw(nn)))**2
if (aaw(n) .lt. aawerr) then
aaw(n) = 0.0
endif
if (bbw(n) .lt. bbwerr) then
bbw(n) = 0.0
endif
c---- End of loop over subset of data.
130 continue
c---- Tested tol.
endif
cc.... Find the 2-D vector from "aa" to "bb".
call aptvdic (aaw, aar, bbw, bbr, ns, tol, abw, abr, abd, nerr)
c.... Find the cross product of the 2-D vectors "aa" and "bb".
call aptvaxc (aaw, aar, bbw, bbr, ns, tol, abm, nerr)
c.... Find the coefficients of the equation of the surface.
c---- Loop over subset of data.
do 140 n = 1, ns
nn = n + n1 - 1
qc(nn) = abm(n) - abr(n) * cenw(nn)**2
qu(nn) = 0.0
qv(nn) = 0.0
qw(nn) = 2.0 * abr(n) * cenw(nn)
quv(nn) = 0.0
qvw(nn) = 0.0
qwu(nn) = 0.0
quu(nn) = abw(n)
qvv(nn) = quu(nn)
qww(nn) = -abr(n)
c---- End of loop over subset of data.
140 continue
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 9904 format (/ 'aptelis results:' /
cbug & (i3,' qc= ',1pe22.14 /
cbug & ' qu,qv,qw= ',1p3e22.14 /
cbug & ' quv,qvw,qwu=',1p3e22.14 /
cbug & ' quu,qvv,qww=',1p3e22.14))
cbug write ( 3, 9904) (n, qc(n), qu(n), qv(n), qw(n),
cbug & quv(n), qvw(n), qwu(n), quu(n), qvv(n), qww(n),
cbug & n = 1, np)
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptelis. (+1 line.)
end
UCRL-WEB-209832