subroutine aptparl (px, py, pz, vx, vy, vz, ax, ay, az,
& bx, by, bz, cx, cy, cz, tol,
& nints, atype, ptx, pty, ptz,
& vtx, vty, vtz, vtlen,
& dprox, dline, tprox, pltx, plty, pltz, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTPARL
c
c subroutine aptparl (px, py, pz, vx, vy, vz, ax, ay, az,
c bx, by, bz, cx, cy, cz, tol,
c nints, atype, ptx, pty, ptz,
c vtx, vty, vtz, vtlen,
c dprox, dline, tprox, pltx, plty, pltz, nerr)
c
c Version: aptparl Updated 1997 April 21 15:00.
c aptparl Originated 1997 April 17 18:20.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: For the particle with initial coordinates p = (px, py, pz),
c initial velocity v = (vx, vy, vz) and constant acceleration
c a = (ax, ay, az), and for the extended line through points
c b = (bx, by, bz) and c = (cx, cy, cz), to find the number nints
c (1, 2 or 3) of proximal or intersection points
c pt(n) = (ptx(n), pty(n), ptz(n)), n = 1, nints, of the particle
c path with the line, and for each such point pt(n), to find the
c velocities vtlen(n) = (vtx(n), vty(n), vtz(n)), the speeds
c (velocity magnitudes) vtlen(n), the path lengths from point p,
c dprox(n), the distances dline(n) from points pt(n) to the line,
c the times from point p to points pt(n), tprox(n), and the
c corresponding proximal points plt = (pltx, plty, pltz) on
c the line.
c
c The particle velocity and coordinates are given by:
c vt = v + a*t,
c pt = p + v*t + 0.5*a*t**2.
c The line direction is the unit vector u = (c - b) / |c - b|.
c
c The distance vector from the particle to the line is given by:
c d = (pt - b) - ((pt - b) dot u) * u
c where is the unit vector of the line, u = (c - b) / |c - b|.
c
c To find the times of minimum distances, differentiate d**2
c with respect to time, and solve the resulting cubic equation for
c t = tprox, to find either one or three real roots.
c c0 + c1 * t * c2 * t**2 + c3 * t**3 = 0,
c where
c c0 = 2 * bpn dot vn
c c1 = 2 * (bpn dot an + vn**2)
c c2 = 3 * vn dot an
c c4 = an**2
c where
c bpl = ((p - b) dot u) * u
c vl = (v dot u) * u
c al = (a dot u) * u
c bpn = p - b - bpl
c vn = v - vl
c an = a - al
c
c The proximal or intersection points pt, and the velocity vt are
c found by substituting the roots tprox in the equations for pt
c and vt. The proximal distance dline is the magnitude of d.
c The proximal points on the line are at t = tprox:
c plt = pt - d = b + bpl + vl * t + 0.5 * al * t**2
c
c See subroutine aptparb for finding dprox.
c
c History:
c
c Input: px, py, pz, vx, vy, vz, ax, ay, az, bx, by, bz, cx, cy, cz, tol.
c
c Output: nints, atype, ptx, pty, ptz, vtx, vty, vtz, vtlen, dprox, dline,
c tprox, pltx, plty, pltz, nerr.
c
c Calls: aptcubs, aptparb, aptvdis, aptvdot, aptvpap, aptvsum, aptvunb
c
c
c Glossary:
c
c ax,ay,az Input The x, y and z components of the constant acceleration
c vector "a".
c
c atype Output The type of result, "proximal", "on line", or
c "unknown", the latter if the solution to the cubic
c equation for time is not accurate.
c
c bx,by,bz Input The x, y and z coordinates of point "b" on the line.
c
c cx,cy,cz Input The x, y and z coordinates of point "c" on the line.
c
c dline Output The distances from the proximal points or intersections
c on the partice path to the line. Size 3.
c
c dprox Output The distance along the particle path to the proximal
c point. Size = 3.
c
c nerr Output Indicates an input error, if nonzero.
c 1 if the line has zero length (points b and c are
c the same).
c
c nints Output The number of proximal points or intersections of the
c particle path and the line.
c
c pltx,y,z Output The x, y and z coordinates of the proximal points or
c intersections on the line. Size 3.
c
c ptx,y,z Output The x, y and z coordinates of the proximal points or
c intersections of the particle path with the line.
c Size 3.
c
c px,py,pz Input The initial x, y and z coordinates of the particle.
c
c tol Input Numerical tolerance limit.
c On Cray computers, recommend 1.e-5 to 1.e-11.
c
c tprox Output The times along the particle path to the proximal
c points or intersections with the line. Size 3.
c
c vtlen Output The magnitudes of velocities vt. Size 3.
c
c vtx,y,z Output The x, y and z components of the particle velocity
c at points pt. Size = 3.
c
c vx,vy,vz Input The initial x, y and z components of the particle
c velocity.
ccend.
c.... Dimensioned arguments.
dimension atype(1) ! Type of point, "proximal" or "on point".
character*8 atype ! Type of point, "proximal" or "on point".
dimension dprox(1) ! Path distance to prox or intersection.
dimension dline(1) ! Path distance to prox or intersection.
dimension pltx(1) ! Coordinate x of prox point on line.
dimension plty(1) ! Coordinate y of prox point on line.
dimension pltz(1) ! Coordinate z of prox point on line.
dimension ptx(1) ! Coordinate x of prox point on path.
dimension pty(1) ! Coordinate y of prox point on path.
dimension ptz(1) ! Coordinate z of prox point on path.
dimension tprox(1) ! Time to proximity or intersection.
dimension vtlen(1) ! Magnitude of particle velocity.
dimension vtx(1) ! Component x of particle velocity.
dimension vty(1) ! Component y of particle velocity.
dimension vtz(1) ! Component z of particle velocity.
c.... Local variables.
dimension timag(3) ! Imaginary part of complex time.
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptparl finding proximal points or intersections' /
cbug & ' of particle path and a line. tol= ',1pe22.14)
cbug 9902 format (' Particle at p with velocity v, acceleration a:' /
cbug & ' px,py,pz= ',1p3e22.14 /
cbug & ' vx,vy,vz= ',1p3e22.14 /
cbug & ' ax,ay,az= ',1p3e22.14 )
cbug 9903 format (' Points b and c:' /
cbug & ' bx,by,bz= ',1p3e22.14 /
cbug & ' cx,cy,cz= ',1p3e22.14 )
cbug write ( 3, 9901) tol
cbug write ( 3, 9902) px, py, pz, vx, vy, vz, ax, ay, az
cbug write ( 3, 9903) bx, by, bz, cx, cy, cz
cbugc***DEBUG ends.
c.... Initialize.
nerr = 0
nints = 0
do 100 n = 1, 3
dprox(n) = -1.e99
pltx(n) = -1.e99
plty(n) = -1.e99
pltz(n) = -1.e99
ptx(n) = -1.e99
pty(n) = -1.e99
ptz(n) = -1.e99
tprox(n) = -1.e99
vtlen(n) = -1.e99
vtx(n) = -1.e99
vty(n) = -1.e99
vtz(n) = -1.e99
100 continue
c.... Test for input errors.
c.... Find the unit vector ubc from point "b" to point "c".
call aptvdis (bx, by, bz, cx, cy, cz, 1, tol,
& bcx, bcy, bcz, bclen, nerr)
if (bclen .eq. 0.0) then ! Line has zero length, no direction.
nerr = 1
go to 999
endif
call aptvunb (bcx, bcy, bcz, 1, tol, ux, uy, uz, ulen, nerr)
c.... Find the vector "bp" from point "b" to point "p".
call aptvdis (bx, by, bz, px, py, pz, 1, tol,
& bpx, bpy, bpz, bplen, nerr)
c.... Find the components of vector "bp" parallel and perpendicular to line.
call aptvpap (bpx, bpy, bpz, bcx, bcy, bcz, 1, tol,
& bplx, bply, bplz, bpllen,
& bpnx, bpny, bpnz, bpnlen, nerr)
c.... Find the point "pl" on the line proximal to point "p".
call aptvsum (0, 1., bx, by, bz, 1., bplx, bply, bplz, 1, tol,
& plx, ply, plz, pllen, nerr)
c.... Find the components of vector "v" parallel and perpendicular to line.
call aptvpap (vx, vy, vz, bcx, bcy, bcz, 1, tol,
& vlx, vly, vlz, vllen,
& vnx, vny, vnz, vnlen, nerr)
c.... Find the components of vector "a" parallel and perpendicular to line.
call aptvpap (ax, ay, az, bcx, bcy, bcz, 1, tol,
& alx, aly, alz, allen,
& anx, any, anz, anlen, nerr)
c.... Find the dot product of bpn and vn.
call aptvdot (bpnx, bpny, bpnz, vnx, vny, vnz, 1, tol,
& bpdotv, nerr)
c.... Find the dot product of bpn and an.
call aptvdot (bpnx, bpny, bpnz, anx, any, anz, 1, tol,
& bpdota, nerr)
c.... Find the dot product of vn and an.
call aptvdot (vnx, vny, vnz, anx, any, anz, 1, tol,
& vdota, nerr)
c.... Find the times of proximity or intersection.
if ((vnlen .eq. 0.0) .and. (anlen .eq. 0.0)) then
nints = 1
tprox(1) = 0.0
elseif (anlen .eq. 0.0) then ! Path is linear.
nints = 1
tprox(1) = -bpdotv / vnlen**2
else ! Path is parabolic.
c.... Find the coefficients of the cubic equation for time.
c0 = 2.0 * bpdotv / anlen**2
c1 = 2.0 * (bpdota + vnlen**2) / anlen**2
c2 = 3.0 * vdota / anlen**2
tolx = tol
if (tol .eq. 0.0) tolx = 1.e-11
call aptcubs (c0, c1, c2, tolx,
& nints, tprox, timag, atype, itrun)
endif ! Tested type of path.
c.... Find the proximal or intersection points on the path and the line.
do 140 n = 1, nints
call aptparb (px, py, pz, vx, vy, vz, ax, ay, az,
& tprox(n), tol,
& ptx(n), pty(n), ptz(n),
& vtx(n), vty(n), vtz(n), vtlen(n),
& dprox(n), ntype)
call aptparb (plx, ply, plz, vlx, vly, vlz, alx, aly, alz,
& tprox(n), tol,
& pltx(n), plty(n), pltz(n),
& vltx, vlty, vltz, vltlen,
& dlprox, ntype)
call aptvdis (ptx(n), pty(n), ptz(n),
& pltx(n), plty(n), pltz(n), 1, tol,
& dx, dy, dz, dline(n), nerr)
if (itrun .eq. 2) then
atype(n) = 'unknown'
elseif (dline(n) .eq. 0.0) then
atype(n) = 'ON LINE'
else
atype(n) = 'proximal'
endif
140 continue
999 continue
cbugc***DEBUG begins.
cbug 9904 format (/ 'aptparl results: nints=',i2 )
cbug 9905 format (/ 'Path point type ',a8 /
cbug & 'Proximal dist ',1pe22.14,2x,a8)
cbug 9906 format ( 'Time, path len',1p2e22.14)
cbug 9907 format ( 'Prox pt(xyz) ',1p3e22.14)
cbug 9908 format ( 'Prox vel(xyz) ',1p3e22.14 /
cbug & ' speed ',1pe22.14 )
cbug 9909 format ( 'Line pt(xyz) ',1p3e22.14)
cbug write ( 3, 9904) nints
cbug do 710 n = 1, nints ! Loop over proximal points.
cbug write ( 3, 9905) atype(n), dline(n)
cbug write ( 3, 9906) tprox(n), dprox(n)
cbug write ( 3, 9907) ptx(n), pty(n), ptz(n)
cbug write ( 3, 9908) vtx(n), vty(n), vtz(n), vtlen(n)
cbug write ( 3, 9909) pltx(n), plty(n), pltz(n)
cbug 710 continue ! End of loop over proximal points.
cbugc***DEBUG ends.
return
c.... End of subroutine aptparl. (+1 line.)
end
UCRL-WEB-209832