subroutine aptaxis (qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz,
& tol, ax, ay, az, rotm, ntype, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTAXIS
c
c call aptaxis (qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz,
c & tol, ax, ay, az, rotm, ntype, nerr)
c
c Version: aptaxis Updated 1999 March 3 17:40.
c aptaxis Originated 1990 November 1 16:00.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: Given the coefficients of the implicit equation for a general
c quadric surface:
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 to find the surface type ntype, the coefficients of the
c standard form of the surface, the 3 by 3 rotation matrix rotm
c which aligns the axes of the surface with the x, y and z axes,
c and the translation vector a = (ax, ay, az) which translates the
c central point of the surface to the origin.
c
c After the transformation, replace x, y, z with
c
c x' = rotm(1,1)*(x-ax) + rotm(1,2)*(y-ay) + rotm(1,3)*(z-az),
c y' = rotm(2,1)*(x-ax) + rotm(2,2)*(y-ay) + rotm(2,3)*(z-az),
c z' = rotm(3,1)*(x-ax) + rotm(3,2)*(y-ay) + rotm(3,3)*(z-az),
c
c and use the new coefficients.
c
c Note: use subroutine aptmopv to transform points and vectors
c with matrix rotm, subroutine apttran to translate points with
c vector "a".
c
c Input: qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz, tol.
c
c Output: qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz,
c ax, ay, az, rotm, nerr.
c
c Calls: aptmopv, aptmprd, aptrois, aptrota, aptrott, aptrotv, apttris,
c aptvaxb
c
c
c Glossary:
c
c ax,ay,az Output The initial center of the quadric surface.
c It may be used to translate an initial point
c (x, y, z) to the point (x', y', z') as
c follows:
c
c x' = x - ax, y' = y = ay, z' = z - az.
c
c nerr Output Indicates an input error, if not zero:
c 1 if all coefficients are zero except qc.
c
c ntype Output Indicates the classification of the quadric surface.
c The standard forms of each type are shown below
c (each coefficient must have the preceding sign).
c -1 = Input error, no classification.
c 0 = Plane. x = 0
c 1 = Coincident planes. x**2 = 0
c 2 = Real parallel planes. -1 + qxx*x**2 = 0
c 3 = Real intersecting planes. x**2 - |qyy|*y**2 = 0
c 4 = Parabolic cylinder. -|qy|*y + x**2 = 0
c 5 = Hyperbolic cylinder.
c 1 + qxx*x**2 - |qyy|*y**2 = 0
c 6 = Real elliptic cylinder.
c -1 + qxx*x**2 + qyy*y**2 = 0
c 7 = Real circular cylinder.
c -1 + qxx * (x**2 + y**2) = 0
c 8 = Hyperbolic paraboloid.
c -|qz|*z + x**2 - |qyy|*y**2 = 0
c 9 = Elliptic paraboloid.
c -|qz|*z + x**2 + qyy*y**2 = 0
c 10 = Circular paraboloid.
c -|qz|*z + x**2 + y**2 = 0
c 11 = Real elliptic cone.
c x**2 + qyy*y**2 - |qzz|*z**2 = 0
c 12 = Real circular cone.
c x**2 + y**2 - |qzz|*z**2 = 0
c 13 = Hyperboloid of one sheet.
c -1 + qxx*x**2 + qyy*y**2 - |qzz|*z**2 = 0
c (circular if qxx = qyy).
c 14 = Hyperboloid of two sheets.
c 1 + qxx*x**2 + qyy*y**2 - |qzz|*z**2 = 0
c (circular if qxx = qyy).
c 15 = Real ellipsoid.
c -1 + qxx*x**2 + qyy*y**2 + qzz*z**2 = 0
c (prolate spheroid if qxx = qyy,
c oblate spheroid if qyy = qzz).
c 16 = Real sphere.
c -1 + qxx * (x**2 + y**2 + z**2) = 0
c 17 = Imaginary parallel planes. 1 + qxx*x**2 = 0
c 18 = Imaginary intersecting planes (x = y = 0).
c x**2 + qyy*y**2 = 0
c 19 = Imaginary elliptic cylinder.
c 1 + qxx*x**2 + qyy*y**2 = 0
c 20 = Imaginary circular cylinder.
c 1 + qxx * (x**2 + y**2) = 0
c 21 = Imaginary elliptic cone (x = y = z = 0).
c x**2 + qyy*y**2 + qzz*z**2 = 0
c 22 = Imaginary circular cone (x = y = z = 0).
c x**2 + y**2 + qzz*z**2 = 0
c 23 = Imaginary ellipsoid.
c 1 + qxx*x**2 + qyy*y**2 + qzz*z**2 = 0
c 24 = Imaginary sphere.
c 1 + qxx * (x**2 + y**2 + z**2) = 0
c
c qc In/Out Constant term in the implicit surface equation.
c Zero only if the surface passes through the origin.
c
c qx,qy,qz In/Out Coefficients of x, y, z in implicit surface equation.
c For a plane in standard form, only qx is non-zero.
c For a second-order quadric surface in standard form,
c qx and qy are zero, but if qzz is zero, qz may be
c non-zero.
c
c q.. In/Out Coefficients of the second-order terms in the 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. For a plane, all are zero.
c For a second-order quadric surface in standard form,
c qxy, qyz and qzx are zero, and qxx, qyy and qzz are
c ordered in equal or decreasing order of absolute
c value, with qxx positive.
c
c rotm Output The rotation matrix operator which aligns the axes of
c the surface with the x, y and z axes. It may be
c used to rotate a point (x, y, z) to the point
c (x', y', z') as follows:
c
c x' = rotm(1,1) * x + rotm(1,2) * y + rotm(1,3) * z,
c y' = rotm(2,1) * x + rotm(2,2) * y + rotm(2,3) * z,
c z' = rotm(3,1) * x + rotm(3,2) * y + rotm(3,3) * z.
c
c Its row vectors are the initial unit vectors along
c the principal axes of the surface, which are mutually
c perpendicular. Two may be arbitrary.
c Size 3 by 3.
c
c tol Input Numerical tolerance limit. If zero, 1.e-11 will be
c used instead.
c On Cray computers, recommend 1.e-11.
c On 32-bit computers, recommend 1.e-6.
c
c History: 1997 June 11 14:40. Removed requirement that qc > 0
c for hyperbolic cylinder.
c 1997 September 17 12:40. Changed standard forms of
c hyperbolic and parabolic cylinders and elliptic and
c hyperbolic paraboloids.
c Eliminated an unnecessary rotation of 180 degrees
c around x axis.
c 1999 February 22 17:00. Made vector b in call to
c aptrotv a unit vector, in case otherwise too small.
c 1999 March 1 14:40. Made test for plane exact, no
c quadric coefficients allowed.
c 1990 March 3 17:40. Added truncation of center point.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
dimension rotm (3,3) ! Coordinate transformation matrix.
c.... Local variables.
common /laptaxis/ bx ! Component x of axis vector "b".
common /laptaxis/ by ! Component y of axis vector "b".
common /laptaxis/ bz ! Component z of axis vector "b".
common /laptaxis/ blen ! Length of axis vector "b".
common /laptaxis/ cx ! Component x of axis vector "c".
common /laptaxis/ cy ! Component y of axis vector "c".
common /laptaxis/ cz ! Component z of axis vector "c".
common /laptaxis/ clen ! Length of axis vector "c".
common /laptaxis/ det ! Invariant of surface.
common /laptaxis/ deterr ! Estimated error in det.
common /laptaxis/ dsum ! Invariant of surface.
common /laptaxis/ dsumerr ! Estimated error in dsum.
common /laptaxis/ eigen (3) ! Eigenvalue of characteristic matrix.
common /laptaxis/ fact ! Intermediate result.
common /laptaxis/ fuz ! A very small number.
common /laptaxis/ i ! Index in rotm, rotms.
common /laptaxis/ j ! Index in rotm, rotms.
common /laptaxis/ n ! Index in eigen.
common /laptaxis/ nswap ! Indicates axes swapped, if 1.
common /laptaxis/ pi ! Mathematical constant, pi.
common /laptaxis/ pp ! Factor in finding eigen.
common /laptaxis/ pperr ! Estimated error in pp.
common /laptaxis/ psq ! Sum qx**2 + qy**2 + qz**2.
common /laptaxis/ qcs ! Value of qc before centering.
common /laptaxis/ qd ! Factor in finding eigen.
common /laptaxis/ qderr ! Estimated error in qd.
common /laptaxis/ qq ! Factor in finding eigen.
common /laptaxis/ qqerr ! Estimated error in qq.
common /laptaxis/ qsq ! Sum qxy**2 + qyz**2 + qzx**2.
common /laptaxis/ rotms (3,3) ! Last cumulative rotation matrix.
common /laptaxis/ rotp (3,3) ! Rotation matrix.
common /laptaxis/ ssq ! Sum qxx**2 + qyy**2 + qzz**2.
common /laptaxis/ theta ! Factor in finding eigen.
common /laptaxis/ thetax ! Angle of rotation around x axis.
common /laptaxis/ thetay ! Angle of rotation around y ayis.
common /laptaxis/ thetaz ! Angle of rotation around z axis.
common /laptaxis/ trace ! Invariant 2.0 * (qxx + qyy + qzz).
common /laptaxis/ tracerr ! Estimated error in trace.
common /laptaxis/ tsq ! Sum of qc**2, psq, qsq, ssq.
cbugc***DEBUG begins.
cbug common /laptaxis/ dets ! Initial values of det.
cbug common /laptaxis/ dsums ! Initial values of dsum.
cbug common /laptaxis/ traces ! Initial values of trace.
cbugc***DEBUG ends.
tolx = tol
if (tol .eq. 0.0) tolx = 1.e-11
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptaxis aligning an implicit 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 9904 format (/ " Invariants: trace, dsum, det=" /
cbug & 2x,1p3e22.14," errors:" / 2x,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.
fuz = 1.e-99 ! A very small number.
pi = 3.14159265358979323 ! Mathematical constant, pi.
nerr = -1
ntype = -1
psq = qx**2 + qy**2 + qz**2
qsq = qxy**2 + qyz**2 + qzx**2
ssq = qxx**2 + qyy**2 + qzz**2
tsq = qc**2 + psq + qsq + ssq
cbugcbugc***DEBUG begins.
cbugcbug 9801 format (/ ' Sums of squares of coefficients:' /
cbugcbug & ' psq,qsq,ssq=',1p3e22.14 /
cbugcbug & ' tsq =',1pe22.14 )
cbugcbug write ( 3, 9801) psq, qsq, ssq, tsq
cbugcbugc***DEBUG ends.
c.... Test for an input error.
cold if ((psq + qsq + ssq) .le. tolx**2 * tsq) then ! Bad coefficients.
if ((psq + qsq + ssq) .eq. 0.0) then ! Bad coefficients.
nerr = 1
cbugcbugc***DEBUG begins.
cbugcbug 9903 format (/ " FATAL. All coefficients but qc are zero.",
cbugcbug & " qc=",1pe22.14)
cbugcbug write ( 3, 9903) qc
cbugcbugc***DEBUG ends.
go to 310
endif ! Tested for bad coefficients.
c.... Initialize the rotation matrix to the identity matrix.
do 110 i = 1, 3
do 105 j = 1, 3
rotm(i,j) = 0.0
105 continue
110 continue
rotm(1,1) = 1.0
rotm(2,2) = 1.0
rotm(3,3) = 1.0
c.... Find invariants of surface.
trace = 2.0 * (qxx + qyy + qzz)
dsum = 4.0 * (qxx * qyy + qyy * qzz + qzz * qxx) -
& (qxy**2 + qyz**2 + qzx**2)
det = 2.0 * (4.0 * qxx * qyy * qzz + qxy * qyz * qzx -
& (qxx * qyz**2 + qyy * qzx**2 + qzz * qxy**2))
if (tolx .gt. 0.0) then ! Truncate small values to zero.
tracerr = 2.0 * tolx * (abs (qxx) + abs (qyy) + abs (qzz))
dsumerr = 4.0 * tolx * (abs (qxx * qyy) +
& abs (qyy * qzz) +
& abs (qzz * qxx)) +
& tolx * (qxy**2 + qyz**2 + qzx**2)
deterr = 2.0 * tolx *
& (4.0 * abs (qxx * qyy * qzz) +
& abs (qxy * qyz * qzx) +
& qyz**2 * abs (qxx) +
& qzx**2 * abs (qyy) +
& qxy**2 * abs (qzz))
if (abs (trace) .le. tracerr) trace = 0.0
if (abs (dsum) .le. dsumerr) dsum = 0.0
if (abs (det) .le. deterr) det = 0.0
endif ! Tested tolx.
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9904) trace, dsum, det, tracerr, dsumerr, deterr
cbugcbugc***DEBUG ends.
c=======================================================================--------
c.... See if the surface is a plane.
c++++ Surface is a plane.
cold if ((qsq + ssq) .le. tolx**2 * tsq) then
if ((trace .eq. 0.0) .and.
& (dsum .eq. 0.0) .and.
& (det .eq. 0.0)) then
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, '(/ " aptaxis: surface is a plane.")')
cbugcbugc***DEBUG ends.
ntype = 0 ! Plane.
c.... Use the point in the plane nearest the origin for the center
c.... of the standard form.
ax = - qx * qc / (psq + fuz)
ay = - qy * qc / (psq + fuz)
az = - qz * qc / (psq + fuz)
c.... Rotate the surface to be perpendicular to the x axis, and
c.... translate the central point to the origin.
call aptrotv (qx, qy, qz, 1., 0., 0., tolx, rotm, nerr)
qc = 0.0
qx = sqrt (psq)
qy = 0.0
qz = 0.0
go to 310
endif ! Surface is a plane.
c=======================================================================--------
c.... Surface is quadric. See if any mix terms (qxy, qyz, qzx) are non-zero.
cbugcbugc***DEBUG begins.
cbugcbug 9802 format (/ ' Test for coefficients qxy, qyz, qzx.' )
cbugcbug write ( 3, 9802)
cbugcbugc***DEBUG ends.
c++++ Non-zero mix coefficients.
if (qsq .gt. tolx**2 * ssq) then
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, '(/ " aptaxis: non-zero mix coefficients.")')
cbugcbug write ( 3, 9908) nerr, ntype, qc, qx, qy, qz,
cbugcbug & qxy, qyz, qzx, qxx, qyy, qzz
cbugcbugc***DEBUG ends.
if ((qzx**2 .le. tolx**2 * qsq) .and.
& (qxy**2 .le. tolx**2 * qsq)) then
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, '(/ " aptaxis: only qyz non-zero.")')
cbugcbugc***DEBUG ends.
c.... Only qyz non-zero. Make qyz zero by rotation around the x axis.
thetax = 0.5 * atan2 (qyz, qzz - qyy)
call aptrota (thetax, 1, 1., 0., 0., tolx, rotm, nerr)
call aptrois (rotm, 0, 0., 0., 0., qc, qx, qy, qz,
& qxy, qyz, qzx, qxx, qyy, qzz, tolx, nerr)
qyz = 0.0
elseif ((qxy**2 .le. tolx**2 * qsq) .and.
& (qyz**2 .le. tolx**2 * qsq)) then
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, '(/ " aptaxis: only qzx non-zero.")')
cbugcbug write ( 3, 9908) nerr, ntype, qc, qx, qy, qz,
cbugcbug & qxy, qyz, qzx, qxx, qyy, qzz
cbugcbugc***DEBUG ends.
c.... Only qzx non-zero. Make qzx zero by rotation around the y axis.
thetay = 0.5 * atan2 (qzx, qxx - qzz)
call aptrota (thetay, 1, 0., 1., 0., tolx, rotm, nerr)
call aptrois (rotm, 0, 0., 0., 0., qc, qx, qy, qz,
& qxy, qyz, qzx, qxx, qyy, qzz, tolx, nerr)
qzx = 0.0
elseif ((qyz**2 .le. tolx**2 * qsq) .and.
& (qzx**2 .le. tolx**2 * qsq)) then
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, '(/ " aptaxis: only qxy non-zero.")')
cbugcbug write ( 3, 9908) nerr, ntype, qc, qx, qy, qz,
cbugcbug & qxy, qyz, qzx, qxx, qyy, qzz
cbugcbugc***DEBUG ends.
c.... Only qxy non-zero. Make qxy zero by rotation around the z axis.
thetaz = 0.5 * atan2 (qxy, qyy - qxx)
call aptrota (thetaz, 1, 0., 0., 1., tolx, rotm, nerr)
call aptrois (rotm, 0, 0., 0., 0., qc, qx, qy, qz,
& qxy, qyz, qzx, qxx, qyy, qzz, tolx, nerr)
qxy = 0.0
c=======================================================================--------
else ! Two or more mix coefficients.
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, '(/ " aptaxis: two or more mix coefficients.")')
cbugcbugc***DEBUG ends.
c.... Find solutions of the characteristic equation. Solve the cubic:
c.... eigen**3 - trace * eigen**2 + dsum * eigen - det = 0.
pp = dsum - trace**2 / 3.0
pperr = dsumerr + 2.0 * tolx * tracerr * abs (trace) / 3.0
if (abs (pp) .le. pperr) pp = 0.0
qq = -2.0 * (trace / 3.0)**3 + trace * dsum / 3.0 - det
qqerr = 2.0 * tracerr * (trace / 3.0)**2 +
& (dsumerr * abs (trace) + tracerr * abs (dsum)) / 3.0 +
& deterr
if (abs (qq) .le. qqerr) qq = 0.0
qd = (pp / 3.0)**3 + (qq / 2.0)**2
qderr = pperr * abs (pp) / 9.0 + qqerr * abs (qq) / 2.0
if (abs (qd) .le. qderr) qd = 0.0
cbugcbugc***DEBUG begins.
cbugcbug 9905 format (/ " Cubic factors: pp, qq, qd=" /
cbugcbug & 2x,1p3e22.14," errors:" / 2x,1p3e22.14)
cbugcbug write ( 3, 9905) pp, qq, qd, pperr, qqerr, qderr
cbugcbugc***DEBUG ends.
if (qd .ge. 0.0) then ! 3 real roots, last 2 same.
fact = sign ((0.5 * abs (qq))**(1.0 / 3.0), qq)
eigen(1) = trace / 3.0 - 2.0 * fact
eigen(2) = trace / 3.0 + fact
eigen(3) = eigen(2)
else ! 3 real roots, all different.
fact = 2.0 * sqrt (-pp / 3.0)
theta = acos (3.0 * qq / (fact * pp))
eigen(1) = trace / 3.0 + fact * cos (theta / 3.0)
eigen(2) = trace / 3.0 - fact * cos ((theta + pi) / 3.0)
eigen(3) = trace / 3.0 - fact * cos ((theta - pi) / 3.0)
endif
cbugcbugc***DEBUG begins.
cbugcbug 9906 format (/ " Eigenvalues of characteristic equation:" /
cbugcbug & " e1,e2,e3=",1p3e22.14)
cbugcbug write ( 3, 9906) eigen(1), eigen(2), eigen(3)
cbugcbugc***DEBUG ends.
c.... Find a principal axis "b" of the surface. Use longest.
blen = 0.0
bx = 1.0
by = 0.0
bz = 0.0
do 120 n = 1, 3 ! Loop over eigenvalues.
call aptvaxb (2.0 * qxx - eigen(n), qxy, qzx,
& qxy, 2.0 * qyy - eigen(n), qyz, 1, 0.0,
& cx, cy, cz, clen, nerr)
if (clen .gt. blen) then
bx = cx
by = cy
bz = cz
blen = clen
endif
call aptvaxb (qxy, 2.0 * qyy - eigen(n), qyz,
& qzx, qyz, 2.0 * qzz - eigen(n), 1, 0.0,
& cx, cy, cz, clen, nerr)
if (clen .gt. blen) then
bx = cx
by = cy
bz = cz
blen = clen
endif
call aptvaxb (qzx, qyz, 2.0 * qzz - eigen(n),
& 2.0 * qxx - eigen(n), qxy, qzx, 1, 0.0,
& cx, cy, cz, clen, nerr)
if (clen .gt. blen) then
bx = cx
by = cy
bz = cz
blen = clen
endif
120 continue ! End of loop over eigenvalues.
c.... Find a unit vector along a principal axis.
call aptvuna (bx, by, bz, 1, 0.0, blen, nerr)
cbugcbugc***DEBUG begins.
cbugcbug
cbugcbug 9907 format (/ " Unit axis vector:" /
cbugcbug & " bx,by,bz=",1p3e22.14)
cbugcbug write ( 3, 9907) bx, by, bz
cbugcbug write ( 3, '(/ " Rotate onto x axis to zero qxy, qzx.")')
cbugcbug write ( 3, 9908) nerr, ntype, qc, qx, qy, qz,
cbugcbug & qxy, qyz, qzx, qxx, qyy, qzz
cbugcbugc***DEBUG ends.
c.... Find rotation matrix to put axis vector "b" on x axis, and
c.... then rotate surface. Makes qxy and qzx zero.
call aptrotv (bx, by, bz, 1., 0., 0., tolx, rotp, nerr)
call aptrois (rotp, 0, 0., 0., 0., qc, qx, qy, qz,
& qxy, qyz, qzx, qxx, qyy, qzz, tolx, nerr)
c.... Find the cumulative rotation operator rotm.
do 130 i = 1, 3
do 125 j = 1, 3
rotms(i,j) = rotm(i,j)
125 continue
130 continue
call aptmprd (3, rotp, rotms, tolx, rotm, nerr)
qxy = 0.0
qzx = 0.0
c.... Test qyz to see if more rotation is needed.
c++++ Rotate qyz to zero.
if (qyz**2 .gt. tolx**2 * qsq) then
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, '(/ " aptaxis: rotating around x axis.")')
cbugcbug write ( 3, 9908) nerr, ntype, qc, qx, qy, qz,
cbugcbug & qxy, qyz, qzx, qxx, qyy, qzz
cbugcbugc***DEBUG ends.
c.... Make qyz zero by rotation around the x axis.
thetax = 0.5 * atan2 (qyz, qzz - qyy)
call aptrota (thetax, 1, 1., 0., 0., tolx, rotp, nerr)
call aptrois (rotp, 0, 0., 0., 0., qc, qx, qy, qz,
& qxy, qyz, qzx, qxx, qyy, qzz, tolx, nerr)
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9908) nerr, ntype, qc, qx, qy, qz,
cbugcbug & qxy, qyz, qzx, qxx, qyy, qzz
cbugcbugc***DEBUG ends.
qyz = 0.0
c.... Find the cumulative rotation operator rotm.
do 140 i = 1, 3
do 135 j = 1, 3
rotms(i,j) = rotm(i,j)
135 continue
140 continue
call aptmprd (3, rotp, rotms, tolx, rotm, nerr)
endif ! Zeroed out qyz.
endif ! Tested number of mix coefficients.
endif ! Non-zero mix coefficients.
c=======================================================================--------
c.... Order qxx, qyy, qxx in equal or decreasing order of absolute value.
c.... This insures that qxx is non-zero.
if ((qzz**2 - qyy**2) .gt. tolx * (qzz**2 + qyy**2)) then
! Swap y and z axes.
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, '(/ " aptaxis: swapping y and z (mag).")')
cbugcbug write ( 3, 9908) nerr, ntype, qc, qx, qy, qz,
cbugcbug & qxy, qyz, qzx, qxx, qyy, qzz
cbugcbugc***DEBUG ends.
call aptrott (0., 1., 0., 0., 0., 1., 0., 0., 1., 0., 1., 0.,
& tolx, rotp, nerr)
call aptrois (rotp, 0, 0., 0., 0., qc, qx, qy, qz,
& qxy, qyz, qzx, qxx, qyy, qzz, tolx, nerr)
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9908) nerr, ntype, qc, qx, qy, qz,
cbugcbug & qxy, qyz, qzx, qxx, qyy, qzz
cbugcbugc***DEBUG ends.
c.... Find the cumulative rotation operator rotm.
do 150 i = 1, 3
do 145 j = 1, 3
rotms(i,j) = rotm(i,j)
145 continue
150 continue
call aptmprd (3, rotp, rotms, tolx, rotm, nerr)
endif ! Swapped y and z axes.
if ((qyy**2 - qxx**2) .gt. tolx * (qyy**2 + qxx**2)) then
! Swap x and y axes.
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, '(/ " aptaxis: swapping x and y (mag).")')
cbugcbug write ( 3, 9908) nerr, ntype, qc, qx, qy, qz,
cbugcbug & qxy, qyz, qzx, qxx, qyy, qzz
cbugcbugc***DEBUG ends.
call aptrott (1., 0., 0., 0., 1., 0., 0., 1., 0., 1., 0., 0.,
& tolx, rotp, nerr)
call aptrois (rotp, 0, 0., 0., 0., qc, qx, qy, qz,
& qxy, qyz, qzx, qxx, qyy, qzz, tolx, nerr)
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9908) nerr, ntype, qc, qx, qy, qz,
cbugcbug & qxy, qyz, qzx, qxx, qyy, qzz
cbugcbugc***DEBUG ends.
c.... Find the cumulative rotation operator rotm.
do 160 i = 1, 3
do 155 j = 1, 3
rotms(i,j) = rotm(i,j)
155 continue
160 continue
call aptmprd (3, rotp, rotms, tolx, rotm, nerr)
endif ! Swapped x and y axes.
if ((qzz**2 - qyy**2) .gt. tolx * (qzz**2 + qyy**2)) then
! Swap y and z axes.
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, '(/ " aptaxis: swapping y and z (mag).")')
cbugcbug write ( 3, 9908) nerr, ntype, qc, qx, qy, qz,
cbugcbug & qxy, qyz, qzx, qxx, qyy, qzz
cbugcbugc***DEBUG ends.
call aptrott (0., 1., 0., 0., 0., 1., 0., 0., 1., 0., 1., 0.,
& tolx, rotp, nerr)
call aptrois (rotp, 0, 0., 0., 0., qc, qx, qy, qz,
& qxy, qyz, qzx, qxx, qyy, qzz, tolx, nerr)
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9908) nerr, ntype, qc, qx, qy, qz,
cbugcbug & qxy, qyz, qzx, qxx, qyy, qzz
cbugcbugc***DEBUG ends.
c.... Find the cumulative rotation operator rotm.
do 170 i = 1, 3
do 165 j = 1, 3
rotms(i,j) = rotm(i,j)
165 continue
170 continue
call aptmprd (3, rotp, rotms, tolx, rotm, nerr)
endif ! Swapped y and z axes.
c=======================================================================--------
c.... If qxx, qyy, qzz are all non-zero, put any with odd sign last.
nswap = 0
if ((qxx * qyy .lt. 0.0) .and. (qyy * qzz .gt. 0.0)) then
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, '(/ " aptaxis: swapping z and x (sign).")')
cbugcbug write ( 3, 9908) nerr, ntype, qc, qx, qy, qz,
cbugcbug & qxy, qyz, qzx, qxx, qyy, qzz
cbugcbugc***DEBUG ends.
c.... Swap z and x axes.
nswap = 1
call aptrott (0., 0., 1., 1., 0., 0., 1., 0., 0., 0., 0., 1.,
& tolx, rotp, nerr)
elseif ((qxx * qyy .lt. 0.0) .and. (qzz * qxx .gt. 0.0)) then
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, '(/ " aptaxis: swapping y and z (sign).")')
cbugcbug write ( 3, 9908) nerr, ntype, qc, qx, qy, qz,
cbugcbug & qxy, qyz, qzx, qxx, qyy, qzz
cbugcbugc***DEBUG ends.
c.... Swap y and z axes.
nswap = 1
call aptrott (0., 1., 0., 0., 0., 1., 0., 0., 1., 0., 1., 0.,
& tolx, rotp, nerr)
endif ! Tested signs of qxx, qyy, qzz, qc.
if (nswap .eq. 1) then ! Complete swap of axes.
call aptrois (rotp, 0, 0., 0., 0., qc, qx, qy, qz,
& qxy, qyz, qzx, qxx, qyy, qzz, tolx, nerr)
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9908) nerr, ntype, qc, qx, qy, qz,
cbugcbug & qxy, qyz, qzx, qxx, qyy, qzz
cbugcbugc***DEBUG ends.
c.... Find the cumulative rotation operator rotm.
do 180 i = 1, 3
do 175 j = 1, 3
rotms(i,j) = rotm(i,j)
175 continue
180 continue
call aptmprd (3, rotp, rotms, tolx, rotm, nerr)
endif ! Completed swap of axes.
c=======================================================================--------
c.... Eliminate qy if equation is linear in y and z.
psq = qx**2 + qy**2 + qz**2
ssq = qxx**2 + qyy**2 + qzz**2
if ((qy**2 .gt. tolx**2 * psq) .and.
& (qz**2 .gt. tolx**2 * psq) .and.
& (qyy**2 .le. tolx**2 * ssq) .and.
& (qzz**2 .le. tolx**2 * ssq)) then
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, '(/ " aptaxis: zero out qy by x axis rotation")')
cbugcbug write ( 3, 9908) nerr, ntype, qc, qx, qy, qz,
cbugcbug & qxy, qyz, qzx, qxx, qyy, qzz
cbugcbugc***DEBUG ends.
c.... Linear in y and z. Make qy zero by rotation around the x axis.
thetax = atan2 (qy, qz)
call aptrota (thetax, 1, 1., 0., 0., tolx, rotp, nerr)
call aptrois (rotp, 0, 0., 0., 0., qc, qx, qy, qz,
& qxy, qyz, qzx, qxx, qyy, qzz, tolx, nerr)
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9908) nerr, ntype, qc, qx, qy, qz,
cbugcbug & qxy, qyz, qzx, qxx, qyy, qzz
cbugcbugc***DEBUG ends.
qy = 0.0
c.... Find the cumulative rotation operator rotm.
do 190 i = 1, 3
do 185 j = 1, 3
rotms(i,j) = rotm(i,j)
185 continue
190 continue
call aptmprd (3, rotp, rotms, tolx, rotm, nerr)
endif ! Linear in y and z.
c=======================================================================--------
c.... Multiply coefficients by -1 if qxx is negative.
if (qxx .lt. 0.0) then ! Make qxx positive.
qxx = -qxx
qyy = -qyy
qzz = -qzz
qx = -qx
qy = -qy
qz = -qz
qc = -qc
endif ! Made qxx positive.
c=======================================================================--------
c.... Swap x and y axes if qxx and qyy are both positive, and qyy > qxx.
if ((qxx .gt. 0.0) .and. (qyy .gt. 0.0) .and. (qyy .gt. qxx)) then
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, '(/ " aptaxis: swapping x and y (mag).")')
cbugcbug write ( 3, 9908) nerr, ntype, qc, qx, qy, qz,
cbugcbug & qxy, qyz, qzx, qxx, qyy, qzz
cbugcbugc***DEBUG ends.
call aptrott (1., 0., 0., 0., 1., 0., 0., 1., 0., 1., 0., 0.,
& tolx, rotp, nerr)
call aptrois (rotp, 0, 0., 0., 0., qc, qx, qy, qz,
& qxy, qyz, qzx, qxx, qyy, qzz, tolx, nerr)
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9908) nerr, ntype, qc, qx, qy, qz,
cbugcbug & qxy, qyz, qzx, qxx, qyy, qzz
cbugcbugc***DEBUG ends.
c.... Find the cumulative rotation operator rotm.
do 200 i = 1, 3
do 195 j = 1, 3
rotms(i,j) = rotm(i,j)
195 continue
200 continue
call aptmprd (3, rotp, rotms, tolx, rotm, nerr)
endif ! Swapped x and y axes.
c=======================================================================--------
c.... Swap y and z axes if a parabolic cylinder is on y axis.
if ((qxx .ne. 0.0) .and. (qyy .eq. 0.0) .and. (qz .ne. 0.0)) then
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, '(/ " aptaxis: swapping y and z (axis).")')
cbugcbug write ( 3, 9902) qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz
cbugcbugc***DEBUG ends.
call aptrott (0., 1., 0., 0., 0., 1., 0., 0., 1., 0., 1., 0.,
& tolx, rotp, nerr)
call aptrois (rotp, 0, 0., 0., 0., qc, qx, qy, qz,
& qxy, qyz, qzx, qxx, qyy, qzz, tolx, nerr)
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9902) qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz
cbugcbugc***DEBUG ends.
c.... Find the cumulative rotation operator rotm.
do 220 i = 1, 3
do 210 j = 1, 3
rotms(i,j) = rotm(i,j)
210 continue
220 continue
call aptmprd (3, rotp, rotms, tolx, rotm, nerr)
endif ! Swapped y and z axes.
c=======================================================================--------
c.... Rotate 180 degrees around x axis if qxx > 0 and qy > 0 or qz > 0.
c.... Orient parabolic cylinders and paraboloids in standard form.
if ((qxx .ne. 0.0) .and. ((qy .gt. 0.0) .or. (qz .gt. 0.0))) then
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, '(/ " aptaxis: rotating 180 deg around x axis.")')
cbugcbug write ( 3, 9902) qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz
cbugcbugc***DEBUG ends.
call aptrota (180., 0, 1., 0., 0., tolx, rotp, nerr)
call aptrois (rotp, 0, 0., 0., 0., qc, qx, qy, qz,
& qxy, qyz, qzx, qxx, qyy, qzz, tolx, nerr)
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9902) qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz
cbugcbugc***DEBUG ends.
c.... Find the cumulative rotation operator rotm.
do 240 i = 1, 3
do 230 j = 1, 3
rotms(i,j) = rotm(i,j)
230 continue
240 continue
call aptmprd (3, rotp, rotms, tolx, rotm, nerr)
endif ! Rotated 180 degrees around x axis.
c=======================================================================--------
c.... Rotate 180 degrees around x axis if qxx, qyy and qzz all non-zero,
c.... and rotm(2,2) = rotm(3,3) = -1 (a rotation with no effect).
if ((qxx .ne. 0.0) .and. (qyy .ne. 0.0) .and. (qzz .ne. 0.0) .and.
& (rotm(2,2) .eq. -1.0) .and. (rotm(3,3) .eq. -1.0)) then
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, '(/ " aptaxis: rotating 180 deg around x axis.")')
cbugcbug write ( 3, 9902) qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz
cbugcbugc***DEBUG ends.
call aptrota (180., 0, 1., 0., 0., tolx, rotp, nerr)
call aptrois (rotp, 0, 0., 0., 0., qc, qx, qy, qz,
& qxy, qyz, qzx, qxx, qyy, qzz, tolx, nerr)
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9902) qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz
cbugcbugc***DEBUG ends.
c.... Find the cumulative rotation operator rotm.
do 260 i = 1, 3
do 250 j = 1, 3
rotms(i,j) = rotm(i,j)
250 continue
260 continue
call aptmprd (3, rotp, rotms, tolx, rotm, nerr)
endif ! Rotated 180 degrees around x axis.
c=======================================================================--------
c.... Make a hyperbolic cylinder intersect the y axis.
c.... Swap x and y axes if qxx > 0, qyy < 0, and qc < 0.
if ((qxx .gt. 0.0) .and. (qyy .lt. 0.0) .and. (qc .lt. 0.0)) then
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, '(/ " aptaxis: swapping x and y (axis).")')
cbugcbug write ( 3, 9902) qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz
cbugcbugc***DEBUG ends.
call aptrott (1., 0., 0., 0., 1., 0., 0., 1., 0., 1., 0., 0.,
& tolx, rotp, nerr)
call aptrois (rotp, 0, 0., 0., 0., qc, qx, qy, qz,
& qxy, qyz, qzx, qxx, qyy, qzz, tolx, nerr)
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9902) qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz
cbugcbugc***DEBUG ends.
c.... Find the cumulative rotation operator rotm.
do 280 i = 1, 3
do 270 j = 1, 3
rotms(i,j) = rotm(i,j)
270 continue
280 continue
call aptmprd (3, rotp, rotms, tolx, rotm, nerr)
qc = -qc
qxx = -qxx
qyy = -qyy
endif ! Swapped x and y axes.
c=======================================================================--------
c.... Find the coordinates of the center of the standard form.
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, '(/ " aptaxis: finding standard center.")')
cbugcbug write ( 3, 9902) qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz
cbugcbugc***DEBUG ends.
ax = 0.0
ay = 0.0
az = 0.0
qcs = qc
qcserr = tolx * qcs
cbugcbugc***DEBUG begins.
cbugcbugc***DEBUG ends.
c++++ Zero out qc.
if (qxx**2 .le. tolx**2 * ssq) then
qxx = 0.0
c++++ Linear in x.
if (qx**2 .gt. tolx**2 * psq) then
ax = - qcs / qx
qcs = 0.0
endif ! Tested qx.
else ! Quadratic in x.
ax = -0.5 * qx / qxx
qcs = qcs + ax * (qx + qxx * ax)
qcserr = qcserr + tolx * abs (ax) * (abs (qx) + abs (qxx * ax))
if (abs (qcs) .le. qcserr) qcs = 0.0
endif ! Tested qxx.
c++++ Zero out qc.
if (qyy**2 .le. tolx**2 * ssq) then
qyy = 0.0
c++++ Linear in y.
if (qy**2 .gt. tolx**2 * psq) then
ay = - qcs / qy
qcs = 0.0
endif ! Tested qy.
else ! Quadratic in y.
ay = -0.5 * qy / qyy
qcs = qcs + ay * (qy + qyy * ay)
qcserr = qcserr + tolx * abs (ay) * (abs (qy) + abs (qyy * ay))
if (abs (qcs) .le. qcserr) qcs = 0.0
endif ! Tested qyy.
c++++ Zero out qc.
if (qzz**2 .le. tolx**2 * ssq) then
qzz = 0.0
c++++ Linear in z.
if (qz**2 .gt. tolx**2 * psq) then
az = - qcs / qz
qcs = 0.0
endif ! Tested qz.
else ! Quadratic in z.
az = -0.5 * qz / qzz
qcs = qcs + az * (qz + qzz * az)
qcserr = qcserr + tolx * abs (az) * (abs (qz) + abs (qzz * az))
if (abs (qcs) .le. qcserr) qcs = 0.0
endif ! Tested qzz.
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9909) ax, ay, az
cbugcbugc***DEBUG ends.
c.... Translate the standard center to the origin.
call apttris (0, ax, ay, az, qc, qx, qy, qz, qxy, qyz, qzx,
& qxx, qyy, qzz, tolx, nerr)
c=======================================================================--------
c.... Find the initial standard center.
call aptmopv (rotm, 1, 0., 0., 0., ax, ay, az, 1, tolx, nerr)
c=======================================================================--------
c.... Find the surface type.
if (qyy .eq. 0.0) then ! Zero qyy, zero qzz.
if ((qz .eq. 0.0) .and.
& (qy .eq. 0.0)) then
! Zero qz.
if (qc .eq. 0.0) then ! Zero qc.
ntype = 1 ! Coincident planes.
elseif (qc .lt. 0.0) then ! Negative qc.
ntype = 2 ! Real parallel planes.
else ! Positive qc.
ntype = 17 ! Imaginary parallel planes.
endif ! Tested qc.
else ! Non-zero qz.
ntype = 4 ! Parabolic cylinder.
endif ! Tested qz.
elseif (qyy .lt. 0.0) then ! Negative qyy, zero qzz.
if (qz .eq. 0.0) then ! Zero qz.
if (qc .eq. 0.0) then ! Zero qc.
ntype = 3 ! Real intersecting planes.
else ! Non-zero qc.
ntype = 5 ! Hyperbolic cylinder.
endif ! Tested qc.
else ! Non-zero qz.
ntype = 8 ! Hyperbolic paraboloid.
endif ! Tested qz.
else ! Positive qyy.
if (qzz .eq. 0.0) then ! Positive qyy, zero qzz.
if (qz .eq. 0.0) then ! Zero qz.
if (qc .lt. 0.0) then ! Negative qc.
if (abs (qyy - qxx) .gt. tolx * qxx) then
ntype = 6 ! Real elliptic cylinder.
else
ntype = 7 ! Real circular cylinder.
endif
elseif (qc .eq. 0.0) then ! Zero qc.
ntype = 18 ! Imaginary intersecting planes.
else ! Positive qc.
if (abs (qyy - qxx) .gt. tolx * qxx) then
ntype = 19 ! Imaginary elliptic cylinder.
else
ntype = 20 ! Imaginary circular cylinder.
endif
endif ! Tested qc.
else ! Non-zero qz.
if (abs (qyy - qxx) .gt. tolx * qxx) then
ntype = 9 ! Elliptic paraboloid.
else
ntype = 10 ! Circular paraboloid.
endif
endif ! Tested qz.
elseif (qzz .lt. 0.0) then ! Positive qyy, negative qzz.
if (qc .eq. 0.0) then ! Zero qc.
if (abs (qyy - qxx) .gt. tolx * qxx) then
ntype = 11 ! Real elliptic cone.
else
ntype = 12 ! Real circular cone.
endif
elseif (qc .lt. 0.0) then ! Negative qc.
ntype = 13 ! Hyperboloid of one sheet.
else ! Positive qc.
ntype = 14 ! Hyperboloid of two sheets.
endif ! Tested qc.
else ! Positive qyy, positive qzz.
if (qc .lt. 0.0) then ! Negative qc.
if ((abs (qyy - qxx) .gt. tolx * qxx) .or.
& (abs (qzz - qyy) .gt. tolx * qyy) .or.
& (abs (qxx - qzz) .gt. tolx * qzz)) then
ntype = 15 ! Real ellipsoid.
else
ntype = 16 ! Real sphere.
endif
elseif (qc .eq. 0.0) then ! Zero qc.
if ((abs (qyy - qxx) .gt. tolx * qxx) .and.
& (abs (qzz - qyy) .gt. tolx * qyy) .and.
& (abs (qxx - qzz) .gt. tolx * qzz)) then
ntype = 21 ! Imaginary elliptic cone.
else
ntype = 22 ! Imaginary circular cone.
endif
else ! Positive qc.
if (abs (qyy - qxx) .gt. tolx * qxx) then
ntype = 23 ! Imaginary ellipsoid.
else
ntype = 24 ! Imaginary sphere.
endif
endif ! Tested qc.
endif ! Tested qzz.
endif ! Tested qyy.
c=======================================================================--------
310 continue
cbugc***DEBUG begins.
cbug 9908 format (/ 'aptaxis results: nerr=',i3,' ntype=',i3 /
cbug & ' qc= ',1pe22.14 /
cbug & ' qx,qy,qz= ',1p3e22.14 /
cbug & ' qxy,qyz,qzx=',1p3e22.14 /
cbug & ' qxx,qyy,qzz=',1p3e22.14)
cbug 9909 format (/ ' ax,ay,az= ',1p3e22.14)
cbug 9910 format (/ ' rotm=',3(/ 14x,1p3e22.14))
cbug write ( 3, 9908) nerr, ntype, qc, qx, qy, qz,
cbug & qxy, qyz, qzx, qxx, qyy, qzz
cbug write ( 3, 9909) ax, ay, az
cbug write ( 3, 9910) ((rotm(i,j), j = 1, 3), i = 1, 3)
cbugc.... Recalculate the invariants of the surface.
cbug traces = trace
cbug dsums = dsum
cbug dets = det
cbug trace = 2.0 * (qxx + qyy + qzz)
cbug dsum = 4.0 * (qxx * qyy + qyy * qzz + qzz * qxx) -
cbug & (qxy**2 + qyz**2 + qzx**2)
cbug det = 2.0 * (4.0 * qxx * qyy * qzz + qxy * qyz * qzx -
cbug & (qxx * qyz**2 + qyy * qzx**2 + qzz * qxy**2))
cbug tracerr = 2.0 * tolx * (abs (qxx) + abs (qyy) + abs (qzz))
cbug dsumerr = 4.0 * tolx * (abs (qxx * qyy) +
cbug & abs (qyy * qzz) +
cbug & abs (qzz * qxx)) +
cbug & tolx * (qxy**2 + qyz**2 + qzx**2)
cbug deterr = 2.0 * tolx *
cbug & (4.0 * abs (qxx * qyy * qzz) +
cbug & abs (qxy * qyz * qzx) +
cbug & qyz**2 * abs (qxx) +
cbug & qzx**2 * abs (qyy) +
cbug & qxy**2 * abs (qzz))
cbug
cbug if (abs (trace) .le. tracerr) trace = 0.0
cbug if (abs (dsum) .le. dsumerr) dsum = 0.0
cbug if (abs (det) .le. deterr) det = 0.0
cbug tracerr = abs (trace) - abs (traces)
cbug dsumerr = dsum - dsums
cbug deterr = abs (det) - abs (dets)
cbug write ( 3, 9904) trace, dsum, det, tracerr, dsumerr, deterr
cbugc***DEBUG ends.
return
c.... End of subroutine aptaxis. (+1 line.)
end
UCRL-WEB-209832