subroutine aptspcy (qc1, qx1, qy1, qz1, qxy1, qyz1, qzx1,
& qxx1, qyy1, qzz1,
& qc2, qx2, qy2, qz2, qxy2, qyz2, qzx2,
& qxx2, qyy2, qzz2, tol,
& rsph, rcyl, ax1, ay1, az1, ax2, ay2, az2,
& bx, by, bz, daxis,
& cx1, cy1, cz1, cx2, cy2, cz2,
& dx, dy, dz, dsurf, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTSPCY
c
c call aptspcy (qc1, qx1, qy1, qz1, qxy1, qyz1, qzx1,
c qxx1, qyy1, qzz1,
c qc2, qx2, qy2, qz2, qxy2, qyz2, qzx2,
c qxx2, qyy2, qzz2, tol,
c rsph, rcyl, ax1, ay1, az1, ax2, ay2, az2,
c bx, by, bz, daxis,
c cx1, cy1, cz1, cx2, cy2, cz2,
c dx, dy, dz, dsurf, nerr)
c
c Version: aptspcy Updated 1995 January 5 12:00.
c aptspcy Originated 1995 January 5 12:00.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: For the sphere specified by the coefficients qc1, ..., qzz1
c and the circular cylinder specified by the coefficients
c qc2, ..., qzz2, to find the radii rsph of the sphere and
c rcyl of the cylinder; to find the center of the sphere
c a1 = (ax1, ay1, az1) and the proximal point a2 = (ax2, ay2, az2)
c on the axis of the cylinder, and the vector b = (bx, by, bz) and
c distance daxis between them; to find the proximal points
c c1 = (cx1, cy1, cz1) on the sphere and c2 = (cx2, cy2, cz2) on
c the circular cylinder of minimum separation or maximum overlap
c of their surfaces, and the vector d = (dx, dy, dz) and distance
c dsurf between them.
c Three other larger distances of separation or overlap, all in
c direction d, may be constructed from rsph, rcyl and daxis:
c daxis - rsph - rcyl = dsurf (min dist to external tangency)
c daxis + rsph + rcyl (max dist to external tangency)
c daxis - rsph + rcyl (distance to internal tangency)
c daxis + rsph - rcyl (distance to internal tangency)
c
c The sphere is defined by the implicit quadric equation:
c
c fs(x,y,z) = qc1 + qx1 * x + qy1 * y + qz1 * z +
c qxy1 * x * y + qyz1 * y * z + qzx1 * z * x +
c qxx1 * x**2 + qyy1 * y**2 + qzz1 * z**2 = 0,
c
c which in standard form, with the center at the origin, is
c fs = -rsph**2 + x**2 + y**2 + z**2 = 0,
c where rsph is the radius.
c
c If the sphere is centered at point a1 = (ax1, ay1, az1), with
c radius rsph, the coefficients of the implicit quadric equation
c are:
c qc1 = -rsph**2 + ax1**2 + ay1**2 + az1**2
c qx1 = -2.0 * ax1 qxy1 = 0.0 qxx1 = 1.0
c qy1 = -2.0 * ay1 qyz1 = 0.0 qyy1 = 1.0
c qz1 = -2.0 * az1 qzx1 = 0.0 qzz1 = 1.0
c
c The circular cylinder is defined by the equation:
c
c fc(x,y,z) = qc2 + qx2 * x + qy2 * y + qz2 * z +
c qxy2 * x * y + qyz2 * y * z + qzx2 * z * x +
c qxx2 * x**2 + qyy2 * y**2 + qzz2 * z**2 = 0,
c
c which in standard form, with the axis through the origin and
c in the direction of the z axis, is
c fc = -rcyl**2 + x**2 + y**2 = 0,
c where rcyl is the radius.
c
c Use APT subroutine aptclis to find the coefficients of the
c implicit quadric equation, if the circular cylinder is
c defined only by a point on the axis, the axis vector, and the
c radius.
c
c Flag nerr indicates any input error.
c
c Input: qc1, qx1, qy1, qz1, qxy1, qyz1, qzx1, qxx1, qyy1, qzz1,
c qc2, qx2, qy2, qz2, qxy2, qyz2, qzx2, qxx2, qyy2, qzz2, tol.
c
c Output: rsph, rcyl, ax1, ay1, az1, ax2, ay2, az2, bx, by, bz, daxis,
c cx1, cy1, cz1, cx2, cy2, cz2, dx, dy, dz, dsurf, nerr.
c
c Calls: aptaxis, aptptln, aptvdis, aptvxun
c
c
c Glossary:
c
c ax1, ... Output The x, y, z coordinates of the center of the sphere.
c
c ax2, ... Output The x, y, z coordinates of the point on the axis of
c the cylinder nearest the center of the sphere.
c
c bx,by,bz Output The x, y, z components of the vector from point
c a1 = (ax1, ay1, az1) to point a2 = (ax2, ay2, az2).
c The length of vector b is daxis. Moving the sphere
c by this amount would put its center on the axis of
c the cylinder. Vector b is parallel to vector d.
c
c cx1, ... Output The x, y, z coordinates of the point of minimum
c distance or maximum overlap on the sphere,
c c1 = (cx1, cy1, cz1).
c
c cx2, ... Output The x, y, z coordinates of the point of minimum
c distance or maximum overlap on the cylinder,
c c2 = (cx2, cy2, cz2).
c
c daxis Output Minimum distance between the center of the sphere and
c the axis of the cylinder.
c
c dsurf Output Distance of minimum separation or maximum overlap
c of the surfaces of the sphere and the surface of the
c cylinder: dsurf = daxis - rsph - rcyl.
c Negative if the sphere and the cylinder overlap.
c Zero if the sphere and the cylinder are externally
c tangent.
c Positive if the sphere and the cylinder do not
c overlap.
c Moving the sphere and the cylinder toward each other
c in the direction of vector d by distance dsurf would
c make them externally tangent. So would the larger
c distance daxis + rsph + rcyl. Movement by the
c following distances would make the sphere and the
c cylinder internally tangent: daxis - rsph + rcyl
c and daxis + rsph - rcyl.
c
c dx,dy,dz Output The x, y, z components of the vector from point
c c1 = (cx1, cy1, cz1) to point c2 = (cx2, cy2, cz2).
c The length of vector d is dsurf. Moving the sphere
c by this amount would make it externally tangent to
c the cylinder. Vector d is parallel to vector b.
c
c nerr Output Indicates an input error, if not zero.
c 1 if the first quadric surface is not a sphere.
c 2 if the second quadric surface is not a circular
c cylinder.
c
c qc1, ... Input Coefficients of the quadric equation for the sphere:
c qc1 + qx1 * x + qy1 * y + qz1 * z +
c qxy1 * x * y + qyz1 * y * z + qzx1 * z * x +
c qxx1 * x**2 + qyy1 * y**2 + qxx1 * z**2 = 0
c Given the center a1 = (ax1, ay1, az1) and the
c radius rsph:
c qc1 = -rsph**2 + ax1**2 + ay1**2 + az1**2
c qx1 = -2.0 * ax1 qxy1 = 0.0 qxx1 = 1.0
c qy1 = -2.0 * ay1 qyz1 = 0.0 qyy1 = 1.0
c qz1 = -2.0 * az1 qzx1 = 0.0 qzz1 = 1.0
c
c qc2, ... Input Coefficients of the quadric equation for the cylinder:
c qc2 + qx2 * x + qy2 * y + qz2 * z +
c qxy2 * x * y + qyz2 * y * z + qzx2 * z * x +
c qxx2 * x**2 + qyy2 * y**2 + qxx2 * z**2 = 0
c
c rcyl Output The radius of the circular cylinder.
c
c rsph Output The radius of the sphere.
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.
common /laptspcy/ rotq(3,3) ! Rotation matrix to align cylinder.
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptspcy finding distance between sphere and',
cbug & ' cylinder.',
cbug & ' tol=',1pe13.5)
cbug 9902 format (
cbug & ' Sphere qc=',1pe22.14 /
cbug & ' qx,qy,qz= ',1p3e22.14 /
cbug & ' qxy,qyz,qzx=',1p3e22.14 /
cbug & ' qxx,qyy,qzz=',1p3e22.14)
cbug 9903 format (
cbug & ' Cylinder qc=',1pe22.14 /
cbug & ' qx,qy,qz= ',1p3e22.14 /
cbug & ' qxy,qyz,qzx=',1p3e22.14 /
cbug & ' qxx,qyy,qzz=',1p3e22.14)
cbug write ( 3, 9901) tol
cbug write ( 3, 9902) qc1, qx1, qy1, qz1, qxy1, qyz1, qzx1,
cbug & qxx1, qyy1, qzz1
cbug write ( 3, 9903) qc1, qx1, qy1, qz1, qxy1, qyz1, qzx1,
cbug & qxx1, qyy1, qzz1
cbugc***DEBUG ends.
c.... Initialize.
nerr = 0
daxis = -1.e99
ax1 = -1.e99
ay1 = -1.e99
az1 = -1.e99
ax2 = -1.e99
ay2 = -1.e99
az2 = -1.e99
bx = -1.e99
by = -1.e99
bz = -1.e99
dsurf = -1.e99
cx1 = -1.e99
cy1 = -1.e99
cz1 = -1.e99
cx2 = -1.e99
cy2 = -1.e99
cz2 = -1.e99
dx = -1.e99
dy = -1.e99
dz = -1.e99
c.... Find the center and radius of the sphere.
ac = qc1
ax = qx1
ay = qy1
az = qz1
axy = qxy1
ayz = qyz1
azx = qzx1
axx = qxx1
ayy = qyy1
azz = qzz1
tolx = tol
if (tolx .le. 0.0) tolx = 1.e-11
call aptaxis (ac, ax, ay, az, axy, ayz, azx, axx, ayy, azz, tolx,
& xcen1, ycen1, zcen1, rotq, ntype, nerraxis)
if ((ntype .ne. 16) .or.
& (nerraxis .ne. 0)) then
nerr = 1
go to 710
endif
ax1 = xcen1
ay1 = ycen1
az1 = zcen1
rsph = sqrt (-ac / axx)
c.... Find the radius of the cylinder, and a line on the axis.
ac = qc2
ax = qx2
ay = qy2
az = qz2
axy = qxy2
ayz = qyz2
azx = qzx2
axx = qxx2
ayy = qyy2
azz = qzz2
tolx = tol
if (tolx .le. 0.0) tolx = 1.e-11
call aptaxis (ac, ax, ay, az, axy, ayz, azx, axx, ayy, azz, tolx,
& xcen2, ycen2, zcen2, rotq, ntype, nerraxis)
if ((ntype .ne. 7) .or.
& (nerraxis .ne. 0)) then
nerr = 2
go to 710
endif
vx2 = rotq(3,1)
vy2 = rotq(3,2)
vz2 = rotq(3,3)
xend2 = xcen2 + vx2
yend2 = ycen2 + vy2
zend2 = zcen2 + vz2
rcyl = sqrt (-ac / axx)
c.... Find the proximal point on the axis of the cylinder, and the distance
c.... from the center of the sphere to the axis of the cylinder.
call aptptln (ax1, ay1, az1,
& xcen2, ycen2, zcen2, xend2, yend2, zend2,
& 1, tol, 0,
& dpmin, fdmin, ax2, ay2, az2, nlim, itrun, nerra)
call aptvdis (ax1, ay1, az1, ax2, ay2, az2, 1, tol,
& bx, by, bz, daxis, nerra)
c.... Find the distance of minimum separation or maximum overlap of the
c.... sphere and the cylinder.
dsurf = daxis - rsph - rcyl
disterr = tol * (abs (daxis) + abs (rsph) + abs (rcyl))
if (abs (dsurf) .le. disterr) dsurf = 0.0
if (daxis .eq. 0.0) then ! Center of sphere on cylinder axis.
c.... The direction of proximity is arbitrary. Any vector perpendicular to
c.... the axis of the cylinder will do. Find the unit vector perpendicular
c to the z axis and the axis of the cylinder. If none, use x vector.
call aptvxun (0., 0., 1., vx2, vy2, vz2, 1, tol,
& vx12, vy12, vz12, vlen, nerra)
if (vlen .eq. 0.0) then ! Cylinder on z axis, use x vector.
vx12 = 1.0
vy12 = 0.0
vz12 = 0.0
endif
cx1 = ax1 + rsph * vx12
cy1 = ay1 + rsph * vy12
cz1 = az1 + rsph * vz12
cx2 = ax2 - rcyl * vx12
cy2 = ay2 - rcyl * vy12
cz2 = az2 - rcyl * vz12
c.... Find the vector "d" between the proximal surface points.
call aptvdis (cx1, cy1, cz1, cx2, cy2, cz2, 1, tol,
& dx, dy, dz, dist, nerra)
go to 710
else ! Center of sphere not on cylinder axis.
cx1 = ax1 + rsph * bx / daxis
cy1 = ay1 + rsph * by / daxis
cz1 = az1 + rsph * bz / daxis
cx2 = ax2 - rcyl * bx / daxis
cy2 = ay2 - rcyl * by / daxis
cz2 = az2 - rcyl * bz / daxis
c.... Find the vector "d" between the proximal points.
call aptvdis (cx1, cy1, cz1, cx2, cy2, cz2, 1, tol,
& dx, dy, dz, dist, nerra)
endif ! Tested whether axes intersect or not.
710 continue
cbugc***DEBUG begins.
cbug 9908 format (/ 'aptspcy results: nerr=',i3 /
cbug & ' rsph,rcyl= ',1p2e22.14 /
cbug & ' ax1,ay1,az1=',1p3e22.14 /
cbug & ' ax2,ay2,az2=',1p3e22.14 /
cbug & ' bx,by,bz= ',1p3e22.14 /
cbug & ' daxis= ',1pe22.14 /
cbug & ' cx1,cy1,cz1=',1p3e22.14 /
cbug & ' cx2,cy2,cz2=',1p3e22.14 /
cbug & ' dx,dy,dz= ',1p3e22.14 /
cbug & ' dsurf= ',1pe22.14 )
cbug write ( 3, 9908) nerr, rsph, rcyl,
cbug & ax1, ay1, az1, ax2, ay2, az2, bx, by, bz, daxis,
cbug & cx1, cy1, cz1, cx2, cy2, cz2, dx, dy, dz, dsurf
cbug dpp = daxis + rsph + rcyl
cbug errd = tol * dpp
cbug dmp = daxis - rsph + rcyl
cbug if (abs (dmp) .le. errd) dmp = 0.0
cbug dpm = daxis + rsph - rcyl
cbug if (abs (dpm) .le. errd) dpm = 0.0
cbug 9909 format (' dpp,dmp,dpm=',1p3e22.14 )
cbug write ( 3, 9909) dpp, dmp, dpm
cbugc***DEBUG ends.
return
c.... End of subroutine aptspcy. (+1 line)
end
UCRL-WEB-209832