subroutine aptrois (smat, inv, ax, ay, az, qc, qx, qy, qz,
& qxy, qyz, qzx, qxx, qyy, qzz, tol, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTROIS
c
c call aptrois (smat, inv, ax, ay, az, qc, qx, qy, qz,
c & qxy, qyz, qzx, qxx, qyy, qzz, tol, nerr)
c
c Version: aptrois Updated 1999 March 15 11:30.
c aptrois Originated 1990 October 16 15: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 general
c coordinate transformation (rotation, reflection, and/or
c inversion) around the invariant point a = (ax, ay, az),
c based on the transformation matrix smat (inv = 0) or its
c inverse (inv = 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 After the transformation, replace x, y, z with
c
c x' = smat(1,1) * x + smat(1,2) * y + smat(1,3) * z
c y' = smat(2,1) * x + smat(2,2) * y + smat(2,3) * z
c z' = smat(3,1) * x + smat(3,2) * y + smat(3,3) * z (inv = 0)
c
c x' = smat(1,1) * x + smat(2,1) * y + smat(3,1) * z
c y' = smat(1,2) * x + smat(2,2) * y + smat(3,2) * z
c z' = smat(1,3) * x + smat(2,3) * y + smat(3,3) * z (inv = 1)
c
c Note: use subroutine aptmopv to transform points and vectors
c with matrix smat.
c
c Input: smat, inv, ax, ay, az, qc, qx, qy, qz, qxy, qyz, qzx,
c qxx, qyy, qzz, tol.
c
c Output: qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz, nerr.
c
c Calls: apttris
c
c Glossary:
c
c ax,ay,az Input The x, y, z coordinates of an invariant point "a".
c
c inv Input The transformation option:
c 0 to use the transformation matrix smat.
c 1 to use the inverse of the transformation matrix
c smat (exchange rows and columns).
c
c nerr Output Indicates an input error, if not zero:
c 1 if inv is not 0 or 1.
c
c qc In/Out Constant term in implicit surface equation.
c Zero only if surface passes through origin.
c Invariant to a general symmetry transformation.
c
c qx,qy,qz In/Out Coefficients of x, y, z in implicit surface equation.
c All zero only if the center of symmetry of the
c surface is at the origin.
c
c q.. In/Out Coefficients of the 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. 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 The sum qxx + qyy + qzz is invariant to a general
c symmetry transformation.
c
c smat Input The coordinate transformation matrix operator.
c It transforms a point (x, y, z) to the point
c (x', y', z') as follows:
c x' = smat(1,1) * x + smat(1,2) * y + smat(1,3) * z
c y' = smat(2,1) * x + smat(2,2) * y + smat(2,3) * z
c z' = smat(3,1) * x + smat(3,2) * y + smat(3,3) * z.
c The rows and columns of smat must be unit vectors,
c with each row orthogonal to the other rows, and each
c column orthogonal to the other columns. For the
c inverse operation, the rows and columns of smat
c are interchanged. See the APT subroutines with
c names beginning with aptinv, aptref, aptrot.
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 transformation matrix.
dimension smat (3,3)
c.... Local variables.
c---- Coefficient of x in equation for x.
common /laptrois/ cxx
c---- Coefficient of y in equation for x.
common /laptrois/ cxy
c---- Coefficient of z in equation for x.
common /laptrois/ cxz
c---- Coefficient of x in equation for y.
common /laptrois/ cyx
c---- Coefficient of y in equation for y.
common /laptrois/ cyy
c---- Coefficient of z in equation for y.
common /laptrois/ cyz
c---- Coefficient of x in equation for z.
common /laptrois/ czx
c---- Coefficient of y in equation for z.
common /laptrois/ czy
c---- Coefficient of z in equation for z.
common /laptrois/ czz
c---- Distance of point "a" from origin.
common /laptrois/ da
c---- Temporary value of qc.
common /laptrois/ qcp
c---- Estimated error in qc.
common /laptrois/ qcperr
c---- Temporary value of qx.
common /laptrois/ qxp
c---- Estimated error in qx.
common /laptrois/ qxperr
c---- Temporary value of qxx.
common /laptrois/ qxxp
c---- Estimated error in qxx.
common /laptrois/ qxxperr
c---- Temporary value of qxy.
common /laptrois/ qxyp
c---- Estimated error in qxy.
common /laptrois/ qxyperr
c---- Temporary value of qy.
common /laptrois/ qyp
c---- Estimated error in qy.
common /laptrois/ qyperr
c---- Temporary value of qyy.
common /laptrois/ qyyp
c---- Estimated error in qyy.
common /laptrois/ qyyperr
c---- Temporary value of qyz.
common /laptrois/ qyzp
c---- Estimated error in qyz.
common /laptrois/ qyzperr
c---- Temporary value of qz.
common /laptrois/ qzp
c---- Estimated error in qz.
common /laptrois/ qzperr
c---- Temporary value of qzx.
common /laptrois/ qzxp
c---- Estimated error in qzx.
common /laptrois/ qzxperr
c---- Temporary value of qzz.
common /laptrois/ qzzp
c---- Estimated error in qzz.
common /laptrois/ qzzperr
cbugc***DEBUG begins.
cbugc---- Index in smat.
cbug common /laptrois/ i
cbugc---- Index in smat.
cbug common /laptrois/ j
cbug 9901 format (/ 'aptrois transforming an implicit surface. inv=',i2,
cbug & ' tol=',1pe13.5)
cbug 9902 format (' smat=',3(/ 14x,1p3e22.14))
cbug 9903 format (/
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) inv, tol
cbug write ( 3, 9902) ((smat(i,j), j = 1, 3), i = 1, 3)
cbug write ( 3, 9903) ax, ay, az, qc, qx, qy, qz,
cbug & qxy, qyz, qzx, qxx, qyy, qzz
cbugc***DEBUG ends.
c.... Initialize.
nerr = 0
c.... Test for input errors.
if ((inv .lt. 0) .or. (inv .gt. 1)) then
nerr = 1
cbugcbugc***DEBUG begins.
cbugcbug 9904 format (/ "inv is not 0 or 1")
cbugcbug write ( 3, 9904)
cbugcbugc***DEBUG ends.
go to 210
endif
c.... Find the distance of point "a" from the origin.
da = sqrt (ax**2 + ay**2 + az**2)
c.... Translate the invariant point "a" to the origin.
c---- Translate point "a" to origin.
if (da .gt. tol) then
call apttris (0, ax, ay, az, qc, qx, qy, qz, qxy, qyz, qzx,
& qxx, qyy, qzz, tol, nerr)
cbugcbugc***DEBUG begins.
cbugcbug 9801 format (/ 'aptrois DEBUG 1:' /
cbugcbug & ' qc= ',1pe22.14 /
cbugcbug & ' qx,qy,qz= ',1p3e22.14 /
cbugcbug & ' qxy,qyz,qzx=',1p3e22.14 /
cbugcbug & ' qxx,qyy,qzz=',1p3e22.14)
cbugcbug write ( 3, 9801) qc, qx, qy, qz,
cbugcbug & qxy, qyz, qzx, qxx, qyy, qzz
cbugcbugc***DEBUG ends.
c---- Tested da.
endif
c.... Find the transformation coefficients.
c---- Operate with matrix smat.
if (inv .eq. 0) then
cxx = smat(1,1)
cxy = smat(1,2)
cxz = smat(1,3)
cyx = smat(2,1)
cyy = smat(2,2)
cyz = smat(2,3)
czx = smat(3,1)
czy = smat(3,2)
czz = smat(3,3)
c---- Operate with the inverse of smat.
elseif (inv .eq. 1) then
cxx = smat(1,1)
cxy = smat(2,1)
cxz = smat(3,1)
cyx = smat(1,2)
cyy = smat(2,2)
cyz = smat(3,2)
czx = smat(1,3)
czy = smat(2,3)
czz = smat(3,3)
c---- Tested inv.
endif
c.... Find the new implicit equation coefficients.
qcp = qc
qxp = cxx * qx + cxy * qy + cxz * qz
qyp = cyy * qy + cyz * qz + cyx * qx
qzp = czz * qz + czx * qx + czy * qy
qxyp = (cxx * cyy + cxy * cyx) * qxy +
& (cxy * cyz + cxz * cyy) * qyz +
& (cxz * cyx + cxx * cyz) * qzx +
& 2.0 * (cxx * cyx * qxx + cxy * cyy * qyy + cxz * cyz * qzz)
qyzp = (cyy * czz + cyz * czy) * qyz +
& (cyz * czx + cyx * czz) * qzx +
& (cyx * czy + cyy * czx) * qxy +
& 2.0 * (cyy * czy * qyy + cyz * czz * qzz + cyx * czx * qxx)
qzxp = (czz * cxx + czx * cxz) * qzx +
& (czx * cxy + czy * cxx) * qxy +
& (czy * cxz + czz * cxy) * qyz +
& 2.0 * (czz * cxz * qzz + czx * cxx * qxx + czy * cxy * qyy)
qxxp = cxx * cxy * qxy + cxy * cxz * qyz + cxz * cxx * qzx +
& cxx**2 * qxx + cxy**2 * qyy + cxz**2 * qzz
qyyp = cyy * cyz * qyz + cyz * cyx * qzx + cyx * cyy * qxy +
& cyy**2 * qyy + cyz**2 * qzz + cyx**2 * qxx
qzzp = czz * czx * qzx + czx * czy * qxy + czy * czz * qyz +
& czz**2 * qzz + czx**2 * qxx + czy**2 * qyy
cbugcbugc***DEBUG begins.
cbugcbug 9802 format (/ 'aptrois DEBUG 2:' /
cbugcbug & ' qc= (all p) ',1pe22.14 /
cbugcbug & ' qx,qy,qz= ',1p3e22.14 /
cbugcbug & ' qxy,qyz,qzx=',1p3e22.14 /
cbugcbug & ' qxx,qyy,qzz=',1p3e22.14)
cbugcbug write ( 3, 9802) qcp, qxp, qyp, qzp,
cbugcbug & qxyp, qyzp, qzxp, qxxp, qyyp, qzzp
cbugcbugc***DEBUG ends.
c.... See if small values are to be truncated to zero.
c---- Test for truncation to zero.
if (tol .gt. 0.0) then
qcperr = (abs (qc)) * tol
qxperr = (abs (cxx * qx) + abs (cxy * qy) + abs (cxz * qz)) *
& 2.0 * tol
qyperr = (abs (cyy * qy) + abs (cyz * qz) + abs (cyx * qx)) *
& 2.0 * tol
qzperr = (abs (czz * qz) + abs (czx * qx) + abs (czy * qy)) *
& 2.0 * tol
qxyperr = (abs (cxx * cyy * qxy) + abs (cxy * cyx * qxy) +
& abs (cxy * cyz * qyz) + abs (cxz * cyy * qyz) +
& abs (cxz * cyx * qzx) + abs (cxx * cyz * qzx) +
& 2.0 * (abs (cxx * cyx * qxx) +
& abs (cxy * cyy * qyy) +
& abs (cxz * cyz * qzz))) * 3.0 * tol
qyzperr = (abs (cyy * czz * qyz) + abs (cyz * czy * qyz) +
& abs (cyz * czx * qzx) + abs (cyx * czz * qzx) +
& abs (cyx * czy * qxy) + abs (cyy * czx * qxy) +
& 2.0 * (abs (cyy * czy * qyy) +
& abs (cyz * czz * qzz) +
& abs (cyx * czx * qxx))) * 3.0 * tol
qzxperr = (abs (czz * cxx * qzx) + abs (czx * cxz * qzx) +
& abs (czx * cxy * qxy) + abs (czy * cxx * qxy) +
& abs (czy * cxz * qyz) + abs (czz * cxy * qyz) +
& 2.0 * (abs (czz * cxz * qzz) +
& abs (czx * cxx * qxx) +
& abs (czy * cxy * qyy))) * 3.0 * tol
qxxperr = (abs (cxx * cxy * qxy) +
& abs (cxy * cxz * qyz) +
& abs (cxz * cxx * qzx) +
& abs (cxx**2 * qxx) +
& abs (cxy**2 * qyy) +
& abs (cxz**2 * qzz)) * 3.0 * tol
qyyperr = (abs (cyy * cyz * qyz) +
& abs (cyz * cyx * qzx) +
& abs (cyx * cyy * qxy) +
& abs (cyy**2 * qyy) +
& abs (cyz**2 * qzz) +
& abs (cyx**2 * qxx)) * 3.0 * tol
qzzperr = (abs (czz * czx * qzx) +
& abs (czx * czy * qxy) +
& abs (czy * czz * qyz) +
& abs (czz**2 * qzz) +
& abs (czx**2 * qxx) +
& abs (czy**2 * qyy)) * 3.0 * tol
cbugcbugc***DEBUG begins.
cbugcbug 9803 format (/ 'aptrois DEBUG 3:' /
cbugcbug & ' qc=(all err)',1pe22.14 /
cbugcbug & ' qx,qy,qz= ',1p3e22.14 /
cbugcbug & ' qxy,qyz,qzx=',1p3e22.14 /
cbugcbug & ' qxx,qyy,qzz=',1p3e22.14)
cbugcbug write ( 3, 9803) qcperr, qxperr, qyperr, qzperr,
cbugcbug & qxyperr, qyzperr, qzxperr, qxxperr, qyyperr, qzzperr
cbugcbugc***DEBUG ends.
if (abs (qcp) .le. qcperr) then
qcp = 0.0
endif
if (abs (qxp) .le. qxperr) then
qxp = 0.0
endif
if (abs (qyp) .le. qyperr) then
qyp = 0.0
endif
if (abs (qzp) .le. qzperr) then
qzp = 0.0
endif
if (abs (qxyp) .le. qxyperr) then
qxyp = 0.0
endif
if (abs (qyzp) .le. qyzperr) then
qyzp = 0.0
endif
if (abs (qzxp) .le. qzxperr) then
qzxp = 0.0
endif
if (abs (qxxp) .le. qxxperr) then
qxxp = 0.0
endif
if (abs (qyyp) .le. qyyperr) then
qyyp = 0.0
endif
if (abs (qzzp) .le. qzzperr) then
qzzp = 0.0
endif
c---- Tested tol.
endif
c.... Store the new coefficients.
qc = qcp
qx = qxp
qy = qyp
qz = qzp
qxy = qxyp
qyz = qyzp
qzx = qzxp
qxx = qxxp
qyy = qyyp
qzz = qzzp
cbugcbugc***DEBUG begins.
cbugcbug 9804 format (/ 'aptrois DEBUG 4:' /
cbugcbug & ' qc= ',1pe22.14 /
cbugcbug & ' qx,qy,qz= ',1p3e22.14 /
cbugcbug & ' qxy,qyz,qzx=',1p3e22.14 /
cbugcbug & ' qxx,qyy,qzz=',1p3e22.14)
cbugcbug write ( 3, 9804) qc, qx, qy, qz,
cbugcbug & qxy, qyz, qzx, qxx, qyy, qzz
cbugcbugc***DEBUG ends.
c.... Translate the origin back to the invariant point "a".
c---- Translate origin to point "a".
if (da .gt. tol) then
call apttris (1, ax, ay, az, qc, qx, qy, qz, qxy, qyz, qzx,
& qxx, qyy, qzz, tol, nerr)
cbugcbugc***DEBUG begins.
cbugcbug 9805 format (/ 'aptrois DEBUG 5:' /
cbugcbug & ' qc= ',1pe22.14 /
cbugcbug & ' qx,qy,qz= ',1p3e22.14 /
cbugcbug & ' qxy,qyz,qzx=',1p3e22.14 /
cbugcbug & ' qxx,qyy,qzz=',1p3e22.14)
cbugcbug write ( 3, 9805) qc, qx, qy, qz,
cbugcbug & qxy, qyz, qzx, qxx, qyy, qzz
cbugcbugc***DEBUG ends.
c---- Tested da.
endif
c.... Truncate small coefficients.
psq = qx**2 + qy**2 + qz**2
qsq = qxy**2 + qyz**2 + qzx**2
ssq = qxx**2 + qyy**2 + qzz**2
if (qx**2 .le. tol**2 * psq) qx = 0.0
if (qy**2 .le. tol**2 * psq) qy = 0.0
if (qx**2 .le. tol**2 * psq) qx = 0.0
if (qxy**2 .le. tol**2 * qsq) qxy = 0.0
if (qyz**2 .le. tol**2 * qsq) qyz = 0.0
if (qzx**2 .le. tol**2 * qsq) qzx = 0.0
if (qxx**2 .le. tol**2 * ssq) qxx = 0.0
if (qyy**2 .le. tol**2 * ssq) qyy = 0.0
if (qzz**2 .le. tol**2 * ssq) qzz = 0.0
cbugc***DEBUG begins.
cbug 9905 format (/ 'aptrois 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, 9905) qc, qx, qy, qz,
cbug & qxy, qyz, qzx, qxx, qyy, qzz
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptrois. (+1 line.)
end
UCRL-WEB-209832