subroutine aptcycy (qc1, qx1, qy1, qz1, qxy1, qyz1, qzx1,
& qxx1, qyy1, qzz1,
& qc2, qx2, qy2, qz2, qxy2, qyz2, qzx2,
& qxx2, qyy2, qzz2, tol,
& rcyl1, rcyl2, ax1, ay1, az1, ax2, ay2, az2,
& bx, by, bz, daxis,
& cx1, cy1, cz1, cx2, cy2, cz2,
& dx, dy, dz, dsurf, ipar, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTCYCY
c
c call aptcycy (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 rcyl1, rcyl2, ax1, ay1, az1, ax2, ay2, az2,
c bx, by, bz, daxis,
c cx1, cy1, cz1, cx2, cy2, cz2,
c dx, dy, dz, dsurf, ipar, nerr)
c
c Version: aptcycy Updated 1995 December 9 16:35.
c aptcycy Originated 1995 January 3 13:40.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: For the two circular cylinders specified by the coefficients
c qc1, ..., qzz1 and qc2, ..., qzz2, to find the radii rcyl1 and
c rcyl2 of the two cylinders; to find the proximal points
c a1 = (ax1, ay1, az1) and a2 = (ax2, ay2, az2) on the two axes,
c and the vector b = (bx, by, bz) and distance daxis between them;
c to find the points c1 = (cx1, cy1, cz1) and c2 = (cx2, cy2, cz2)
c of minimum separation or maximum overlap of the two cylinders,
c and the vector d = (dx, dy, dz) and distance dsurf between them;
c and to find (indicator ipar) if the two axes are parallel or
c not, and if the two cylinders are coincident or not. Three
c other larger distances of separation or overlap, all in
c direction d, may be constructed from rcyl1, rcyl2 and daxis:
c daxis - rcyl1 - rcyl2 = dsurf (min dist to external tangency)
c daxis + rcyl1 + rcyl2 (max dist to external tangency)
c daxis - rcyl1 + rcyl2 (distance to internal tangency)
c daxis + rcyl1 - rcyl2 (distance to internal tangency)
c
c The circular cylinders are defined by the implicit quadric
c equations:
c
c f1(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 f2(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 f = -r**2 + x**2 + y**2 = 0,
c where r is the radius.
c
c Use APT subroutine aptclis to find the coefficients of the
c implicit quadric equation, if either 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: rcyl1, rcyl2, ax1, ay1, az1, ax2, ay2, az2, bx, by, bz, daxis,
c cx1, cy1, cz1, cx2, cy2, cz2, dx, dy, dz, dsurf, ipar, nerr.
c
c Calls: aptaxis, aptlnln, aptmopv, apttran, aptvdis, aptvxun
c
c
c Glossary:
c
c ax1, ... Output The x, y, z coordinates of the point on the axis of
c the first cylinder nearest the axis of the second
c cylinder.
c
c ax2, ... Output The x, y, z coordinates of the point on the axis of
c the second cylinder nearest the axis of the first
c cylinder.
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 first
c cylinder by this amount would make its axis
c intersect the axis of the second cylinder.
c 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 first cylinder,
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 second cylinder,
c c2 = (cx2, cy2, cz2).
c
c daxis Output Minimum distance between the axes of the two cylinders.
c Zero if they intersect. If daxis = 0 and ipar = 1,
c the two axes coincide, so the two cylinders are
c concentric. If daxis = 0 and ipar = 2, the two
c cylinders coincide.
c
c dsurf Output Distance of minimum separation or maximum overlap
c of the surfaces of the two cylinders.
c dsurf = daxis - rcyl1 - rcyl2
c Negative if the two cylinders overlap.
c Zero if the two cylinders are externally tangent.
c Positive if the two cylinders do not overlap.
c Moving the cylinders toward each other in the
c direction of vector d by distance dsurf would make
c them externally tangent. So would the larger
c distance daxis + rcyl1 + rcyl2. Movement by the
c following distances would make the two circular
c cylinders internally tangent: daxis - rcyl1 + rcyl2
c and daxis + rcyl1 - rcyl2.
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 first
c cylinder by this amount would make it externally
c tangent to the second cylinder.
c Vector d is parallel to vector b.
c
c ipar Output Zero if the axes of the two cylinders are not parallel.
c 1 if the axes of the two cylinders are parallel.
c If ipar = 1 and daxis = 0, the two axes coincide,
c so the cylinders are concentric.
c 2 if the two cylinders coincide.
c
c nerr Output Indicates an input error, if not zero.
c 1 if either of the quadric surfaces is not a circular
c cylinder.
c
c qc1, ... Input Coefficients of the quadric equation for the first
c cylinder:
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
c qc2, ... Input Coefficients of the quadric equation for the second
c 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 rcyl1 Output The radius of the first circular cylinder.
c
c rcyl2 Output The radius of the second circular cylinder.
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 /laptcycy/ rotq(3,3) ! Rotation matrix to align cylinder.
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptcycy finding distance between two cylinders.',
cbug & ' tol=',1pe13.5)
cbug 9902 format (
cbug & ' Cyl 1: 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 & ' Cyl 2: 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 radius of the first cylinder, and a line on the axis.
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. 7) .or.
& (nerraxis .ne. 0)) then
nerr = 1
go to 710
endif
vx1 = rotq(3,1)
vy1 = rotq(3,2)
vz1 = rotq(3,3)
xend1 = xcen1 + vx1
yend1 = ycen1 + vy1
zend1 = zcen1 + vz1
rcyl1 = sqrt (-ac / axx)
c.... Find the radius of the second 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 = 1
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
rcyl2 = sqrt (-ac / axx)
c.... Find the proximal points on each axis, and the distance of separation.
call aptlnln (xcen1, ycen1, zcen1, xend1, yend1, zend1,
& xcen2, ycen2, zcen2, xend2, yend2, zend2,
& 1, tol, dpmin, fracab, fraccd,
& ax1, ay1, az1, ax2, ay2, az2, itrun, ipar, 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.... circular cylinders.
dsurf = daxis - rcyl1 - rcyl2
disterr = tol * (abs (daxis) + abs (rcyl1) + abs (rcyl2))
if (abs (dsurf) .le. disterr) dsurf = 0.0
if (daxis .eq. 0.0) then ! Cylinder axes intersect or coincide.
if (ipar .eq. 1) then ! Cylinder axes coincide.
rcylerr = tol * (rcyl1 + rcyl2)
if (abs (rcyl2 - rcyl1) .le. rcylerr) then ! Cylinders coincide.
ipar = 2
endif
c.... Find two arbitrary proximal points on the circle of proximity.
cx1 = rcyl1
cy1 = 0.0
cz1 = 0.0
cx2 = -rcyl2
cy2 = 0.0
cz2 = 0.0
call aptmopv (rotq, 1, 0., 0., 0., cx1, cy1, cz1, 1, tol,
& nerr)
call aptmopv (rotq, 1, 0., 0., 0., cx2, cy2, cz2, 1, tol,
& nerr)
dcenx1 = -xcen1
dceny1 = -ycen1
dcenz1 = -zcen1
call apttran (dcenx1, dceny1, dcenz1, cx1, cy1, cz1, 1, tol,
& nerr)
call apttran (dcenx1, dceny1, dcenz1, cx2, cy2, cz2, 1, tol,
& nerr)
c.... Find the vector "d" between the arbitrary proximal surface points.
call aptvdis (cx1, cy1, cz1, cx2, cy2, cz2, 1, tol,
& dx, dy, dz, dist, nerra)
go to 710
endif
c.... Find the unit vector perpendicular to both axes.
call aptvxun (vx1, vy1, vz1, vx2, vy2, vz2, 1, tol,
& vx12, vy12, vz12, vlen, nerra)
cx1 = ax1 + rcyl1 * vx12
cy1 = ay1 + rcyl1 * vy12
cz1 = az1 + rcyl1 * vz12
cx2 = ax2 - rcyl2 * vx12
cy2 = ay2 - rcyl2 * vy12
cz2 = az2 - rcyl2 * 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 ! Cylinder axes do not intersect.
cx1 = ax1 + rcyl1 * bx / daxis
cy1 = ay1 + rcyl1 * by / daxis
cz1 = az1 + rcyl1 * bz / daxis
cx2 = ax2 - rcyl2 * bx / daxis
cy2 = ay2 - rcyl2 * by / daxis
cz2 = az2 - rcyl2 * 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 (/ 'aptcycy results: ipar=',i3,13x,' nerr=',i3 /
cbug & ' rcyl1,rcyl2=',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) ipar, nerr, rcyl1, rcyl2,
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 aptcycy. (+1 line)
end
UCRL-WEB-209832