subroutine apttris (nopt, ax, ay, az, qc, qx, qy, qz,
& qxy, qyz, qzx, qxx, qyy, qzz, tol, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTTRIS
c
c call apttris (nopt, ax, ay, az, qc, qx, qy, qz,
c & qxy, qyz, qzx, qxx, qyy, qzz, tol, nerr)
c
c Version: apttris Updated 1999 March 15 12:15.
c apttris Originated 1990 October 10 18:00.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To find the new coefficients of the equation for a general
c second-order surface in xyz space, resulting from a
c translation of the point a = (ax, ay, az) to the origin
c (nopt = 0), or the origin to point "a" (nopt = 1).
c
c The equation for the surface is initially as follows:
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 translation is equivalent to the change of coordinates:
c
c x' = x - ax, y' = y - ay, z' = z - az (nopt = 0), or
c x' = x + ax, y' = y + ay, z' = z + az (nopt = 1).
c
c The equation for the translated surface is as follows:
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 Input: nopt, ax, ay, az, qc, qx, qy, qz, qxy, qyz, qzx,
c qxx, qyy, qzz, tol.
c
c Output: qc, qx, qy, qz, nerr.
c
c Glossary:
c
c ax,ay,az Input The x, y, z coordinates of point "a", or the
c translation vector "a".
c
c nerr Output Indicates an input error, if not zero:
c 1 if nopt is not 0 or 1.
c
c nopt Input The translation option:
c 0 to translate point "a" to the origin (subtract
c the coordinates of point "a" from all coordinates).
c 1 to translate the origin to point "a" (add
c the coordinates of point "a" to all coordinates).
c
c qc In/Out Constant term in implicit surface equation.
c Changed by a translation.
c Zero only if the surface passes through the origin.
c qc' = qc + f(ax,ay,az) (nopt = 0), or
c qc' = qc + f(-ax,-ay,-az) (nopt = 1).
c
c qx,qy,qz In/Out Coefficients of x, y, z in implicit surface equation.
c May be changed by a translation.
c All zero only if the center of symmetry of the
c surface is at the origin.
c (qx',qy',qz') = (qx,qy,qz) + J * (ax, ay, az), or
c (qx',qy',qz') = (qx,qy,qz) - J * (ax, ay, az),
c for nopt = 0 or 1, respectively, where J is the
c 3 by 3 matrix of second derivatives of f(x,y,z).
c
c q.. Input Coefficients of second-order terms in implicit
c surface equation: qxy, qyz, qzx, qxx, qyy, qzz.
c Coefficients of x*y, y*z, z*x, x**2, y**2, z**2,
c respectively. Invariant to translation of the
c origin. All zero for a planar surface.
c The coefficients qxy, qyz, qzx are all zero only if
c the symmetry axes of the surface are aligned with the
c x, y and z axes, or the surface is a plane.
c
c tol Input Numerical tolerance limit.
c On computers with 64-bit floating point numbers,
c recommend 1.e-5 to 1.e-11.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments. (none.)
c.... Local variables.
c---- Square of length of vector "a".
common /lapttris/ dasq
c---- Sum of squares of coeff.
common /lapttris/ dtsq
c---- Sum of squares of coeff.
common /lapttris/ dxsq
c---- Sum of squares of coeff.
common /lapttris/ dysq
c---- Sum of squares of coeff.
common /lapttris/ dzsq
c---- Temporary value of qc.
common /lapttris/ qcp
c---- Estimated error in qc.
common /lapttris/ qcperr
c---- Temporary value of qx.
common /lapttris/ qxp
c---- Estimated error in qx.
common /lapttris/ qxperr
c---- Temporary value of qy.
common /lapttris/ qyp
c---- Estimated error in qy.
common /lapttris/ qyperr
c---- Temporary value of qz.
common /lapttris/ qzp
c---- Estimated error in qz.
common /lapttris/ qzperr
cbugc***DEBUG begins.
cbug 9901 format (/ 'apttris translating an implicit surface. nopt=',i2,
cbug & ' tol=',1pe13.5 /
cbug & ' ax,ay,az= ',1p3e22.14 /
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) nopt, tol, ax, ay, az, qc, qx, qy, qz,
cbug & qxy, qyz, qzx, qxx, qyy, qzz
cbugc***DEBUG ends.
c.... Initialize.
nerr = 0
dasq = ax**2 + ay**2 + az**2
c.... Test for input errors.
if ((nopt .lt. 0) .or. (nopt .gt. 1)) then
nerr = 1
cbugc***DEBUG begins.
cbug 9902 format (/ "nopt is not in the range from 0 to 1")
cbug write ( 3, 9902)
cbugc***DEBUG ends.
go to 210
endif
c.... Find the new coefficients of the implicit second-order surface.
c---- Subtract "a" from all coordinates.
if (nopt .eq. 0) then
if (dasq .le. tol**2) go to 210
qcp = qc + qx * ax + qy * ay + qz * az +
& qxy * ax * ay + qyz * ay * az + qzx * az * ax +
& qxx * ax**2 + qyy * ay**2 + qzz * az**2
qxp = qx + qxy * ay + qzx * az + 2.0 * qxx * ax
qyp = qy + qyz * az + qxy * ax + 2.0 * qyy * ay
qzp = qz + qzx * ax + qyz * ay + 2.0 * qzz * az
c---- Add "a" to all coordinates.
elseif (nopt .eq. 1) then
if (dasq .le. tol**2) go to 210
qcp = qc - qx * ax - qy * ay - qz * az +
& qxy * ax * ay + qyz * ay * az + qzx * az * ax +
& qxx * ax**2 + qyy * ay**2 + qzz * az**2
qxp = qx - qxy * ay - qzx * az - 2.0 * qxx * ax
qyp = qy - qyz * az - qxy * ax - 2.0 * qyy * ay
qzp = qz - qzx * ax - qyz * ay - 2.0 * qzz * az
c---- Tested nopt.
endif
c.... See if small values are to be truncated to zero.
c---- Test for truncation to zero.
if (tol .gt. 0.0) then
qcperr = tol * (abs (qc) +
& 2.0 * (abs (qx * ax) + abs (qy * ay) + abs (qz * az)) +
& 3.0 * (abs (qxy * ax * ay) +
& abs (qyz * ay * az) +
& abs (qzx * az * ax) +
& abs (qxx * ax**2) +
& abs (qyy * ay**2) +
& abs (qzz * az**2)))
qxperr = tol * (abs (qx) +
& 2.0 * (abs (qxy * ay) + abs (qzx * az) +
& abs (2.0 * qxx * ax)))
qyperr = tol * (abs (qy) +
& 2.0 * (abs (qyz * az) + abs (qxy * ax) +
& abs (2.0 * qyy * ay)))
qzperr = tol * (abs (qz) +
& 2.0 * (abs (qzx * ax) + abs (qyz * ay) +
& abs (2.0 * qzz * az)))
if (abs (qcp) .lt. qcperr) then
qcp = 0.0
endif
if (abs (qxp) .lt. qxperr) then
qxp = 0.0
endif
if (abs (qyp) .lt. qyperr) then
qyp = 0.0
endif
if (abs (qzp) .lt. qzperr) then
qzp = 0.0
endif
c---- Tested tol.
endif
c.... Store the new coefficients.
qc = qcp
qx = qxp
qy = qyp
qz = qzp
cbugc***DEBUG begins.
cbug 9903 format (/ 'apttris results:' /
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, 9903) qc, qx, qy, qz,
cbug & qxy, qyz, qzx, qxx, qyy, qzz
cbugc***DEBUG ends.
210 return
c.... End of subroutine apttris. (+1 line.)
end
UCRL-WEB-209832