subroutine aptqupt (qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz,
& tol, aqp, qpx, qpy, qpz, nqp, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTQUPT
c
c call aptqupt (qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz,
c tol, aqp, qpx, qpy, qpz, nqp, nerr)
c
c Version: aptqupt Updated 1994 November 9 12:20.
c aptqupt Originated 1994 November 9 12:20.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: For the quadric surface specified by the coordinates qc to qzz,
c to find up to 6 points of the surface which are extreme points
c in the x, y and z directions, up to 6 points which are
c intercepts on the x, y and z axes, and up to 6 points which are
c intercepts on the symmetry axes (x', y' and z') of the quadric
c surface.
c
c The implicit equation of the quadric curve is:
c
c F(x,y,z) = qc + qx * x + qy * y +
c qxy * x * y + qyz * y * z + qzx * z * x +
c qxx * x**2 + qyy * y**2 + qzz * z**2 = 0 .
c
c Results may be unpredictable if the coefficients of the quadric
c surface differ by many orders of magnitude.
c
c Method: See subroutine aptqext for a description of the method of
c finding extrema.
c
c Any intercepts on the x, y and z axes or the symmetry axes are
c solutions of the equations (in original or standard form):
c
c qc + qx * x + qxx * x**2 = 0
c qc + qy * y + qyy * y**2 = 0
c qc + qz * z + qzz * z**2 = 0
c
c Input: qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz, tol.
c
c Output: nqp, aqp, qpx, qpy, qpz, nerr.
c
c Calls: aptaxis, aptmopv, aptqext, aptqrts, apttran
c
c
c Glossary:
c
c aqp Output Type of quadric surface point. Array size must be at
c least 18.
c Minimum: "Min x", "Min y" or "Min z".
c Maximum: "Max x", "Max y" or "Max z".
c Aligned planes: "Const x", "Const y" or "Const z".
c Intersection of two planes, perpendicular to a major
c axis:
c "Inter x", "Inter y" or "Inter z".
c Saddle point: "Saddle x", Saddle y" or "Saddle z".
c Singular isolated point: "Point x", "Point y" or
c "Point z".
c X axis intercept: "Int x 1" or "Int x 2".
c Y axis intercept: "Int y 1" or "Int y 2".
c Z axis intercept: "Int z 1" or "Int z 2".
c X' axis intercept: "Int x' 1" or "Int x' 2".
c Y' axis intercept: "Int y' 1" or "Int y' 2".
c Z' axis intercept: "Int z' 1" or "Int z' 2".
c If none, return "unknown".
c
c qpx,y,z Output The x, y and z coordinates of an extreme point, a
c saddle point, or an intersection of two planes, for
c which the x, y or z coordinate is a maximum or a
c minimum, or an intercept on the x, y or z axis, or
c the x', y' or z' axis. Array size must be at least
c 18. If none, return (-1.e99,-1.e99,-1.e99).
c
c nerr Output Indicates an input error, if not zero:
c 1 if all coefficients, except possible qc, are zero,
c or the quadric surface type is imaginary.
c
c nqp Output Number of extreme and intercept points found.
c Never more than 18.
c
c qc In Constant term in the equation of the quadric curve.
c
c qx,qy,qz In Coefficients of x, y and z, respectively, in the
c equation of the quadric curve.
c
c qxx In Coefficient of x**2 in the equation of the quadric
c curve.
c
c qxy In Coefficient of x * y in the equation of the quadric
c curve.
c
c qyy In Coefficient of y**2 in the equation of the quadric
c curve.
c
c qyz In Coefficient of y * z in the equation of the quadric
c curve.
c
c qzx In Coefficient of z * x in the equation of the quadric
c curve.
c
c qzz In Coefficient of z**2 in the equation of the quadric
c curve.
c
c tol Input Numerical tolerance limit. With 64-bit floating
c point numbers, 1.e-5 to 1.e-11 is recommended.
c On 32-bit computers, recommend 1.e-6.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
dimension aqp(1) ! Type of quadric surface point.
character*8 aqp ! Type of quadric surface point.
dimension qpx(1) ! Coordinate x of a point in quadric nq.
dimension qpy(1) ! Coordinate y of a point in quadric nq.
dimension qpz(1) ! Coordinate z of a point in quadric nq.
c.... Local variables.
common /lqupt/ ex(6) ! Coordinate x of extreme point.
common /lqupt/ ey(6) ! Coordinate y of extreme point.
common /lqupt/ ez(6) ! Coordinate z of extreme point.
common /lqupt/ adir(6) ! Type of extreme point.
character*8 adir ! Type of extreme point.
common /lqupt/ rotq(3,3) ! Standard form rotation operator.
c=======================================================================********
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptqupt finding points on a quadric surface.' /
cbug & ' tol=',1pe13.5 )
cbug 9902 format (
cbug & ' qc= ',1pe22.14 /
cbug & ' qx,qy,qz= ',1p3e22.14 /
cbug & ' qxy,qyz,qzx=',1p3e22.14 /
cbug & ' qxx,qyy,qzz=',1p3e22.14)
cbug write ( 3, 9901) tol
cbug write ( 3, 9902) qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz
cbugc***DEBUG ends.
c.... Initialize.
nerr = 0
nqp = 0
do 110 n = 1, 18
aqp(n) = 'unknown'
qpx(n) = -1.e99
qpy(n) = -1.e99
qpz(n) = -1.e99
110 continue
c=======================================================================********
c.... Do the principal axis transformation on the quadric surface.
ac = qc
ax = qx
ay = qy
az = qz
axy = qxy
ayz = qyz
azx = qzx
axx = qxx
ayy = qyy
azz = qzz
tolx = tol
if (tolx .le. 0.0) tolx = 1.e-11
call aptaxis (ac, ax, ay, az, axy, ayz, azx, axx, ayy, azz,
& tol, cenx, ceny, cenz, rotq, ntype,
& nerra)
if (ntype .gt. 16) then
nerr = 1
go to 210
endif
if (nerra .ne. 0) then
nerr = 1
go to 210
endif
c.... See if the quadric surface is already in standard form.
kstd = 0
if ((rotq(1,1) .eq. 1.0) .and.
& (rotq(2,2) .eq. 1.0) .and.
& (rotq(3,3) .eq. 1.0) .and.
& (cenx .eq. 0.0) .and.
& (ceny .eq. 0.0) .and.
& (cenz .eq. 0.0)) then
kstd = 1
endif
c=======================================================================********
c.... Find any extreme points on the quadric surface.
call aptqext (qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz,
& tol, ntyps, cenx, ceny, cenz, np, ex, ey, ez,
& adir, nerra)
if (np .gt. 0) then
do 120 nex = 1, np
nqp = nqp + 1
aqp(nqp) = adir(nex)
qpx(nqp) = ex(nex)
qpy(nqp) = ey(nex)
qpz(nqp) = ez(nex)
120 continue
endif
c=======================================================================********
c.... Find any intercepts of the original form on the x axis.
call aptqrts (0, qxx, qx, qc, qq, tol,
& nroots, xint1, xint2, itrun)
if (nroots .ge. 1) then
nqp = nqp + 1
aqp(nqp) = 'Int x 1'
qpx(nqp) = xint1
qpy(nqp) = 0.0
qpz(nqp) = 0.0
endif
if (nroots .eq. 2) then
nqp = nqp + 1
aqp(nqp) = 'Int x 2'
qpx(nqp) = xint2
qpy(nqp) = 0.0
qpz(nqp) = 0.0
endif
c.... Find any intercepts of the original form on the y axis.
call aptqrts (0, qyy, qy, qc, qq, tol,
& nroots, yint1, yint2, itrun)
if (nroots .ge. 1) then
nqp = nqp + 1
aqp(nqp) = 'Int y 1'
qpx(nqp) = 0.0
qpy(nqp) = yint1
qpz(nqp) = 0.0
endif
if (nroots .eq. 2) then
nqp = nqp + 1
aqp(nqp) = 'Int y 2'
qpx(nqp) = 0.0
qpy(nqp) = yint2
qpz(nqp) = 0.0
endif
c.... Find any intercepts of the original form on the z axis.
call aptqrts (0, qzz, qz, qc, qq, tol,
& nroots, zint1, zint2, itrun)
if (nroots .ge. 1) then
nqp = nqp + 1
aqp(nqp) = 'Int z 1'
qpx(nqp) = 0.0
qpy(nqp) = 0.0
qpz(nqp) = zint1
endif
if (nroots .eq. 2) then
nqp = nqp + 1
aqp(nqp) = 'Int z 2'
qpx(nqp) = 0.0
qpy(nqp) = 0.0
qpz(nqp) = zint2
endif
c=======================================================================********
c.... Skip stuff for standard form if quadric is already standard.
if (kstd .eq. 1) go to 210
c=======================================================================********
c.... Find intercepts of standard form on x' axis, original coordinates.
dx = -cenx
dy = -ceny
dz = -cenz
call aptqrts (0, axx, 0., ac, qq, tol,
& nroots, xint1, xint2, itrun)
if (nroots .ge. 1) then
yr = 0.0
zr = 0.0
call aptmopv (rotq, 1, 0., 0., 0., xint1, yr, zr, 1, tol,
& nerra)
call apttran (dx, dy, dz, xint1, yr, zr, 1, tol, nerra)
nqp = nqp + 1
aqp(nqp) = 'Int x'' 1'
qpx(nqp) = xint1
qpy(nqp) = yr
qpz(nqp) = zr
endif
if (nroots .eq. 2) then
yr = 0.0
zr = 0.0
call aptmopv (rotq, 1, 0., 0., 0., xint2, yr, zr, 1, tol,
& nerr)
call apttran (dx, dy, dz, xint2, yr, zr, 1, tol, nerra)
nqp = nqp + 1
aqp(nqp) = 'Int x'' 2'
qpx(nqp) = xint2
qpy(nqp) = yr
qpz(nqp) = zr
endif
c=======================================================================********
c.... Find intercepts of standard form on y' axis, original coordinates.
call aptqrts (0, ayy, 0., ac, qq, tol,
& nroots, yint1, yint2, itrun)
if (nroots .ge. 1) then
xr = 0.0
zr = 0.0
call aptmopv (rotq, 1, 0., 0., 0., xr, yint1, zr, 1, tol,
& nerra)
call apttran (dx, dy, dz, xr, yint1, zr, 1, tol, nerra)
nqp = nqp + 1
aqp(nqp) = 'Int y'' 1'
qpx(nqp) = xr
qpy(nqp) = yint1
qpz(nqp) = zr
endif
if (nroots .eq. 2) then
xr = 0.0
zr = 0.0
call aptmopv (rotq, 1, 0., 0., 0., xr, yint2, zr, 1, tol,
& nerra)
call apttran (dx, dy, dz, xr, yint2, zr, 1, tol, nerra)
nqp = nqp + 1
aqp(nqp) = 'Int y'' 2'
qpx(nqp) = xr
qpy(nqp) = yint2
qpz(nqp) = zr
endif
c=======================================================================********
c.... Find intercepts of standard form on z' axis, original coordinates.
call aptqrts (0, azz, az, ac, qq, tol,
& nroots, zint1, zint2, itrun)
if (nroots .ge. 1) then
xr = 0.0
yr = 0.0
call aptmopv (rotq, 1, 0., 0., 0., xr, yr, zint1, 1, tol,
& nerra)
call apttran (dx, dy, dz, xr, yr, zint1, 1, tol, nerra)
nqp = nqp + 1
aqp(nqp) = 'Int z'' 1'
qpx(nqp) = xr
qpy(nqp) = yr
qpz(nqp) = zint1
endif
if (nroots .eq. 2) then
xr = 0.0
yr = 0.0
call aptmopv (rotq, 1, 0., 0., 0., xr, yr, zint2, 1, tol,
& nerra)
call apttran (dx, dy, dz, xr, yr, zint2, 1, tol, nerra)
nqp = nqp + 1
aqp(nqp) = 'Int z'' 2'
qpx(nqp) = xr
qpy(nqp) = yr
qpz(nqp) = zint2
endif
c=======================================================================********
210 continue
cbugc***DEBUG begins.
cbug 9910 format (/ 'aptqupt results: nqp=',i2,' nerr=',i2)
cbug 9912 format (a8,' xyz=',1p3e22.14)
cbug write ( 3, 9910) nqp, nerr
cbug do 991 n = 1, nqp
cbug write ( 3, 9912) aqp(n), qpx(n), qpy(n), qpz(n)
cbug 991 continue
cbugc***DEBUG ends.
return
c.... End of subroutine aptqupt. (+1 line)
end
UCRL-WEB-209832