subroutine aptrkis (ax, ay, az, bx, by, bz, qc, qx, qy, qz,
& qxy, qyz, qzx, qxx, qyy, qzz, dintmn, dintmx,
& np, tol, nints, cx, cy, cz, dint,
& sx, sy, sz, costh, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTRKIS
c
c call aptrkis (ax, ay, az, bx, by, bz, qc, qx, qy, qz,
c & qxy, qyz, qzx, qxx, qyy, qzz, dintmn, dintmx,
c & np, tol, nints, cx, cy, cz, dint,
c & sx, sy, sz, costh, nerr)
c
c Version: aptrkis Updated 1994 October 13 11:00.
c aptrkis Originated 1990 February 23 10:00.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To find, for each of np sets of input data, the distance dint
c to the intersection c = (cx, cy, cz) of the linear track through
c point a = (ax, ay, az) with direction vector b = (bx, by, bz),
c and the general second order surface for which the equation is
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 and for which dint is between dintmn and dintmx. If two such
c intersections occur, the one with the smaller magnitude of dint
c will be returned. If no such intersection is found, nints will
c be 0, and dint and the coordinates of point "c" will be very
c large. Also, to find the vector "s" normal to the surface at
c point "c", and the angle costh between vector "b" and vector
c "s". Flag nerr indicates any input error. If dint is less
c than tol times the distance from the origin to point "a", dint
c will be truncated to zero, and "c" will equal "a".
c
c The vector normal to the surface is s = (df/dx, df/dy, df/dz).
c The sign of the direction of the intersection is determined by
c the dot product b * s.
c
c History: 1990 March 14. Changed tol to 0.0 in call to unit vector
c subroutine. Allows small magnitudes.
c 1994 February 15. Truncated dint to 0.0 if small compared with
c the distance from the origin to point "a".
c 1994 October 13. Added arguments sx, sy, sz and costh.
c
c Input: ax, ay, az, bx, by, bz, qc, qx, qy, qz,
c qxy, qyz, qzx, qxx, qyy, qzz, dintmn, dintmx, np, tol.
c
c Output: nints, cx, cy, cz, dint, sx, sy, sz, costh, nerr.
c
c Calls: aptqrtv, aptvadd, aptvunb
c
c
c Glossary:
c
c ax,ay,az Input The x, y, z coordinates of point "a". Size np.
c
c bx,by,bz Input The x, y, z components of direction vector "b".
c If the magnitude is too small, nints will be -1.
c
c costh Output The cosine of the angle between direction vector "b"
c and the normal vector "s" at the intersection point
c "c". Size np.
c
c cx,cy,cz Output The x, y, z coordinates of the point of intersection
c of the line through point "a" with unit direction
c vector "b", and the surface, if nints = 1. Size np.
c
c dint Output The distance of the point of intersection "c" from
c point "a", if nints = 1. Size np.
c Positive if in the same direction as vector "b".
c Acceptable only if between dintmn and dintmx.
c
c dintmn Input The minimum allowable value of dint. Size np.
c
c dintmx Input The maximum allowable value of dint. Size np.
c
c nerr Output Indicates an input error, if not 0.
c 1 if np is not positive.
c
c nints Output Indicates the type of intersection found:
c 1 if the track through point "a" in direction "b"
c intersects the surface at a distance dint between
c dintmn and dintmx.
c 0 if no such intersection was found.
c -1 if vector "b" is too short, based on tol.
c Size np.
c
c np Input Size of arrays.
c
c q.. Input Coefficients of the implicit equation of a second-order
c surface in xyz space (qc, qx, qy, qz, qxy, qyz, qzx,
c qxx, qyy, qzz). Size np.
c
c sx,sy,sz Output The x, y z components of the normal vector at the
c intersection point "c" on the quadric surface.
c
c tol Input Numerical tolerance limit.
c On Cray computers, recommend 1.e-5 to 1.e-11.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c---- Coordinate x of point "a".
dimension ax (1)
c---- Coordinate y of point "a".
dimension ay (1)
c---- Coordinate z of point "a".
dimension az (1)
c---- Cosine of angle between "b" and "s".
dimension costh (1)
c---- Coordinate x of point "c".
dimension cx (1)
c---- Coordinate y of point "c".
dimension cy (1)
c---- Coordinate z of point "c".
dimension cz (1)
c---- Component x of vector "b".
dimension bx (1)
c---- Component y of vector "b".
dimension by (1)
c---- Component z of vector "b".
dimension bz (1)
c---- Distance from point "a" to point "c".
dimension dint (1)
c---- Minimum allowed value of dint.
dimension dintmn (1)
c---- Maximum allowed value of dint.
dimension dintmx (1)
c---- Track intersects surface, if 1.
dimension nints (1)
c---- Coefficient of term 1.
dimension qc (1)
c---- Coefficient of term x.
dimension qx (1)
c---- Coefficient of term x**2.
dimension qxx (1)
c---- Coefficient of term xy.
dimension qxy (1)
c---- Coefficient of term y.
dimension qy (1)
c---- Coefficient of term y**2.
dimension qyy (1)
c---- Coefficient of term yz.
dimension qyz (1)
c---- Coefficient of term z.
dimension qz (1)
c---- Coefficient of term zx.
dimension qzx (1)
c---- Coefficient of term z**2.
dimension qzz (1)
c---- Component x of vector "s".
dimension sx (1)
c---- Component y of vector "s".
dimension sy (1)
c---- Component z of vector "s".
dimension sz (1)
c.... Local variables.
c---- Coefficient of dint**2 in quadratic.
common /laptrkis/ aa (64)
c---- Coefficient of dint in quadratic.
common /laptrkis/ bb (64)
c---- A very big number.
common /laptrkis/ big
c---- Magnitude of vector "b".
common /laptrkis/ blen (64)
c---- Constant term in quadratic.
common /laptrkis/ cc (64)
c---- First real root, if nroots .gt. 0.
common /laptrkis/ dint1 (64)
c---- Second real root, if nroots = 2.
common /laptrkis/ dint2 (64)
c---- Truncation error in qq, if 1.
common /laptrkis/ itrun (64)
c---- Index in arrays.
common /laptrkis/ n
c---- First index of subset of data.
common /laptrkis/ n1
c---- Last index of subset of data.
common /laptrkis/ n2
c---- 1 if dint1 acceptable.
common /laptrkis/ nexit1
c---- 1 if dint2 acceptable.
common /laptrkis/ nexit2
c---- Index in external array.
common /laptrkis/ nn
c---- Number of real roots of quadratic.
common /laptrkis/ nroots (64)
c---- Size of current subset of data.
common /laptrkis/ ns
c---- Value of bb**2 - 4.0 * aa * cc.
common /laptrkis/ qq (64)
c---- Component x of unit vector "u".
common /laptrkis/ ux (64)
c---- Component y of unit vector "u".
common /laptrkis/ uy (64)
c---- Component z of unit vector "u".
common /laptrkis/ uz (64)
c---- Magnitude of a vector.
common /laptrkis/ vlen (64)
cbugc***DEBUG begins.
cbug common /laptrkis/ ddistx (64)
cbug common /laptrkis/ ddisty (64)
cbug common /laptrkis/ ddistz (64)
cbug common /laptrkis/ ddist (64)
cbug 9901 format (/ 'aptrkis finding intersection of the line through' /
cbug & ' point a with direction vector b, with the second-order' /
cbug & ' surface with the coefficients of the given terms:' /
cbug & (i3,' ax,ay,az=',1p3e22.14 /
cbug & ' bx,by,bz=',1p3e22.14 /
cbug & ' 1 =',1pe22.14 /
cbug & ' x, y, z =',1p3e22.14 /
cbug & " xy,yz,zx=",1p3e22.14 /
cbug & " xx,yy,zz=",1p3e22.14 /
cbug & " dintmn,x=",1p2e22.14))
cbug write ( 3, 9901) (n, ax(n), ay(n), az(n), bx(n), by(n), bz(n),
cbug & qc(n), qx(n), qy(n), qz(n), qxy(n), qyz(n), qzx(n),
cbug & qxx(n), qyy(n), qzz(n), dintmn(n), dintmx(n), n = 1, np)
cbugc***DEBUG ends.
c.... Initialize.
c---- A very big number.
big = 1.e+99
nerr = 0
c.... Test for input errors.
if (np .le. 0) then
nerr = 1
go to 210
endif
c.... Set up the indices of the first subset of data.
n1 = 1
n2 = min (np, 64)
110 ns = n2 - n1 + 1
c.... Find the unit vector parallel to direction vector "b".
call aptvunb (bx(n1), by(n1), bz(n1), ns, 0.,
& ux, uy, uz, blen, nerr)
c.... Find the coefficients of the equation for dint.
c---- Loop over subset of data.
do 120 n = 1, ns
nn = n + n1 - 1
aa(n) = ux(n) * (qxx(nn) * ux(n) + qxy(nn) * uy(n)) +
& uy(n) * (qyy(nn) * uy(n) + qyz(nn) * uz(n)) +
& uz(n) * (qzz(nn) * uz(n) + qzx(nn) * ux(n))
bb(n) = ux(n) * (qx(nn) + 2.0 * qxx(nn) * ax(nn) +
& qxy(nn) * ay(nn) + qzx(nn) * az(nn)) +
& uy(n) * (qy(nn) + 2.0 * qyy(nn) * ay(nn) +
& qyz(nn) * az(nn) + qxy(nn) * ax(nn)) +
& uz(n) * (qz(nn) + 2.0 * qzz(nn) * az(nn) +
& qzx(nn) * ax(nn) + qyz(nn) * ay(nn))
cc(n) = qc(nn) +
& ax(nn) * (qx(nn) + qxx(nn) * ax(nn) + qxy(nn) * ay(nn)) +
& ay(nn) * (qy(nn) + qyy(nn) * ay(nn) + qyz(nn) * az(nn)) +
& az(nn) * (qz(nn) + qzz(nn) * az(nn) + qzx(nn) * ax(nn))
c---- End of loop over subset of data.
120 continue
c.... Find the roots of the equation for dint.
call aptqrtv (0, aa, bb, cc, qq, ns, tol,
& nroots, dint1, dint2, itrun, nerr)
c.... Find the nearest value of dint between dintmn and dintmx, if any.
c---- Loop over local arrays.
do 130 n = 1, ns
nn = n + n1 - 1
if ((nroots(n) .ge. 1) .and.
& (dint1(n) .ge. dintmn(nn)) .and.
& (dint1(n) .le. dintmx(nn)) ) then
nexit1 = 1
else
nexit1 = 0
endif
if (nexit1 .ne. 1) then
dint1(n) = -big
endif
if ((nroots(n) .eq. 2) .and.
& (abs (dint2(n)) .lt.
& abs (dint1(n))) .and.
& (dint2(n) .ge. dintmn(nn)) .and.
& (dint2(n) .le. dintmx(nn)) ) then
nexit2 = 1
else
nexit2 = 0
endif
if (nexit1 .eq. 1) then
nints(nn) = 1
else
nints(nn) = 0
endif
if (nexit2 .eq. 1) then
nints(nn) = 1
endif
if (blen((n)) .le. tol) then
nints(nn) = -1
endif
if (nexit1 .eq. 1) then
dint(nn) = dint1(n)
else
dint(nn) = -big
endif
if (nexit2 .eq. 1) then
dint(nn) = dint2(n)
endif
dinterr2 = tol**2 * (ax(nn)**2 + ay(nn)**2 + az(nn)**2)
if (dint(nn)**2 .lt. dinterr2) dint(nn) = 0.0
c---- End of loop over local arrays.
130 continue
c.... Find the intersection points "c".
call aptvadd (ax(n1), ay(n1), az(n1), 1.0, dint(n1),
& ux(n1), uy(n1), uz(n1), ns, tol,
& cx(n1), cy(n1), cz(n1), vlen, nerr)
c.... Find the normal vectors "s".
do 140 n = 1, ns
nn = n + n1 - 1
sx(nn) = qx(nn) + 2.0 * qxx(nn) * cx(nn) + qxy(nn) * cy(nn) +
& qzx(nn) * cz(nn)
sy(nn) = qy(nn) + qxy(nn) * cx(nn) + 2.0 * qyy(nn) * cy(nn) +
& qyz(nn) * cz(nn)
sz(nn) = qz(nn) + qzx(nn) * cx(nn) + qyz(nn) * cy(nn) +
& 2.0 * qzz(nn) * cz(nn)
140 continue
c.... Find the cosines of the angles between "b" and "s".
call aptvang (bx(n1), by(n1), bz(n1), sx(n1), sy(n1), sz(n1),
& 1, tol, costh(n1), nerr)
cbugc***DEBUG begins.
cbug call aptvdis (ax(n1), ay(n1), az(n1), cx(n1), cy(n1), cz(n1),
cbug & ns, tol, dbugx, dbugy, dbugz, ddist, nerra)
cbugc***DEBUG ends.
c.... See if all data subsets are done.
c---- Do another subset of data.
if (n2 .lt. np) then
n1 = n2 + 1
n2 = min (np, n1 + 63)
go to 110
endif
210 continue
cbugc***DEBUG begins.
cbug 9902 format (/ 'aptrkis results: nerr=',i3 /
cbug & (i3,' dint= ',1pe22.14,' nints=',i2 /
cbug & ' cx,cy,cz=',1p3e22.14 /
cbug & ' sx,sy,sz=',1p3e22.14 /
cbug & ' costh= ',1pe22.14 ))
cbug 9903 format (i3,' ddist= ',1pe22.14)
cbug write ( 3, 9902) nerr, (n, dint(n), nints(n),
cbug & cx(n), cy(n), cz(n), sx(n), sy(n), sz(n), costh(n), n = 1, np)
cbug npx = np
cbug if (npx .gt. 64) npx = 64
cbug call aptvdis (ax, ay, az, cx, cy, cz,
cbug & npx, tol, ddistx, ddisty, ddistz, ddist, nerra)
cbug write ( 3, 9903) (n, ddist(n), n = 1, npx)
cbugc***DEBUG ends.
return
c.... End of subroutine aptrkis. (+1 line.)
end
UCRL-WEB-209832