subroutine aptcniz (ar, az, br, bz, tol, qc, qx, qy, qz,
& qxy, qyz, qzx, qxx, qyy, qzz, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTCNIZ
c
c call aptcniz (ar, az, br, bz, tol, qc, qx, qy, qz,
c & qxy, qyz, qzx, qxx, qyy, qzz, nerr)
c
c Version: aptcniz Updated 1990 December 12 14:30.
c aptcniz Originated 1990 December 7 11:50.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To find the equation of a circular cone, symmetric around
c the z axis, that passes through the two distinct points
c a = (ar, az) and b = (br, bz), where ar = sqrt (ax**2 + ay**2)
c and br = sqrt (bx**2 + by**2). Points "a" and "b" are on the
c same side of the vertex of the cone if ar and br have the same
c sign, otherwise are on opposite sides. The coordinates (x,y,z)
c are along the three major orthogonal axes.
c
c The equation for the conical surface is as follows:
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 where qx = qy = qxy = qyz = qzx = 0, qxx = qyy, and
c qzz has the opposite sign of qxx.
c
c Degenerate case:
c
c If points "a" and "b" coincide, within the limits determined
c by tol, all coefficients will be -1.e99 and nerr = 1.
c
c Other special cases:
c
c If az = bz, the surface will be a z plane.
c
c If ar = br, the surface will be a right circular cylinder.
c
c If ar = br = 0, the surface will be a line on the z axis,
c and nerr = -1.
c
c Note: Use apttris and aptrois to translate and rotate the cone
c into any desired position and orientation.
c Use aptcois to find the equation of a cone on a given axis,
c with a given vertex angle and vertex axial coordinate.
c Use aptcyis to find the equation of a cylinder symmetric
c Around a major axis.
c Use aptelis to find the equation of an ellipsoid symmetric
c around a major axis.
c
c Input: ar, az, br, bz, tol.
c
c Output: qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz, nerr.
c
c Calls: aptvdic, aptvaxc
c
c
c Glossary:
c
c ar, az Input Coordinates r and z of first point in conical surface.
c If points "a" and "b" are on the same side of the
c vertex of the cone, the signs of ar and br must be
c the same. Otherwise, the signs must be different.
c
c br, bz Input Coordinates r and z of second point in conical surface.
c If points "a" and "b" are on the same side of the
c vertex of the cone, the signs of ar and br must be
c the same. Otherwise, the signs must be different.
c
c nerr Output Indicates an input error, if not zero:
c -1 if ar = br = 0, but az and bz are different, so
c there is no cone, but r = 0 and x**2 + y**2 = 0.
c 1 if points "a" and "b" are coincident.
c Indeterminate.
c
c qc Output Constant term in quadric cone equation.
c Zero only if the vertex of the cone is at the origin.
c qc = -(az * br - ar * bz)**2.
c
c qx,qy,qz Output Coefficients of x, y, z in quadric cone equation.
c qx = qy = 0.
c qz = 2.0 * (az * br - ar * bz) * (br - ar).
c
c q.. Output Coefficients of second-order terms in quadric
c cone 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 qxy = qyz = qzx = 0.
c qxx = qyy = (bz - az)**2.
c qzz = -(br - ar)**2.
c
c tol Input Numerical tolerance limit.
c On computers with 64-bit floating point numbers,
c recommend 1.e-5 to 1.e-11.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptcniz finding the quadric coefficients for a cone.',
cbug & ' tol=',1pe13.5 /
cbug & ' Cone passes through the points:')
cbug 9902 format (
cbug & ' ar,az= ',1p2e22.14 /
cbug & ' br,bz= ',1p2e22.14 )
cbug write ( 3, 9901) tol
cbug write ( 3, 9902) ar, az, br, bz
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 2-D vector from point "a" to point "b".
call aptvdic (ar, az, br, bz, 1, tol, abr, abz, ablen, nerra)
if (ablen .eq. 0.0) then
nerr = 1
go to 210
endif
if ((ar .eq. 0.0) .and. (br .eq. 0.0)) then
nerr = -1
endif
c.... Find the cross product of the 2-D vectors "a" and "b".
call aptvaxc (ar, az, br, bz, 1, tol, ablen, nerra)
c.... Find the coefficients of the equation of the cone.
qc = -ablen**2
qx = 0.0
qy = 0.0
qz = -2.0 * ablen * abr
qxy = 0.0
qyz = 0.0
qzx = 0.0
qxx = abz**2
qyy = qxx
qzz = -abr**2
210 continue
cbugc***DEBUG begins.
cbug 9904 format (/ 'aptcniz results. nerr=',i3 /
cbug & ' qc= ',1pe22.14 /
cbug & ' qx,qy,qz= ',1p3e22.14 /
cbug & ' qxy,qyz,qzx=',1p3e22.14 /
cbug & ' qxx,qyy,qzz=',1p3e22.14 )
cbug write ( 3, 9904) nerr, qc, qx, qy, qz,
cbug & qxy, qyz, qzx, qxx, qyy, qzz
cbugc***DEBUG ends.
return
c.... End of subroutine aptcniz. (+1 line.)
end
UCRL-WEB-209832