subroutine aptqtan (ax, ay, az, anx, any, anz, bx, by, bz,
& bnx, bny, bnz, cx, cy, cz, cnx, cny, cnz, tol,
& qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz,
& nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTQTAN
c
c call aptqtan (ax, ay, az, anx, any, anz, bx, by, bz,
c & bnx, bny, bnz, cx, cy, cz, cnx, cny, cnz, tol,
c & qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz, nerr)
c
c Version: aptqtan Updated 1999 March 18 12:00.
c aptqtan Originated 1999 February 23 17:45.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To find the coefficients of the implicit equation of any unique
c quadric surface tangent to three planes at
c point a = (ax, ay, az) with normal vector an = (anx, any, anz),
c point b = (bx, by, bz) with normal vector bn = (bnx, bny, bnz),
c point c = (cx, cy, cz) with normal vector cn = (cnx, cny, cnz).
c The equation of the quadric surface is:
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 normal vector is n(x, y, z) = (nx, ny,nz), where
c
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 If any points coincide, or a quadric surface fit is impossible
c or not unique, the method fails, and nerr = 1 is returned.
c
c Input: ax, ay, az, anx, any, anz, bx, by, bz, bnx, bny, bnz,
c cx, cy, cz, cnx, cny, cnz, tol.
c
c Output: qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz, nerr.
c
c Calls: aptsolv, apttris
c
c
c Glossary:
c
c ax,ay,az Input The x, y and z coordinates of point "a".
c
c anx,y,z Input The x, y and z components of the normal vector "an"
c at point "a".
c
c bx,by,bz Input The x, y and z coordinates of point "b".
c
c bnx,y,z Input The x, y and z components of the normal vector "bn"
c at point "b".
c
c cx,cy,cz Input The x, y end z coordinates of point "c".
c
c cnx,y,z Input The x, y anf z components of the normal vector "cn"
c at point "c".
c
c nerr Output Error indicator. 0 if a good solution was obtained.
c 1 if the matrix of coefficients "ele" is singular.
c The determinant of the matrix "ele" is zero.
c No unique solution exists.
c 2 if any residual "dr" exceeds 1.e-6 times the
c maximum d value. Results may be usable.
c 3 if any residual "solr" exceeds 1.e-6 times the
c maximum sol value. Results may be usable.
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
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.
c
c tol Input Numerical tolerance limit.
c On computers with 48-bit floating-point numbers,
c recommend 1.e-5 to 1.e-11.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Local variables.
common /lhyperb/ ele(12,12) ! Coefficients of equations (row, column).
common /lhyperb/ ein(12,12) ! Inverse of matrix ele (row, column).
common /lhyperb/ d(12) ! Right-hand side constants.
common /lhyperb/ dm(12) ! Right-hand side constants check.
common /lhyperb/ dr(12) ! Right-hand side constants residuals.
common /lhyperb/ work(12,13) ! Working matrix.
common /lhyperb/ sol(12) ! Solution.
common /lhyperb/ solm(12) ! Solution check.
common /lhyperb/ solr(12) ! Solution residuals.
c***DEBUG begins.
9901 format (/ 'aptqtan finding tangent quadric. tol=',1pe13.5)
9902 format (/ 'ax, ay, az= ',1p3e22.14 /
& 'anx,any,anz= ',1p3e22.14 /
& 'bx, by, bz= ',1p3e22.14 /
& 'bnx,bny,bnz= ',1p3e22.14 /
& 'cx, cy, cz= ',1p3e22.14 /
& 'cnx,cny,cnz= ',1p3e22.14 )
write ( 3, 9901) tol
write ( 3, 9902) ax, ay, az, anx, any, anz,
& bx, by, bz, bnx, bny, bnz,
& cx, cy, cz, cnx, cny, cnz
c***DEBUG ends.
c.... Initialize.
nerr = 0
qc = -1.e99
qx = -1.e99
qy = -1.e99
qz = -1.e99
qxy = -1.e99
qyz = -1.e99
qzx = -1.e99
qxx = -1.e99
qyy = -1.e99
qzz = -1.e99
bmult = -1.e99
cmult = -1.e99
c.... See if the normal vectors have nonzero magnitudes.
call aptvusz (anx, any, anz, tol, aux, auy, auz, vlena)
call aptvusz (bnx, bny, bnz, tol, bux, buy, buz, vlenb)
call aptvusz (cnx, cny, cnz, tol, cux, cuy, cuz, vlenc)
if ((vlena .eq. 0.0) .or.
& (vlenb .eq. 0.0) .or.
& (vlenc .eq. 0.0)) then
nerr = 1
go to 210
endif
c.... Terms 1-12 of each equation (second index) are coefficients of
c.... qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy and qzz, bmult and cmult,
c.... where na = (anx, any, anz) * 1.0,
c.... nb = -(bnx, bny, bnz) * bmult,
c.... nc = -(cnx, cny, cnz) * cmult.
c-----------------------------------------------------------------------********
c.... Find the coefficients for point "a", equation 1.
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.
ele(1,1) = 1.0 ! Coefficient of qc.
ele(1,2) = ax ! Coefficient of qx.
ele(1,3) = ay ! Coefficient of qy.
ele(1,4) = az ! Coefficient of qz.
ele(1,5) = ax * ay ! Coefficient of qxy.
ele(1,6) = ay * az ! Coefficient of qyz.
ele(1,7) = az * ax ! Coefficient of qzx.
ele(1,8) = ax**2 ! Coefficient of qxx.
ele(1,9) = ay**2 ! Coefficient of qyy.
ele(1,10) = az**2 ! Coefficient of qzz.
ele(1,11) = 0.0 ! Coefficient of bmult.
ele(1,12) = 0.0 ! Coefficient of cmult.
d(1) = 0.0 ! Right side of equation.
c.... Find the coefficients for normal vector component "anx", equation 2.
c.... nx = qx + 2 * qxx * x + qxy * y + qzx * z = anx.
ele(2,1) = 0.0 ! Coefficient of qc.
ele(2,2) = 1.0 ! Coefficient of qx.
ele(2,3) = 0.0 ! Coefficient of qy.
ele(2,4) = 0.0 ! Coefficient of qz.
ele(2,5) = ay ! Coefficient of qxy.
ele(2,6) = 0.0 ! Coefficient of qyz.
ele(2,7) = az ! Coefficient of qzx.
ele(2,8) = 2.0 * ax ! Coefficient of qxx.
ele(2,9) = 0.0 ! Coefficient of qyy.
ele(2,10) = 0.0 ! Coefficient of qzz.
ele(2,11) = 0.0 ! Coefficient of bmult.
ele(2,12) = 0.0 ! Coefficient of cmult.
d(2) = anx ! Right side of equation.
c.... Find the coefficients for normal vector component "any", equation 3.
c.... ny = qy + qxy * x + 2 * qyy * y + qyz * z = any.
ele(3,1) = 0.0 ! Coefficient of qc.
ele(3,2) = 0.0 ! Coefficient of qx.
ele(3,3) = 1.0 ! Coefficient of qy.
ele(3,4) = 0.0 ! Coefficient of qz.
ele(3,5) = ax ! Coefficient of qxy.
ele(3,6) = az ! Coefficient of qyz.
ele(3,7) = 0.0 ! Coefficient of qzx.
ele(3,8) = 0.0 ! Coefficient of qxx.
ele(3,9) = 2.0 * ay ! Coefficient of qyy.
ele(3,10) = 0.0 ! Coefficient of qzz.
ele(3,11) = 0.0 ! Coefficient of bmult.
ele(3,12) = 0.0 ! Coefficient of cmult.
d(3) = any ! Right side of equation.
c.... Find the coefficients for normal vector component "anz", equation 4.
c.... nz = qz + qzx * x + qyz * y + 2 * qzz * z = anz.
ele(4,1) = 0.0 ! Coefficient of qc.
ele(4,2) = 0.0 ! Coefficient of qx.
ele(4,3) = 0.0 ! Coefficient of qy.
ele(4,4) = 1.0 ! Coefficient of qz.
ele(4,5) = 0.0 ! Coefficient of qxy.
ele(4,6) = ay ! Coefficient of qyz.
ele(4,7) = ax ! Coefficient of qzx.
ele(4,8) = 0.0 ! Coefficient of qxx.
ele(4,9) = 0.0 ! Coefficient of qyy.
ele(4,10) = 2.0 * az ! Coefficient of qzz.
ele(4,11) = 0.0 ! Coefficient of bmult.
ele(4,12) = 0.0 ! Coefficient of cmult.
d(4) = anz ! Right side of equation.
c-----------------------------------------------------------------------********
c.... Find the coefficients for point "b", equation 5.
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.
ele(5,1) = 1.0 ! Coefficient of qc.
ele(5,2) = bx ! Coefficient of qx.
ele(5,3) = by ! Coefficient of qy.
ele(5,4) = bz ! Coefficient of qz.
ele(5,5) = bx * by ! Coefficient of qxy.
ele(5,6) = by * bz ! Coefficient of qyz.
ele(5,7) = bz * bx ! Coefficient of qzx.
ele(5,8) = bx**2 ! Coefficient of qxx.
ele(5,9) = by**2 ! Coefficient of qyy.
ele(5,10) = bz**2 ! Coefficient of qzz.
ele(5,11) = 0.0 ! Coefficient of bmult.
ele(5,12) = 0.0 ! Coefficient of cmult.
d(5) = 0.0 ! Right side of equation.
c.... Find the coefficients for normal vector component "bnx", equation 6.
c.... nx = qx + 2 * qxx * x + qxy * y + qzx * z + bmult * bnx = 0.0
ele(6,1) = 0.0 ! Coefficient of qc.
ele(6,2) = 1.0 ! Coefficient of qx.
ele(6,3) = 0.0 ! Coefficient of qy.
ele(6,4) = 0.0 ! Coefficient of qz.
ele(6,5) = by ! Coefficient of qxy.
ele(6,6) = 0.0 ! Coefficient of qyz.
ele(6,7) = bz ! Coefficient of qzx.
ele(6,8) = 2.0 * bx ! Coefficient of qxx.
ele(6,9) = 0.0 ! Coefficient of qyy.
ele(6,10) = 0.0 ! Coefficient of qzz.
ele(6,11) = bnx ! Coefficient of bmult.
ele(6,12) = 0.0 ! Coefficient of cmult.
d(6) = 0.0 ! Right side of equation.
c.... Find the coefficients for normal vector component "bny", equation 7.
c.... ny = qy + qxy * x + 2 * qyy * y + qyz * z + bmult * bny = 0.0
ele(7,1) = 0.0 ! Coefficient of qc.
ele(7,2) = 0.0 ! Coefficient of qx.
ele(7,3) = 1.0 ! Coefficient of qy.
ele(7,4) = 0.0 ! Coefficient of qz.
ele(7,5) = bx ! Coefficient of qxy.
ele(7,6) = bz ! Coefficient of qyz.
ele(7,7) = 0.0 ! Coefficient of qzx.
ele(7,8) = 0.0 ! Coefficient of qxx.
ele(7,9) = 2.0 * by ! Coefficient of qyy.
ele(7,10) = 0.0 ! Coefficient of qzz.
ele(7,11) = bny ! Coefficient of bmult.
ele(7,12) = 0.0 ! Coefficient of cmult.
d(7) = 0.0 ! Right side of equation.
c.... Find the coefficients for normal vector component "bnz", equation 8.
c.... ny = qy + qxy * x + 2 * qyy * y + qyz * z + bmult * bnz = 0.0
ele(8,1) = 0.0 ! Coefficient of qc.
ele(8,2) = 0.0 ! Coefficient of qx.
ele(8,3) = 0.0 ! Coefficient of qy.
ele(8,4) = 1.0 ! Coefficient of qz.
ele(8,5) = 0.0 ! Coefficient of qxy.
ele(8,6) = by ! Coefficient of qyz.
ele(8,7) = bx ! Coefficient of qzx.
ele(8,8) = 0.0 ! Coefficient of qxx.
ele(8,9) = 0.0 ! Coefficient of qyy.
ele(8,10) = 2.0 * bz ! Coefficient of qzz.
ele(8,11) = bnz ! Coefficient of bmult.
ele(8,12) = 0.0 ! Coefficient of cmult.
d(8) = 0.0 ! Right side of equation.
c-----------------------------------------------------------------------********
c.... Find the coefficients for point "c", equation 9.
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.
ele(9,1) = 1.0 ! Coefficient of qc.
ele(9,2) = cx ! Coefficient of qx.
ele(9,3) = cy ! Coefficient of qy.
ele(9,4) = cz ! Coefficient of qz.
ele(9,5) = cx * cy ! Coefficient of qxy.
ele(9,6) = cy * cz ! Coefficient of qyz.
ele(9,7) = cz * cx ! Coefficient of qzx.
ele(9,8) = cx**2 ! Coefficient of qxx.
ele(9,9) = cy**2 ! Coefficient of qyy.
ele(9,10) = cz**2 ! Coefficient of qzz.
ele(9,11) = 0.0 ! Coefficient of bmult.
ele(9,12) = 0.0 ! Coefficient of cmult.
d(9) = 0.0 ! Right side of equation.
c.... Find the coefficients for normal vector component "cnx", equation 10.
c.... nx = qx + 2 * qxx * x + qxy * y + qzx * z + cmult * cnx = 0.0
ele(10,1) = 0.0 ! Coefficient of qc.
ele(10,2) = 1.0 ! Coefficient of qx.
ele(10,3) = 0.0 ! Coefficient of qy.
ele(10,4) = 0.0 ! Coefficient of qz.
ele(10,5) = cy ! Coefficient of qxy.
ele(10,6) = 0.0 ! Coefficient of qyz.
ele(10,7) = cz ! Coefficient of qzx.
ele(10,8) = 2.0 * cx ! Coefficient of qxx.
ele(10,9) = 0.0 ! Coefficient of qyy.
ele(10,10) = 0.0 ! Coefficient of qzz.
ele(10,11) = 0.0 ! Coefficient of bmult.
ele(10,12) = cnx ! Coefficient of cmult.
d(10) = 0.0 ! Right side of equation.
c.... Find the coefficients for normal vector component "cny", equation 11.
c.... ny = qy + qxy * x + 2 * qyy * y + qyz * z + cmult * cny = 0.0
ele(11,1) = 0.0 ! Coefficient of qc.
ele(11,2) = 0.0 ! Coefficient of qx.
ele(11,3) = 1.0 ! Coefficient of qy.
ele(11,4) = 0.0 ! Coefficient of qz.
ele(11,5) = cx ! Coefficient of qxy.
ele(11,6) = cz ! Coefficient of qyz.
ele(11,7) = 0.0 ! Coefficient of qzx.
ele(11,8) = 0.0 ! Coefficient of qxx.
ele(11,9) = 2.0 * cy ! Coefficient of qyy.
ele(11,10) = 0.0 ! Coefficient of qzz.
ele(11,11) = 0.0 ! Coefficient of bmult.
ele(11,12) = cny ! Coefficient of cmult.
d(11) = 0.0 ! Right side of equation.
c.... Find the coefficients for normal vector component "cnz", equation 12.
c.... ny = qy + qxy * x + 2 * qyy * y + qyz * z + cmult * cnz = 0.0
ele(12,1) = 0.0 ! Coefficient of qc.
ele(12,2) = 0.0 ! Coefficient of qx.
ele(12,3) = 0.0 ! Coefficient of qy.
ele(12,4) = 1.0 ! Coefficient of qz.
ele(12,5) = 0.0 ! Coefficient of qxy.
ele(12,6) = cy ! Coefficient of qyz.
ele(12,7) = cx ! Coefficient of qzx.
ele(12,8) = 0.0 ! Coefficient of qxx.
ele(12,9) = 0.0 ! Coefficient of qyy.
ele(12,10) = 2.0 * cz ! Coefficient of qzz.
ele(12,11) = 0.0 ! Coefficient of bmult.
ele(12,12) = cnz ! Coefficient of cmult.
d(12) = 0.0 ! Right side of equation.
c-----------------------------------------------------------------------********
c***DEBUG begins.
9820 format (12f6.2)
write ( 3, 9820) ((ele(mm,nn), nn = 1, 12), mm = 1, 12)
write ( 3, 9820) (d(mm), mm = 1, 12)
c***DEBUG ends.
c.... Solve the 12 simultaneous equations.
call aptsolv (ele, d, 12, work, 12, tol, sol, ein, solm, solr,
& dm, dr, nerr)
if (nerr .eq. 1) then
go to 210
endif
c.... Store data for the quadric surface.
qc = sol(1)
qx = sol(2)
qy = sol(3)
qz = sol(4)
qxy = sol(5)
qyz = sol(6)
qzx = sol(7)
qxx = sol(8)
qyy = sol(9)
qzz = sol(10)
bmult = sol(11)
cmult = sol(12)
210 continue
c***DEBUG begins.
9910 format (/ 'aptqtan results. nerr = ',i2)
9911 format (/ 'nerr =',i2,'. 1 = singular matrix.',
& ' 2 = large q residuals. 3 = large s residuals.')
9912 format (/ 'unknowns qc =',1pe22.14 /
& 'qx, qy, qz =',1p3e22.14 /
& 'qxy, qyz, qzx=',1p3e22.14 /
& 'qxx, qyy, qzz=',1p3e22.14 )
9913 format (/ 'bmult, cmult =',1p2e22.14 )
write ( 3, 9910) nerr
if (nerr .ne. 0) then
write ( 3, 9911) nerr
endif
write ( 3, 9912) qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz
write ( 3, 9913) bmult, cmult
c***DEBUG ends.
return
c.... End of subroutine aptqtan. (+1 line.)
end
UCRL-WEB-209832