subroutine aptqper (ac, ax, ay, az, axy, ayz, azx, axx, ayy, azz,
& bc, bx, by, bz, bxy, byz, bzx, bxx, byy, bzz,
& tol,
& cc, cx, cy, cz, cxy, cyz, czx, cxx, cyy, czz,
& nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTQPER
c
c call aptqper (ac, ax, ay, az, axy, ayz, azx, axx, ayy, azz,
c & bc, bx, by, bz, bxy, byz, bzx, bxx, byy, bzz, tol,
c & cc, cx, cy, cz, cxy, cyz, czx, cxx, cyy, czz, nerr)
c
c Version: aptqper Updated 2000 May 24 14:00.
c aptqper Originated 2000 May 24 14:00.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: For the two families of quadric surfaces represented by
c quadric surfaces A(x,y,z) and B(x,y,z), with arbitrary values
c of the constant terms ac and bc:
c
c A = ac + ax * x + ay * y + az * z
c + axy * x * y + ayz * y * z + azx * z * x
c + axx * x**2 + ayy * y**2 + azz * z**2 = 0
c B = bc + bx * x + by * y + bz * z
c + bxy * x * y + byz * y * z + bzx * z * x
c + bxx * x**2 + byy * y**2 + bzz * z**2 = 0
c
c to find the quadric surface C(x,y,z) representing the locus
c of points where the families A and B are perpendicular to
c each other:
c
c C = cc + cx * x + cy * y + cz * z
c + cxy * x * y + cyz * y * z + czx * z * x
c + cxx * x**2 + cyy * y**2 + czz * z**2 = 0
c
c Method: the two quadric surfaces have the normal vectors:
c
c AN = (anx, any, anz) = grad (A)
c anx = ax + 2 * axx * x + axy * y + azx * z
c any = ay + axy * x + 2 * ayy * y + ayz * z
c anz = az + azx * x + ayz * y + 2 * azz * z
c
c BN = (bnx, bny, bnz) = grad (B)
c bnx = bx + 2 * bxx * x + bxy * y + bzx * z
c bny = by + bxy * x + 2 * byy * y + byz * z
c bnz = bz + bzx * x + byz * y + 2 * bzz * z
c
c Then C(x,y,z) = AN dot BN = 0, where the normal vectors are
c perpendicular.
c
c Input: ac, ax, ay, az, axy, ayz, azx, axx, ayy, azz,
c bc, bx, by, bz, bxy, byz, bzx, bxx, byy, bzz, tol.
c
c Output: cc, cx, cy, cz, cxy, cyz, czx, cxx, cyy, czz, nerr.
c
c Calls: aptaxis (source, binary in aptaxis (source, binary in Compass
c Cluster East machine, ~edwards/work/apt/src)
c
c Glossary:
c
c ac, ... Input Coefficients of the implicit equation for
c quadric surface A:
c ac + ax * x + ay * y + az * z +
c axy * x * y + ayz * y * z + azx * z * x +
c axx * x**2 + ayy * y**2 + azz * z**2 = 0
c
c bc, ... Input Coefficients of the implicit equation for
c quadric surface B:
c bc + bx * x + by * y + bz * z +
c bxy * x * y + byz * y * z + bzx * z * x +
c bxx * x**2 + byy * y**2 + bzz * z**2 = 0
c
c cc, ... Output Coefficients of the implicit equation for
c quadric surface C:
c cc + cx * x + cy * y + cz * z +
c cxy * x * y + cyz * y * z + czx * z * x +
c cxx * x**2 + cyy * y**2 + czz * z**2 = 0
c
c nerr Output Indicates an input error, if not zero:
c 1 if the equation of quadric surface C is bad.
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 variables.
dimension rotm (3,3) ! Coordinate transformation matrix.
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptqper finding quadric surface where two',
cbug & ' quadric surface families' /
cbug & ' are perpendicular.',' tol=',1pe13.5)
cbug 9903 format ('Quadric ac =',1pe22.14 /
cbug & ' ax,ay,az =',1p3e22.14 /
cbug & ' axy,ayz,azx=',1p3e22.14 /
cbug & ' axx,ayy,azz=',1p3e22.14 )
cbug 9904 format ('Quadric bc =',1pe22.14 /
cbug & ' bx,by,bz =',1p3e22.14 /
cbug & ' bxy,byz,bzx=',1p3e22.14 /
cbug & ' bxx,byy,bzz=',1p3e22.14 )
cbug write ( 3, 9901) tolx
cbug write ( 3, 9903) ac, ax, ay, az, axy, ayz, azx, axx, ayy, azz
cbug write ( 3, 9904) bc, bx, by, bz, bxy, byz, bzx, bxx, byy, bzz
cbugc***DEBUG ends.
c.... Initialize.
nerr = 0
c.... Find the coefficients of quadric surface C.
cc = ax * bx + ay * by + az * bz
cx = 2.0 * axx * bx + 2.0 * ax * bxx +
& ay * bxy + axy * by + az * bzx + azx * bz
cy = 2.0 * ayy * by + 2.0 * ay * byy +
& az * byz + ayz * bz + ax * bxy + axy * bx
cz = 2.0 * azz * bz + 2.0 * az * bzz +
& ax * bzx + azx * bx + ay * byz + ayz * by
cxy = 2.0 * axx * bxy + 2.0 * axy * bxx +
& 2.0 * ayy * bxy + 2.0 * axy * byy +
& azx * byz + ayz * bzx
cyz = 2.0 * ayy * byz + 2.0 * ayz * byy +
& 2.0 * azz * byz + 2.0 * ayz * bzz +
& axy * bzx + azx * bxy
czx = 2.0 * azz * bzx + 2.0 * azx * bzz +
& 2.0 * axx * bzx + 2.0 * azx * bxx +
& ayz * bxy + axy * byz
cxx = 4.0 * axx * bxx + axy * bxy + azx * bzx
cyy = 4.0 * ayy * byy + ayz * byz + axy * bxy
czz = 4.0 * azz * bzz + azx * bzx + ayz * byz
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9992) cc, cx, cy, cz, cxy, cyz, czx, cxx, cyy, czz
cbugcbugc***DEBUG ends.
c.... Replace coefficients with zero, if within truncation limits.
if (tol .gt. 0.0) then ! Truncate small coefficients to zero.
ccerr = abs (ax * bx) + abs (ay * by) + abs (az * bz)
cxerr = 2.0 * abs (axx * bx) + 2.0 * abs (ax * bxx) +
& abs (ay * bxy) + abs (axy * by) +
& abs (az * bzx) + abs (azx * bz)
cyerr = 2.0 * abs (ayy * by) + 2.0 * abs (ay * byy) +
& abs (az * byz) + abs (ayz * bz) +
& abs (ax * bxy) + abs (axy * bx)
czerr = 2.0 * abs (azz * bz) + 2.0 * abs (az * bzz) +
& abs (ax * bzx) + abs (azx * bx) +
& abs (ay * byz) + abs (ayz * by)
cxyerr = 2.0 * abs (axx * bxy) + 2.0 * abs (axy * bxx) +
& 2.0 * abs (ayy * bxy) + 2.0 * abs (axy * byy) +
& abs (azx * byz) + abs (ayz * bzx)
cyzerr = 2.0 * abs (ayy * byz) + 2.0 * abs (ayz * byy) +
& 2.0 * abs (azz * byz) + 2.0 * abs (ayz * bzz) +
& abs (axy * bzx) + abs (azx * bxy)
czxerr = 2.0 * abs (azz * bzx) + 2.0 * abs (azx * bzz) +
& 2.0 * abs (axx * bzx) + 2.0 * abs (azx * bxx) +
& abs (ayz * bxy) + abs (axy * byz)
cxxerr = 4.0 * abs (axx * bxx) + abs (axy * bxy) + abs (azx * bzx)
cyyerr = 4.0 * abs (ayy * byy) + abs (ayz * byz) + abs (axy * bxy)
czzerr = 4.0 * abs (azz * bzz) + abs (azx * bzx) + abs (ayz * byz)
if (abs (cc) .le. tol * ccerr ) cc = 0.0
if (abs (cx) .le. tol * cxerr ) cx = 0.0
if (abs (cy) .le. tol * cyerr ) cy = 0.0
if (abs (cz) .le. tol * czerr ) cz = 0.0
if (abs (cxy) .le. tol * cxyerr) cxy = 0.0
if (abs (cyz) .le. tol * cyzerr) cyz = 0.0
if (abs (czx) .le. tol * czxerr) czx = 0.0
if (abs (cxx) .le. tol * cxxerr) cxx = 0.0
if (abs (cyy) .le. tol * cyyerr) cyy = 0.0
if (abs (czz) .le. tol * czzerr) czz = 0.0
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9992) cc, cx, cy, cz, cxy, cyz, czx, cxx, cyy, czz
cbugcbugc***DEBUG ends.
endif ! Truncated small coefficients to zero.
c=======================================================================********
c.... See if quadric surface C is real, not imaginary or impossible.
dc = cc
dx = cx
dy = cy
dz = cz
dxy = cxy
dyz = cyz
dzx = czx
dxx = cxx
dyy = cyy
dzz = czz
call aptaxis (dc, dx, dy, dz, dxy, dyz, dzx, dxx, dyy, dzz,
& tol, cenx, ceny, cenz, rotm, ntype, nerra)
cbugcbugc***DEBUG begins.
cbugcbug 9998 format ('DEBUG ntype=',i3,' nerra=',i3)
cbugcbug write ( 3, 9998) ntype, nerra
cbugcbugc***DEBUG ends.
if ((nerra .ne. 0) .or.
& (ntype .lt. 0) .or.
& (ntype .gt. 16)) then
nerr = 1 ! Imaginary or impossible surface.
endif
c=======================================================================********
999 continue
cbugc***DEBUG begins.
cbug 9991 format (/ 'aptqper results. nerr=',i2)
cbug 9992 format ('Quadric cc =',1pe22.14 /
cbug & ' cx,cy,cz =',1p3e22.14 /
cbug & ' cxy,cyz,czx=',1p3e22.14 /
cbug & ' cxx,cyy,czz=',1p3e22.14 )
cbug write ( 3, 9991) nerr
cbug write ( 3, 9992) cc, cx, cy, cz, cxy, cyz, czx, cxx, cyy, czz
cbugc***DEBUG ends.
return
c.... End of subroutine aptqper. (+1 line)
end
UCRL-WEB-209832