subroutine aptqhyp (ax, ay, az, bx, by, bz, cx, cy, cz,
& dx, dy, dz, tol, qc, qx, qy, qz,
& qxy, qyz, qzx, qxx, qyy, qzz, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTQHYP
c
c call aptqhyp (ax, ay, az, bx, by, bz, cx, cy, cz,
c & dx, dy, dz, tol, qc, qx, qy, qz,
c & qxy, qyz, qzx, qxx, qyy, qzz, nerr)
c
c Version: aptqhyp Updated 1999 March 3 18:00.
c aptqhyp Originated 1999 February 22 18:00.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To find the quadric surface (a plane or hyperbolic paraboloid)
c that passes through the four vertices a = (ax, ay, az),
c b = (bx, by, bz), c = (cx, cy, cz) and d = (dx, dy, dz) of a
c quadrangle, and also contains the four edges of the quadrangle
c and its center point (at the intersection of the two lines
c joining the midpoints of opposite sides).
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 fewer than three of the vertex points are distinct, the
c method fails, and nerr = 1 is returned.
c
c Input: ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz, tol.
c
c Output: qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz, nerr.
c
c Calls: aptpfit, aptsolv, apttris, aptvdis
c
c
c Glossary:
c
c ax,ay,az Input The x, y and z coordinates of a point at vertex a of
c a quadrangle.
c
c bx,by,bz Input The x, y and z coordinates of a point at vertex b of
c a quadrangle.
c
c cx,cy,cz Input The x, y and z coordinates of a point at vertex c of
c a quadrangle.
c
c dx,dy,dz Input The x, y and z coordinates of a point at vertex d of
c a quadrangle.
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 2 if any right-hand side residual "dr" exceeds 1.e-6
c times the maximum d value. Results may be usable.
c 3 if any solution residual "xyzr" exceeds 1.e-6 times
c the maximum xyz value. Results may be usable.
c 4 if three or more vertex points coincide.
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, the x, y z components of the normal
c vector.
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. Zero for a plane.
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 (/ 'aptqhyp finding hyperboloid. 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 write ( 3, 9901) tol
cbug write ( 3, 9902) ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz
cbugc***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
c.... Find the coordinates of the nine points to be fit.
c.... Points 1-4 are the original vertices.
c.... Points 5-8 are the midpoints of the edges.
c.... Point 9 is the center of the quadrangle.
c.... Add a random vector to insure that coefficient qc is not zero.
avgx = 0.25 * (abs (ax) + abs (bx) + abs (cx) + abs (dx))
avgy = 0.25 * (abs (ay) + abs (by) + abs (cy) + abs (dy))
avgz = 0.25 * (abs (az) + abs (bz) + abs (cz) + abs (dz))
ntry = 0
110 ntry = ntry + 1
xran = ranf () * avgx
yran = ranf () * avgy
zran = ranf () * avgz
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 = 0.5 * (x1 + x2)
y5 = 0.5 * (y1 + y2)
z5 = 0.5 * (z1 + z2)
x6 = 0.5 * (x2 + x3)
y6 = 0.5 * (y2 + y3)
z6 = 0.5 * (z2 + z3)
x7 = 0.5 * (x3 + x4)
y7 = 0.5 * (y3 + y4)
z7 = 0.5 * (z3 + z4)
x8 = 0.5 * (x4 + x1)
y8 = 0.5 * (y4 + y1)
z8 = 0.5 * (z4 + z1)
x9 = 0.25 * (x1 + x2 + x3 + x4)
y9 = 0.25 * (y1 + y2 + y3 + y4)
z9 = 0.25 * (z1 + z2 + z3 + z4)
cbugcbugc***DEBUG begins.
cbugcbug 9110 format ('aptqhyp 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 qyy.
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)
c.... See if a good fit was obtained.
if (nerr .ne. 0) then ! No good solution.
if (ntry .ge. 10) then ! Did 10 tries.
c.... See if any three of the vertex points are unique.
c.... Find the edge lengths.
call aptvdis (ax, ay, az, bx, by, bz, 1, tol,
& abx, aby, abz, ablen, nerra)
call aptvdis (bx, by, bz, cx, cy, cz, 1, tol,
& bcx, bcy, bcz, bclen, nerra)
call aptvdis (cx, cy, cz, dx, dy, dz, 1, tol,
& cdx, cdy, cdz, cdlen, nerra)
call aptvdis (dx, dy, dz, ax, ay, az, 1, tol,
& dax, day, daz, dalen, nerra)
call aptvdis (ax, ay, az, cx, cy, cz, 1, tol,
& acx, acy, acz, aclen, nerra)
call aptvdis (bx, by, bz, dx, dy, dz, 1, tol,
& bdx, bdy, bdz, bdlen, nerra)
nver = 1 ! Point "a".
x1 = ax
y1 = ay
z1 = az
if (ablen .ne. 0.0) then
nver = 2 ! Points "a" and "b".
x2 = bx
y2 = by
z2 = bz
if ((aclen .ne. 0.0) .and.
& (bclen .ne. 0.0)) then
nver = 3 ! Points "a", "b" and "c".
x3 = cx
y3 = cy
z3 = cz
elseif ((dalen .ne. 0.0) .and.
& (bdlen .ne. 0.0)) then
nver = 3 ! Points "a", "b" and "d".
x3 = dx
y3 = dy
z3 = dz
endif
elseif (aclen .ne. 0.0) then
nver = 2 ! Points "a" and "c".
x2 = cx
y2 = cy
z2 = cz
if ((dalen .ne. 0.0) .and.
& (cdlen .ne. 0.0)) then
nver = 3 ! Points "a", "c" and "d".
x3 = dx
y3 = dy
z3 = dz
endif
endif
if (nver .eq. 3) then ! Three distinct points exist.
nerr = 0
c.... Find the coefficients of the equation of the plane.
call aptpfit (x1, y1, z1, x2, y2, z2, x3, y3, z3, tol,
& sc, sx, sy, sz, nerra)
if (nerra .ne. 0) then ! Still no good solution.
nerr = 4
go to 210
endif
qc = sc
qx = sx
qy = sy
qz = sz
qxy = 0.0
qyz = 0.0
qzx = 0.0
qxx = 0.0
qyy = 0.0
qzz = 0.0
endif ! Three distinct points exist.
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.... Normalize the coefficients.
fact = 1.0
if (qc .ne. 0.0) then
fact = -1.0 / qc
elseif (qx .ne. 0.0) then
fact = 1.0 / qx
elseif (qy .ne. 0.0) then
fact = 1.0 / qy
elseif (qz .ne. 0.0) then
fact = 1.0 / qz
elseif (qxy .ne. 0.0) then
fact = 1.0 / qxy
elseif (qyz .ne. 0.0) then
fact = 1.0 / qyz
elseif (qzx .ne. 0.0) then
fact = 1.0 / qzx
elseif (qxx .ne. 0.0) then
fact = 1.0 / qxx
elseif (qyy .ne. 0.0) then
fact = 1.0 / qyy
elseif (qzz .ne. 0.0) then
fact = 1.0 / qzz
endif
qc = qc * fact
qx = qx * fact
qy = qy * fact
qz = qz * fact
qxy = qxy * fact
qyz = qyz * fact
qzx = qzx * fact
qxx = qxx * fact
qyy = qyy * fact
qzz = qzz * fact
210 continue ! End of loop over rows.
cbugc***DEBUG begins.
cbug 9910 format (/ 'aptqhyp 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 aptqhyp. (+1 line.)
end
UCRL-WEB-209832