subroutine aptplis (kpl, apl, qc, qx, qy, qz, qxy, qyz, qzx,
& qxx, qyy, qzz, np, tol, sc, sx, sy, sz,
& sxy, syz, szx, sxx, syy, szz, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTPLIS
c
c call aptplis (kpl, apl, qc, qx, qy, qz, qxy, qyz, qzx,
c & qxx, qyy, qzz, np, tol, sc, sx, sy, sz,
c & sxy, syz, szx, sxx, syy, szz, nerr)
c
c Version: aptplis Updated 1990 December 14 11:30.
c aptplis Originated 1990 December 14 11:30.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To find the intersections of a major plane, at x = apl,
c y = apl, or z = apl, for kpl = 1, 2 or 3, respectively,
c with each of the np quadric surfaces given by the equation:
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 The intersection is represented by the equations:
c
c u = apl (u = x, y or z for kpl = 1, 2 or 3, respectively),
c
c g(x, y, z) = sc + sx * x + sy * y + sz * z +
c sxy * x * y + syz * y * z + szx * z * x +
c sxx * x**2 + syy * y**2 + szz * z**2 = 0 ,
c
c where the coefficents of terms containing u are zero.
c
c General cases:
c
c The intersection equation may have no real solutions (no
c intersection), a tangent point or line solution, or
c be the equation of a real circle, ellipse, one or two
c hyperbolas, a parabola, a line, two coincident lines,
c two parallel lines, or two intersecting lines, in the major
c plane.
c
c Indeterminate case:
c
c If the coefficients of the equation of the quadric surface are
c all zero, the coefficients of the intersection equation will
c also all be zero.
c
c Input: kpl, apl, qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz, tol.
c
c Output: sc, sx, sy, sz, sxy, syz, szx, sxx, syy, szz, nerr.
c
c Calls: apttris
c
c Glossary:
c
c apl Input Axial coordinate of the major plane intersecting the
c quadric surfaces.
c
c kpl Input Indicates orientation of the major plane intersecting
c the quadric surfaces:
c 1 for the plane perpendicular to x axis at x = apl.
c 2 for the plane perpendicular to y axis at y = apl.
c 3 for the plane perpendicular to z axis at z = apl.
c
c nerr Output Indicates an input error, if not zero:
c 1 if np is not positive.
c 2 if kpl is not 1, 2 or 3.
c
c np Input Number of quadric surfaces. Size of arrays qc, qx, qy,
c qxy, qyz, qzx, qxx, qyy, qzz, sc, sx, sy, sz,
c sxy, syz, szx, sxx, syy, szz.
c
c qc In/Out Constant term in the equation of the quadric surface.
c Size np.
c
c qx,qy,qz In/Out Coefficients of x, y, z in the equation of the
c quadric surface. Size np.
c
c q.. In/Out Coefficients of the second-order terms in the equation
c of the quadric surface:
c qxy, qyz, qzx, qxx, qyy, qzz. Size np.
c Coefficients of x*y, y*z, z*x, x**2, y**2, z**2,
c respectively.
c
c sc In/Out Constant term in the equation of the intersection
c curve. Size np.
c
c sx,sy,sz In/Out Coefficients of x, y, z in the equation of the
c intersection curve. Size np.
c
c s.. In/Out Coefficients of the second-order terms in the equation
c of the intersection curve:
c sxy, syz, szx, sxx, syy, szz. Size np.
c Coefficients of x*y, y*z, z*x, x**2, y**2, z**2,
c respectively.
c
c tol Input Numerical tolerance limit.
c If zero, 1.e-11 will be used instead.
c On Cray computers, recommend 1.e-11.
c On 32-bit computers, recommend 1.e-6.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c---- Coefficient of term 1.
dimension qc (1)
c---- Coefficient of term x.
dimension qx (1)
c---- Coefficient of term y.
dimension qy (1)
c---- Coefficient of term z.
dimension qz (1)
c---- Coefficient of term x * y.
dimension qxy (1)
c---- Coefficient of term y * z.
dimension qyz (1)
c---- Coefficient of term z * x.
dimension qzx (1)
c---- Coefficient of term x**2.
dimension qxx (1)
c---- Coefficient of term y**2.
dimension qyy (1)
c---- Coefficient of term z**2.
dimension qzz (1)
c---- Coefficient of term 1.
dimension sc (1)
c---- Coefficient of term x.
dimension sx (1)
c---- Coefficient of term y.
dimension sy (1)
c---- Coefficient of term z.
dimension sz (1)
c---- Coefficient of term x * y.
dimension sxy (1)
c---- Coefficient of term y * z.
dimension syz (1)
c---- Coefficient of term z * x.
dimension szx (1)
c---- Coefficient of term x**2.
dimension sxx (1)
c---- Coefficient of term y**2.
dimension syy (1)
c---- Coefficient of term z**2.
dimension szz (1)
c.... Local variables.
c---- Index in an external array.
common /laptplis/ n
c---- Estimated error in sc.
common /laptplis/ scerr
c---- Estimated error in sx.
common /laptplis/ sxerr
c---- Estimated error in sy.
common /laptplis/ syerr
c---- Estimated error in sz.
common /laptplis/ szerr
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptplis finding the intersection of a major plane',
cbug & ' with a quadric surface.' /
cbug & ' tol=',1pe13.5 /
cbug & ' kpl=',i3,' apl=',1pe22.14)
cbug 9902 format (
cbug & i3,' 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, kpl, apl
cbug write ( 3, 9902) (n, qc(n), qx(n), qy(n), qz(n),
cbug & qxy(n), qyz(n), qzx(n), qxx(n), qyy(n), qzz(n), n = 1, np)
cbugc***DEBUG ends.
c.... Initialize.
nerr = 0
c.... Test for input errors.
if (np .le. 0) then
nerr = 1
cbugc***DEBUG begins.
cbug 9903 format (/ " aptplis FATAL. Bad np=",i3)
cbug write ( 3, 9903) np
cbugc***DEBUG ends.
go to 210
endif
if ((kpl .lt. 1) .or. (kpl .gt. 3)) then
nerr = 2
cbugc***DEBUG begins.
cbug 9904 format (/ " aptplis FATAL. Bad kpl=",i3)
cbug write ( 3, 9904) kpl
cbugc***DEBUG ends.
go to 210
c---- Tested kpl.
endif
c.... Find the orientation of the major plane, then find intersections.
c---- Plane at x = apl.
if (kpl .eq. 1) then
c---- Loop over external arrays.
do 120 n = 1, np
sc(n) = qc(n) + qx(n) * apl + qxx(n) * apl**2
sx(n) = 0.0
sy(n) = qy(n) + qxy(n) * apl
sz(n) = qz(n) + qzx(n) * apl
sxy(n) = 0.0
syz(n) = qyz(n)
szx(n) = 0.0
sxx(n) = 0.0
syy(n) = qyy(n)
szz(n) = qzz(n)
c---- End of loop over external arrays.
120 continue
c---- Test for truncation error.
tolx = tol
if (tol .eq. 0.0) tolx = 1.e-11
c---- Loop over external arrays.
do 130 n = 1, np
scerr = tolx * (abs (qc(n)) + abs (qx(n) * apl) +
& abs (qxx(n)) * apl**2)
syerr = tolx * (abs (qy(n)) + abs (qxy(n) * apl))
szerr = tolx * (abs (qz(n)) + abs (qzx(n) * apl))
if (abs (sc(n)) .lt. scerr) then
sc(n) = 0.0
endif
if (abs (sy(n)) .lt. syerr) then
sy(n) = 0.0
endif
if (abs (sz(n)) .lt. szerr) then
sz(n) = 0.0
endif
c---- End of loop over external arrays.
130 continue
c---- Plane at y = apl.
elseif (kpl .eq. 2) then
c---- Loop over external arrays.
do 140 n = 1, np
sc(n) = qc(n) + qy(n) * apl + qyy(n) * apl**2
sx(n) = qx(n) + qxy(n) * apl
sy(n) = 0.0
sz(n) = qz(n) + qyz(n) * apl
sxy(n) = 0.0
syz(n) = 0.0
szx(n) = qzx(n)
sxx(n) = qxx(n)
syy(n) = 0.0
szz(n) = qzz(n)
c---- End of loop over external arrays.
140 continue
c---- Test for truncation error.
c---- Loop over external arrays.
do 150 n = 1, np
scerr = tolx * (abs (qc(n)) + abs (qy(n) * apl) +
& abs (qyy(n)) * apl**2)
sxerr = tolx * (abs (qx(n)) + abs (qxy(n) * apl))
szerr = tolx * (abs (qz(n)) + abs (qyz(n) * apl))
if (abs (sc(n)) .lt. scerr) then
sc(n) = 0.0
endif
if (abs (sx(n)) .lt. sxerr) then
sx(n) = 0.0
endif
if (abs (sz(n)) .lt. szerr) then
sz(n) = 0.0
endif
c---- End of loop over external arrays.
150 continue
c---- Plane at z = apl.
else
c---- Loop over external arrays.
do 160 n = 1, np
sc(n) = qc(n) + qz(n) * apl + qzz(n) * apl**2
sx(n) = qx(n) + qzx(n) * apl
sy(n) = qy(n) + qyz(n) * apl
sz(n) = 0.0
sxy(n) = qxy(n)
syz(n) = 0.0
szx(n) = 0.0
sxx(n) = qxx(n)
syy(n) = qyy(n)
szz(n) = 0.0
c---- End of loop over external arrays.
160 continue
c---- Test for truncation error.
c---- Loop over external arrays.
do 170 n = 1, np
scerr = tolx * (abs (qc(n)) + abs (qz(n) * apl) +
& abs (qzz(n)) * apl**2)
sxerr = tolx * (abs (qx(n)) + abs (qzx(n) * apl))
syerr = tolx * (abs (qy(n)) + abs (qyz(n) * apl))
if (abs (sc(n)) .lt. scerr) then
sc(n) = 0.0
endif
if (abs (sx(n)) .lt. sxerr) then
sx(n) = 0.0
endif
if (abs (sy(n)) .lt. syerr) then
sy(n) = 0.0
endif
c---- End of loop over external arrays.
170 continue
c---- Tested kpl.
endif
cbugc***DEBUG begins.
cbug 9905 format (/ 'aptplis results:' /
cbug & i3,' sc= ',1pe22.14 /
cbug & ' sx,sy,sz= ',1p3e22.14 /
cbug & ' sxy,syz,szx=',1p3e22.14 /
cbug & ' sxx,syy,szz=',1p3e22.14)
cbug write ( 3, 9905) (n, sc(n), sx(n), sy(n), sz(n),
cbug & sxy(n), syz(n), szx(n), sxx(n), syy(n), szz(n), n = 1, np)
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptplis. (+1 line.)
end
UCRL-WEB-209832