subroutine aptplcy (ax, ay, az, bx, by, bz,
& qc, qx, qy, qz, qxy, qyz, qzx,
& qxx, qyy, qzz, tol,
& rcyl, cx, cy, cz, dx, dy, dz,
& ipar, daxis, dpax, dpay, dpaz, dsurf1, dsurf2,
& dpcx, dpcy, dpcz, pcx, pcy, pcz,
& ppx, ppy, ppz, pinx, piny, pinz, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTPLCY
c
c call aptplcy (ax, ay, az, bx, by, bz,
c & qc, qx, qy, qz, qxy, qyz, qzx,
c & qxx, qyy, qzz, tol,
c & rcyl, cx, cy, cz, dx, dy, dz,
c & ipar, daxis, dpax, dpay, dpaz, dsurf1, dsurf2,
c & dpcx, dpcy, dpcz, pcx, pcy, pcz,
c & ppx, ppy, ppz, pinx, piny, pinz, nerr)
c
c Version: aptplcy Updated 1995 January 13 14:20.
c aptplcy Originated 1995 January 10 13:50.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: For the plane through the point a = (ax, ay, az) with normal
c vector b = (bx, by, bz), and the circular cylinder specified
c by the implicit quadric equation coefficients qc, ..., qzz:
c
c to find the radius rcyl of the cylinder,
c a point c = (cx, cy, cz) on its axis, and
c its axis vector d = (dx, dy, dz);
c
c if the axis of the cylinder is parallel (ipar = 1) to the plane,
c to find its distance from the plane da = abs (daxis);
c to find the minimum and maximum distances dsurf1 = da - rcyl
c and dsurf2 = da + rcyl from the plane to the surface of the
c cylinder;
c to find a point pp = (ppx, ppy, ppz) on the plane proximal to
c the cylinder, and the vector dpa = (dpax, dpay, dpaz) from point
c "pp" to point "c";
c to find a point pc = (pcx, pcy, pcz) on the cylinder proximal to
c the plane, and the vector dpc = (dpcx, dpcy, dpcz) from point
c "pp" to point "pc";
c if da does not exceed rcyl, to find the one (da = rcyl) or
c two (da < rcyl) points pin(1) = (pinx(1), piny(1), pinz(1)) and
c pin(2) = (pinx(2), piny(2), pinz(2)) on the line(s) of
c intersection of the plane with the surface of the cylinder
c (such lines will have the direction vector "d");
c
c if the axis of the cylinder is not parallel to the plane
c (ipar = 0), to find the point of intersection
c pp = (ppx, ppy, ppz) of the axis with the plane, and the four
c points pin(n) = (pinx(n), piny(n), pinz(n)), n = 1, 4, on the
c major and minor axes of the elliptical intersection of the
c plane with the circular cylinder.
c
c If the plane is defined only by its implicit quadric equation,
c
c f(x,y,z) = qc1 + qx1 * x + qy1 * y + qz1 * z = 0,
c
c a point "a" and the normal vector "b" may be found from:
c
c b = (bx, by, bz) = (qx1, qy1, qz1)
c b**2 = qx1**2 + qy1**2 + qz1**2
c a = (ax, ay, az) = -(qc1 / b**2) * (qx1, qy1, qz1)
c
c The circular cylinder is defined by the implicit quadric
c equation:
c
c g(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 which in standard form, with the axis through the origin and
c in the direction of the z axis, is
c g = -rcyl**2 + x**2 + y**2 = 0,
c where rcyl is the radius.
c
c If the circular cylinder is defined only by a point on the axis,
c the axis vector, and the radius, use APT subroutine aptclis to
c find the coefficients of the implicit quadric equation.
c
c Flag nerr indicates any input error.
c
c Input: ax, ay, az, bx, by, bz, qc, qx, qy, qz, qxy, qyz, qzx,
c qxx, qyy, qzz, tol.
c
c Output: rcyl, cx, cy, cz, dx, dy, dz, ipar, daxis, dpax, dpay, dpaz,
c dsurf1, dsurf2, dpcx, dpcy, dpcz, pcx, pcy, pcz,
c ppx, ppy, ppz, pinx, piny, pinz, nerr)
c
c Calls: aptaxis, aptlnpl, aptptpl, aptvang, aptvxun
c
c
c Glossary:
c
c ax,ay,az Input The x, y, z coordinates of point "a" in the plane.
c
c bx,by,bz Input The x, y, z components of the normal vector "b" of the
c plane.
c
c cx,cy,cz Output The x, y, z coordinates of point "c" on the axis of
c the circular cylinder.
c
c daxis Output Minimum distance between the plane and the axis of the
c cylinder. Zero if they intersect (ipar = 0), or the
c axis of the cylinder is entirely within the plane
c (ipar = 1). If positive, in the direction of the
c normal vector "b". If negative, in the opposite
c direction. The length of vector "dpa".
c
c dpax,y,z Output The smallest vector from the plane to the axis of the
c cylinder. If ipar = 0, a null vector.
c If ipar = 1, the vector from point "pp" on the plane
c to point "c" on the axis of the cylinder, with length
c daxis.
c Parallel (daxis > 0) or antiparallel (daxis < 0) to
c vector "b".
c
c dpcx,y,z Output The smallest normal vector from the plane to the
c surface of the cylinder, with length dsurf1.
c If ipar = 0, a null vector.
c If ipar = 1, the vector from point "pp" on the plane
c to point "pc" on the surface of the cylinder.
c The largest normal vector from the plane to
c the surface of the cylinder is 2 * "dpa" - "dpc",
c with length dsurf2.
c
c dsurf1 Output If ipar = 0, dsurf1 = 0.
c If ipar = 1, the minimum distance from the plane
c to the surface of the cylinder,
c dsurf1 = abs (daxis) - rcyl, negative if
c the plane and the cylinder overlap, zero if they are
c tangent, positive if they do not overlap.
c The length of vector "dpa".
c Moving the plane and the cylinder toward each other
c in the direction of vector "b" by distance dsurf1
c would make them tangent.
c
c dsurf2 Output If ipar = 0, dsurf2 = 0. If ipar = 1, the maximum
c distance from the plane to the surface of the
c cylinder, dsurf2 = abs (daxis) + rcyl.
c Moving the plane and the cylinder toward each other
c in the direction of vector "b" by distance dsurf2
c would make them tangent. Zero if ipar = 0.
c
c dx,dy,dz Input The x, y, z components of the unit axis vector "d" of
c the circular cylinder. If ipar = 1, this is the
c direction vector of the lines through points "pp",
c "pc" and "pin(1)" through "pin(4)".
c
c ipar Output Indicates that the axis of the circular cylinder is
c (ipar = 1) or is not (ipar = 0) parallel to the
c plane. See daxis.
c
c nerr Output Indicates an input error, if not zero.
c 1 if vector "b" has zero length.
c 2 if the implicit quadric equation coefficients do
c not define a circular cylinder.
c
c pcx,y,z Output If ipar = 1, the x, y, z coordinates of point "pc" on
c the line on the cylindrical surface nearest point
c "pp" on the plane, at distance dsurf1 from point "pp"
c in direction "dpa". The line has direction "d".
c The line on the cylindrical surface furthest from
c point "pp" is 2 * "c" - "pc", at distance dsurf2 from
c point "pp".
c
c pinx(1),... Output If ipar = 0, the points pin(1) through pin(4) are
c four points on the intersection of the plane with the
c surface of the circular cylinder. The points are on
c the ends of the major and minor axes of the
c elliptical intersection curve. Note that pinx, piny
c and pinz are arrays with a minimum size of 4.
c Point "pin(n)" has x, y, z coordinates pinx(n),
c piny(n), pinz(n), for n = 1, 4.
c
c pinx(1),... Output If ipar = 1, points pin(1) and pin(2) are points on
c any line(s) of intersection of the plane with the
c surface of the circular cylinder. The direction of
c each line is vector "d". There are no such lines if
c daxis > rcyl. Note that pinx, piny and pinz are
c arrays with a minimum size of 4.
c Point "pin(n)" has x, y, z coordinates pinx(n),
c piny(n), pinz(n), for n = 1, 4.
c
c ppx,y,z Output If ipar = 1, the x, y, z coordinates of the point "pp"
c on the plane nearest point "c" on the axis of the
c cylinder. If ipar = 0, the x, y, z coordinates of
c point "pp", the intersection of the plane with the
c axis of the cylinder, which is also the center of the
c ellipse of intersection of the plane with the surface
c of the cylinder.
c
c qc, ... Input Coefficients of the quadric equation for the cylinder:
c qc + qx * x + qy * y + qz * z +
c qxy * x * y + qyz * y * z + qzx * z * x +
c qxx * x**2 + qyy * y**2 + qxx * z**2 = 0
c
c rcyl Output The radius of the 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.... Dimensioned arguments.
dimension pinx(4) ! Coordinate x of a point.
dimension piny(4) ! Coordinate y of a point.
dimension pinz(4) ! Coordinate z of a point.
c.... Local variables.
common /laptplcy/ rotq(3,3) ! Rotation matrix to align cylinder.
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptplcy finding distance between plane and',
cbug & ' cylinder.',
cbug & ' tol=',1pe13.5)
cbug 9902 format (
cbug & ' ax, ay, az= ',1p3e22.14 /
cbug & ' bx, by, bz= ',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) ax, ay, az, bx, by, bz
cbug write ( 3, 9903) qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz
cbugc***DEBUG ends.
c.... Initialize.
nerr = 0
rcyl = -1.e99
cx = -1.e99
cy = -1.e99
cz = -1.e99
dx = -1.e99
dy = -1.e99
dz = -1.e99
pcx = -1.e99
pcy = -1.e99
pcz = -1.e99
daxis = -1.e99
dsurf1 = -1.e99
dsurf2 = -1.e99
ppx = -1.e99
ppy = -1.e99
ppz = -1.e99
dpax = -1.e99
dpay = -1.e99
dpaz = -1.e99
dpcx = -1.e99
dpcy = -1.e99
dpcz = -1.e99
do 110 n = 1, 4
pinx(n) = -1.e99
piny(n) = -1.e99
pinz(n) = -1.e99
110 continue
c.... Test for input errors.
call aptvunb (bx, by, bz, 1, tol, ubx, uby, ubz, blen, nerra)
if (blen .eq. 0.0) then ! Normal vector of plane is null.
nerr = 1
cbugc***DEBUG begins.
cbug 9904 format (/ 'aptplcy ERROR. Normal vector of plane is null.')
cbug write ( 3, 9904)
cbugc***DEBUG ends.
go to 710
endif
c.... Find the radius and axis of the cylinder, and a second point on the axis.
sc = qc
sx = qx
sy = qy
sz = qz
sxy = qxy
syz = qyz
szx = qzx
sxx = qxx
syy = qyy
szz = qzz
tolx = tol
if (tolx .le. 0.0) tolx = 1.e-11
call aptaxis (sc, sx, sy, sz, sxy, syz, szx, sxx, syy, szz, tolx,
& cx, cy, cz, rotq, ntype, nerraxis)
if ((ntype .ne. 7) .or.
& (nerraxis .ne. 0)) then
cbugc***DEBUG begins.
cbug 9905 format (/ 'aptplcy ERROR. Z axis of cylinder is null.')
cbug write ( 3, 9905)
cbugc***DEBUG ends.
nerr = 2
go to 710
endif
dx = rotq(3,1)
dy = rotq(3,2)
dz = rotq(3,3)
cx2 = cx + dx
cy2 = cy + dy
cz2 = cz + dz
rcyl = sqrt (-sc / sxx)
c.... Find a unit vector perpendicular to the normal and the cylinder axis.
call aptvxun (bx, by, bz, dx, dy, dz, 1, tol,
& ux, uy, uz, vlen, nerra)
if (vlen .eq. 0.0) then ! Normal is parallel to cylinder axis.
call aptvxun (bx, by, bz, 0.0, 0.0, 1.0, 1, tol,
& ux, uy, uz, vlen, nerra)
if (vlen .eq. 0.0) then ! Normal is parallel to z axis.
call aptvxun (bx, by, bz, 1.0, 0.0, 0.0, 1, tol,
& ux, uy, uz, vlen, nerra)
endif
endif
c.... Find any intersection of the axis of the cylinder with the plane.
call aptlnpl (cx, cy, cz, cx2, cy2, cz2,
& ax, ay, az, bx, by, bz, 1, tol,
& daxis, dint, fracps, ppx, ppy, ppz, ipar, nerra)
if (ipar .eq. 1) then ! Plane is parallel to cylinder axis.
cbugc***DEBUG begins.
cbug if (daxis .eq. 0.0) then
cbug 9906 format (/ 'aptplcy. Plane contains axis of cylinder.')
cbug write ( 3, 9906)
cbug else
cbug 9907 format (/ 'aptplcy. Plane is parallel to axis of cylinder.')
cbug write ( 3, 9907)
cbug endif
cbugc***DEBUG ends.
c.... Find the distances of minimum and maximum separation or overlap of the
c.... plane and the cylinder, when they are parallel.
dsurf1 = abs (daxis) - rcyl
dsurf2 = abs (daxis) + rcyl
disterr = tol * dsurf2
if (abs (dsurf1) .le. disterr) dsurf1 = 0.0
cbugc***DEBUG begins.
cbug if (dsurf1 .lt. 0.0) then
cbug if (daxis .ne. 0.0) then
cbug 9908 format (/ 'aptplcy. Plane overlaps the cylinder.')
cbug write ( 3, 9908)
cbug endif
cbug elseif (dsurf1 .eq. 0.0) then
cbug 9909 format (/ 'aptplcy. Plane is tangent to the cylinder.')
cbug write ( 3, 9909)
cbug else
cbug 9910 format (/ 'aptplcy. Plane does not contact the cylinder.')
cbug write ( 3, 9910)
cbug endif
cbugc***DEBUG ends.
c.... Find the point "pp" on the plane proximal to point "c" on the axis.
call aptptpl (cx, cy, cz, ax, ay, az, bx, by, bz, 1, tol,
& dpmin, ppx, ppy, ppz, itrun, nerra)
call aptvdis (ppx, ppy, ppz, cx, cy, cz, 1, tol,
& dpax, dpay, dpaz, dpa, nerra)
c.... Find the point "pc" on the cylinder proximal to point "a".
if (dpa .eq. 0.0) then ! The plane contains the axis.
pcx = ppx + dsurf1 * bx / blen
pcy = ppy + dsurf1 * by / blen
pcz = ppz + dsurf1 * bz / blen
else
pcx = ppx + dsurf1 * dpax / dpa
pcy = ppy + dsurf1 * dpay / dpa
pcz = ppz + dsurf1 * dpaz / dpa
endif
c.... Find the vector "dpc" from point "pp" to point "pc".
call aptvdis (ppx, ppy, ppz, pcx, pcy, pcz, 1, tol,
& dpcx, dpcy, dpcz, dpc, nerra)
c.... Find the two points on the intersection lines on the plane.
if (dsurf1 .le. 0.0) then ! Plane and cylinder surface intersect.
distpl = sqrt (-dsurf1 * dsurf2)
pinx(1) = ppx - distpl * ux
piny(1) = ppy - distpl * uy
pinz(1) = ppz - distpl * uz
pinx(2) = ppx + distpl * ux
piny(2) = ppy + distpl * uy
pinz(2) = ppz + distpl * uz
endif ! Plane and cylinder surface intersect.
go to 710
else ! Plane crosses axis of cylinder.
cbugc***DEBUG begins.
cbug 9911 format (/ 'aptplcy. Plane not parallel to axis of cylinder.')
cbug write ( 3, 9911)
cbugc***DEBUG ends.
dsurf1 = 0.0
dsurf2 = 0.0
dpax = 0.0
dpay = 0.0
dpaz = 0.0
dpcx = 0.0
dpcy = 0.0
dpcz = 0.0
c.... Find the two points on minor axis of elliptical curve of intersection.
pinx(1) = ppx - rcyl * ux
piny(1) = ppy - rcyl * uy
pinz(1) = ppz - rcyl * uz
pinx(2) = ppx + rcyl * ux
piny(2) = ppy + rcyl * uy
pinz(2) = ppz + rcyl * uz
c.... Find the two points on major axis of elliptical curve of intersection.
c.... Find the major axis of elliptical curve of intersection.
call aptvxun (bx, by, bz, ux, uy, uz, 1, tol,
& umx, umy, umz, vlen, nerra)
call aptvang (bx, by, bz, dx, dy, dz, 1, tol, costh, nerr)
dmaj = rcyl / abs (costh)
pinx(3) = ppx - dmaj * umx
piny(3) = ppy - dmaj * umy
pinz(3) = ppz - dmaj * umz
pinx(4) = ppx + dmaj * umx
piny(4) = ppy + dmaj * umy
pinz(4) = ppz + dmaj * umz
endif ! Tested whether parallel or not.
710 continue
cbugc***DEBUG begins.
cbug 9912 format (/ 'aptplcy results: ipar=',i3,12x,'nerr=',i3 /
cbug & ' rcyl= ',1pe22.14 /
cbug & ' cx,cy,cz= ',1p3e22.14 /
cbug & ' dx,dy,dz= ',1p3e22.14 /
cbug & ' daxis= ',1pe22.14 /
cbug & ' dpax,y,z= ',1p3e22.14 /
cbug & ' dsurf1,2= ',1p2e22.14 /
cbug & ' dpcx,y,z= =',1p3e22.14 /
cbug & ' pcx,pcy,pcz=',1p3e22.14 /
cbug & ' ppx,ppy,ppz=',1p3e22.14 /
cbug & ' pinx,y,z 1=',1p3e22.14 /
cbug & ' pinx,y,z 2=',1p3e22.14 /
cbug & ' pinx,y,z 3=',1p3e22.14 /
cbug & ' pinx,y,z 4=',1p3e22.14 )
cbug write ( 3, 9912) ipar, nerr, rcyl, cx, cy, cz, dx, dy, dz,
cbug & daxis, dpax, dpay, dpaz, dsurf1, dsurf2, dpcx, dpcy, dpcz,
cbug & pcx, pcy, pcz, ppx, ppy, ppz,
cbug & (pinx(n), piny(n), pinz(n), n = 1, 4)
cbugc***DEBUG ends.
return
c.... End of subroutine aptplcy. (+1 line)
end
UCRL-WEB-209832