subroutine aptqupr (ax, ay, az, qc, qx, qy, qz, qxy, qyz, qzx,
& qxx, qyy, qzz, tol, sc, sx, sy, sz,
& sxy, syz, szx, sxx, syy, szz, ntype, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTQUPR
c
c call aptqupr (ax, ay, az, qc, qx, qy, qz, qxy, qyz, qzx,
c & qxx, qyy, qzz, tol, sc, sx, sy, sz,
c & sxy, syz, szx, sxx, syy, szz, ntype, nerr)
c
c Version: aptqupr Updated 1995 January 31 14:40.
c aptqupr Originated 1995 January 31 14:40.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: Find any projection in the direction a = (ax, ay, az) of the
c quadric surface with implicit equation:
c
c F(x,y,z) = qc + qx * x + qy * y + qz * z +
c qxy * x * y + qyz * y * z + qzx * z * x +
c qxx * x**2 + qyy * y**2 + qzz * z**2 = 0
c
c The result, if any, will be the cylindrical quadric surface with
c its cylindrical axis in the direction of vector "a", tangent to
c the original quadric surface, and the implicit equation with
c coefficients sc, ..., szz, of type ntype.
c
c The original quadric surface has the normal vector:
c
c N(x,y,z) = (NX, NY, NZ), where
c NX = qx + 2 * qxx * x + qxy * y + qzx * z
c NY = qy + qxy * x + 2 * qyy * y + qyz * z
c NZ = qz + qzx * x + qyz * y + 2 * qzz * z
c
c The curve of tangency of the original quadric surface and the
c projected cylindrical quadric surface in in the plane determined
c by setting the dot product of vector "a" and the normal vector
c N to zero, if such a plane exists:
c
c a dot N = ax * NX + ay * NY + az * NZ = 0
c
c Any input error, or nonexistance of a result, will be indicated
c by a nonzero value of nerr.
c
c Input: ax, ay, az, qc, qx, qy, qz, qxy, qyz, qzx,
c qxx, qyy, qzz, tol.
c
c Output: sc, sx, sy, sz, sxy, syz, szx, sxx, syy, szz, ntype, nerr.
c
c Calls: aptaxis, aptrois, aptrotv
c
c
c Glossary:
c
c ax,ay,az Input The x, y, z components of vector "a".
c
c nerr Output Indicates an input error, if not zero:
c 1 if the vector "a" is null.
c 2 if the coefficients qc, ..., qzz do not define
c a real quadric surface, or if there is no real
c projected outline.
c
c ntype Output The type of projected cylindrical quadric surface:
c The following are acceptable (nerr = 0):
c 0 = Plane.
c 1 = Coincident planes.
c 2 = Real parallel planes.
c 3 = Real intersecting planes.
c 4 = Parabolic cylinder.
c 5 = Hyperbolic cylinder.
c 6 = Real elliptic cylinder.
c 7 = Real circular cylinder.
c The following are not acceptable (nerr = 2):
c 17 = Imaginary parallel planes.
c 18 = Imaginary intersecting planes (x = y = 0).
c 19 = Imaginary elliptic cylinder.
c 20 = Imaginary circular cylinder.
c
c qc, ... Input Coefficients of the implicit equation for the quadric
c surface:
c qc + qx * x + qy * y + qz * z +
c qxy * x * y + qyz * y * z + qzx * z * x +
c qxx * x**2 + qyy * y**2 + qxx * z**2 = 0
c
c sc, ... Output Coefficients of the implicit equation for the projected
c cylindrical quadric surface with its cylindrical axis
c parallel to vector "a", and tangent to the original
c quadric surface:
c sc + sx * x + sy * y + sz * z +
c sxy * x * y + syz * y * z + szx * z * x +
c sxx * x**2 + syy * y**2 + sxx * z**2 = 0
c
c tol Input Numerical tolerance limit. With 64-bit floating point
c numbers, 1.e-5 to 1.e-11 is recommended.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Local arrays.
dimension rotm(3,3) ! Rotation tensor, "a" -> (0,0,1)
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptqupr finding projection "a"',
cbug & ' of a quadric surface.',
cbug & ' tol=',1pe13.5)
cbug 9902 format (
cbug & ' ax, ay, az= ',1p3e22.14 )
cbug 9903 format (
cbug & ' Quadric 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) ax, ay, az
cbug write ( 3, 9903) qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz
cbugc***DEBUG ends.
c.... Initialize.
sc = -1.e99
sx = -1.e99
sy = -1.e99
sz = -1.e99
sxy = -1.e99
syz = -1.e99
szx = -1.e99
sxx = -1.e99
syy = -1.e99
szz = -1.e99
c.... Find the rotation tensor to rotate vector "a" to the z axis.
call aptrotv (ax, ay, az, 0.0, 0.0, 1.0, tol, rotm, nerr)
if (nerr .ne. 0) go to 999
c.... Rotate the original quadric surface to make vector "a" a z-vector.
tc = qc
tx = qx
ty = qy
tz = qz
txy = qxy
tyz = qyz
tzx = qzx
txx = qxx
tyy = qyy
tzz = qzz
call aptrois (rotm, 0, 0.0, 0.0, 0.0,
& tc, tx, ty, tz, txy, tyz, tzx, txx, tyy, tzz, tol,
& nerr)
c.... The projection of the temporary quadric surface tc, ..., tzz in the
c.... z direction contains the intersection of the plane NZ = 0 with the
c.... temporary quadric surface.
c.... NZ = tz + tzx * x + tyz * y + 2.0 * tzz * z = 0
c.... z = -0.5 * (tz + tzx * x + tyz * y) / tzz (tzz not zero).
c.... See if any projected outline exists on any z plane.
if ((tz .eq. 0.0) .and.
& (tzx .eq. 0.0) .and.
& (tyz .eq. 0.0)) then
sc = tc
sx = tx
sy = ty
sz = 0.0
sxy = txy
syz = 0.0
szx = 0.0
sxx = txx
syy = tyy
szz = 0.0
elseif (tzz .ne. 0.0) then
sc = tc - 0.25 * tz**2 / tzz
sx = tx - 0.5 * tz * tzx / tzz
sy = ty - 0.5 * tz * tyz / tzz
sz = 0.0
sxy = txy - 0.5 * tzx * tyz / tzz
syz = 0.0
szx = 0.0
sxx = txx - 0.25 * tzx**2 / tzz
syy = tyy - 0.25 * tyz**2 / tzz
szz = 0.0
else ! No projection exists.
nerr = 2
go to 999
endif
c.... Rotate the projected quadric surface back to the original direction.
call aptrois (rotm, 1, 0.0, 0.0, 0.0,
& sc, sx, sy, sz, sxy, syz, szx, sxx, syy, szz, tol,
& nerr)
c.... Check to see if the projection is a real cylindrical quadric surface.
tc = sc
tx = sx
ty = sy
tz = sz
txy = sxy
tyz = syz
tzx = szx
txx = sxx
tyy = syy
tzz = szz
call aptaxis (tc, tx, ty, tz, txy, tyz, tzx, txx, tyy, tzz, tol,
& ax, ay, az, rotm, ntype, nerr)
if ((nerr .ne. 0) .or.
& (ntype .lt. 0) .or.
& (ntype .gt. 7)) then
nerr = 2
endif
999 continue
cbugc***DEBUG begins.
cbug 9920 format ('aptqupr results. nerr=',i2,' ntype=',i2)
cbug 9940 format ('Quadric cylinder projection in direction "a":' /
cbug & ' Quadric sc= ',1pe20.12 /
cbug & ' sx,sy,sz= ',1p3e20.12 /
cbug & ' sxy,syz,szx= ',1p3e20.12 /
cbug & ' sxx,syy,szz= ',1p3e20.12)
cbug write ( 3, 9920) nerr, ntype
cbug write ( 3, 9940) sc, sx, sy, sz, sxy, syz, szx, sxx, syy, szz
cbugc***DEBUG ends.
return
c.... End of subroutine aptqupr. (+1 line)
end
UCRL-WEB-209832