subroutine aptqfit (ax, ay, az, bx, by, bz, cx, cy, cz,
& dx, dy, dz, ex, ey, ez, fx, fy, fz,
& gx, gy, gz, hx, hy, hz, px, py, pz, tol,
& qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz,
& nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTQFIT
c
c call aptqfit (ax, ay, az, bx, by, bz, cx, cy, cz,
c & dx, dy, dz, ex, ey, ez, fx, fy, fz,
c & gx, gy, gz, hx, hy, hz, px, py, pz, tol,
c & qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz, nerr)
c
c Version: aptqfit Updated 1999 February 23 17:45.
c aptqfit 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 any unique quadric surface through the nine points
c a = (ax, ay, az), b = (bx, by, bz), c = (cx, cy, cz),
c d = (dx, dy, dz), e = (ex, ey, ez), f = (fx, fy, fz),
c g = (gx, gy, gz), h = (hx, hy, hz), p = (px, py, pz).
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 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, bx, by, bz, cx, cy, cz, dx, dy, dz,
c ex, ey, ez, fx, fy, fz, gx, gy, gz, hx, hy, hz,
c px, py, pz, 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 bx,by,bz Input The x, y and z coordinates of point "b".
c
c cx,cy,cz Input The x, y and z coordinates of point "c".
c
c dx,dy,dz Input The x, y and z coordinates of point "d".
c
c ex,ey,ez Input The x, y end z coordinates of point "e".
c
c fx,fy,fz Input The x, y and z coordinates of point "f".
c
c gx,gy,gz Input The x, y and z goordinates of point "g".
c
c hx,hy,hz Input The x, y anh z coordinates of point "h".
c
c px,py,pz Input The x, y anp z coordinates of point "p".
c
c nerr Output Error indicator. 0 if a good solution was obtained.
c 1 if the matrix of coefficients "abc" is singular.
c The determinant of the matrix "abc" 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 "xyzr" exceeds 1.e-6 times the
c maximum xyz 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 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.
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 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/ abc(9,9) ! Coefficients of equations (row, column).
common /lhyperb/ cba(9,9) ! Inverse of matrix abc (row, column).
common /lhyperb/ d(9) ! Right-hand side constants.
common /lhyperb/ dm(9) ! Right-hand side constants check.
common /lhyperb/ dr(9) ! Right-hand side constants residuals.
common /lhyperb/ work(9,10) ! Working matrix.
common /lhyperb/ xyz(9) ! Solution.
common /lhyperb/ xyzm(9) ! Solution check.
common /lhyperb/ xyzr(9) ! Solution residuals.
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptqfit finding quadric surface fit. tol=',1pe13.5)
cbug 9902 format (/ 'ax, ay, az= ',1p3e22.14 /
cbug & 'bx, by, bz= ',1p3e22.14 /
cbug & 'cx, cy, cz= ',1p3e22.14 /
cbug & 'dx, dy, dz= ',1p3e22.14 /
cbug & 'ex, ey, ez= ',1p3e22.14 /
cbug & 'fx, fy, fz= ',1p3e22.14 /
cbug & 'gx, gy, gz= ',1p3e22.14 /
cbug & 'hx, hy, hz= ',1p3e22.14 /
cbug & 'px, py, pz= ',1p3e22.14 )
cbug write ( 3, 9901) tol
cbug write ( 3, 9902) ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz,
cbug & ex, ey, ez, fx, fy, fz, gx, gy, gz, hx, hy, hz,
cbug & px, py, pz
cbugc***DEBUG ends.
c.... Initialize.
nerr = 0
c.... Find the coordinates of the nine points to be fit.
c.... Add a random vector to try to avoid a surface through the origin.
ntry = 0
110 ntry = ntry + 1
xran = ranf ()
yran = ranf ()
zran = ranf ()
x1 = ax + xran
y1 = ay + yran
z1 = az + zran
x2 = bx + xran
y2 = by + yran
z2 = bz + zran
x3 = cx + xran
y3 = cy + yran
z3 = cz + zran
x4 = dx + xran
y4 = dy + yran
z4 = dz + zran
x5 = ex + xran
y5 = ey + yran
z5 = ez + zran
x6 = fx + xran
y6 = fy + yran
z6 = fz + zran
x7 = gx + xran
y7 = gy + yran
z7 = gz + zran
x8 = hx + xran
y8 = hy + yran
z8 = hz + zran
x9 = px + xran
y9 = py + yran
z9 = pz + zran
cbugcbugc***DEBUG begins.
cbugcbug 9110 format ('aptqfit DEBUG ',1p3e20.12)
cbugcbug write ( 3, 9110) x1, y1, z1, x2, y2, z2, x3, y3, z3,
cbugcbug & x4, y4, z4, x5, y5, z5, x6, y6, z6,
cbugcbug & x7, y7, z7, x8, y8, z8, x9, y9, z9
cbugcbugc***DEBUG ends.
c.... Find the coefficients of qx.
abc(1,1) = x1
abc(2,1) = x2
abc(3,1) = x3
abc(4,1) = x4
abc(5,1) = x5
abc(6,1) = x6
abc(7,1) = x7
abc(8,1) = x8
abc(9,1) = x9
c.... Find the coefficients of qy.
abc(1,2) = y1
abc(2,2) = y2
abc(3,2) = y3
abc(4,2) = y4
abc(5,2) = y5
abc(6,2) = y6
abc(7,2) = y7
abc(8,2) = y8
abc(9,2) = y9
c.... Find the coefficients of qz.
abc(1,3) = z1
abc(2,3) = z2
abc(3,3) = z3
abc(4,3) = z4
abc(5,3) = z5
abc(6,3) = z6
abc(7,3) = z7
abc(8,3) = z8
abc(9,3) = z9
c.... Find the coefficients of qxy.
abc(1,4) = x1 * y1
abc(2,4) = x2 * y2
abc(3,4) = x3 * y3
abc(4,4) = x4 * y4
abc(5,4) = x5 * y5
abc(6,4) = x6 * y6
abc(7,4) = x7 * y7
abc(8,4) = x8 * y8
abc(9,4) = x9 * y9
c.... Find the coefficients of qyz.
abc(1,5) = y1 * z1
abc(2,5) = y2 * z2
abc(3,5) = y3 * z3
abc(4,5) = y4 * z4
abc(5,5) = y5 * z5
abc(6,5) = y6 * z6
abc(7,5) = y7 * z7
abc(8,5) = y8 * z8
abc(9,5) = y9 * z9
c.... Find the coefficients of qzx.
abc(1,6) = z1 * x1
abc(2,6) = z2 * x2
abc(3,6) = z3 * x3
abc(4,6) = z4 * x4
abc(5,6) = z5 * x5
abc(6,6) = z6 * x6
abc(7,6) = z7 * x7
abc(8,6) = z8 * x8
abc(9,6) = z9 * x9
c.... Find the coefficients of qxx.
abc(1,7) = x1 * x1
abc(2,7) = x2 * x2
abc(3,7) = x3 * x3
abc(4,7) = x4 * x4
abc(5,7) = x5 * x5
abc(6,7) = x6 * x6
abc(7,7) = x7 * x7
abc(8,7) = x8 * x8
abc(9,7) = x9 * x9
c.... Find the coefficients of qyy.
abc(1,8) = y1 * y1
abc(2,8) = y2 * y2
abc(3,8) = y3 * y3
abc(4,8) = y4 * y4
abc(5,8) = y5 * y5
abc(6,8) = y6 * y6
abc(7,8) = y7 * y7
abc(8,8) = y8 * y8
abc(9,8) = y9 * y9
c.... Find the coefficients of qzz.
abc(1,9) = z1 * z1
abc(2,9) = z2 * z2
abc(3,9) = z3 * z3
abc(4,9) = z4 * z4
abc(5,9) = z5 * z5
abc(6,9) = z6 * z6
abc(7,9) = z7 * z7
abc(8,9) = z8 * z8
abc(9,9) = z9 * z9
c.... Find the right-side constants. Assume qc = -1.
d(1) = 1.0
d(2) = 1.0
d(3) = 1.0
d(4) = 1.0
d(5) = 1.0
d(6) = 1.0
d(7) = 1.0
d(8) = 1.0
d(9) = 1.0
c.... Solve the nine simultaneous equations.
call aptsolv (abc, d, 9, work, 9, tol, xyz, cba, xyzm, xyzr,
& dm, dr, nerr)
if (nerr .eq. 1) then
if (ntry .ge. 10) then
go to 210
endif
nerr = 0
go to 110
endif
c.... Store data for the quadric surface.
qc = -1.0
qx = xyz(1)
qy = xyz(2)
qz = xyz(3)
qxy = xyz(4)
qyz = xyz(5)
qzx = xyz(6)
qxx = xyz(7)
qyy = xyz(8)
qzz = xyz(9)
c.... Remove the random vector.
rx = -xran
ry = -yran
rz = -zran
call apttris (1, rx, ry, rz, qc, qx, qy, qz,
& qxy, qyz, qzx, qxx, qyy, qzz, tol, nerr)
c.... Renormalize the coefficients to qc = -1.
if (qc .ne. 0.0) then
qx = -qx / qc
qy = -qy / qc
qz = -qz / qc
qxy = -qxy / qc
qyz = -qyz / qc
qzx = -qzx / qc
qxx = -qxx / qc
qyy = -qyy / qc
qzz = -qzz / qc
qc = -1.0
endif
210 continue ! End of loop over rows.
cbugc***DEBUG begins.
cbug 9910 format (/ 'aptqfit results. nerr = ',i2)
cbug 9911 format (/ 'nerr =',i2,'. 1 = singular matrix.',
cbug & ' 2 = large q residuals. 3 = large s residuals.')
cbug 9912 format (/ 'unknowns qc =',1pe22.14 /
cbug & 'qx, qy, qz =',1p3e22.14 /
cbug & 'qxy, qyz, qzx=',1p3e22.14 /
cbug & 'qxx, qyy, qzz=',1p3e22.14 )
cbug write ( 3, 9910) nerr
cbug if (nerr .ne. 0) then
cbug write ( 3, 9911) nerr
cbug endif
cbug write ( 3, 9912) qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz
cbugc***DEBUG ends.
return
c.... End of subroutine aptqfit. (+1 line.)
end
UCRL-WEB-209832