subroutine aptpars (px, py, pz, vx, vy, vz, ax, ay, az,
& bx, by, bz, tol,
& nints, atype, ptx, pty, ptz,
& vtx, vty, vtz, vtlen, dint, dprox, tint)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTPARS
c
c call aptpars (px, py, pz, vx, vy, vz, ax, ay, az,
c & bx, by, bz, tol,
c & nints, atype, ptx, pty, ptz,
c & vtx, vty, vtz, vtlen, dint, dprox, tint)
c
c Version: aptpars Updated 1997 April 16 16:20.
c aptpars Originated 1997 April 16 16: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 point b = (bx, by, bz), to find
c the number nints (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 point, and for each such point pt(n), to find the
c velocity vtlen(n) = (vtx(n), vty(n), vtz(n)), the speed
c (velocity magnitude) vtlen(n), the path length from point p,
c dint(n), the distance dprox(n) from point pt(n) to point b,
c and the time tint(n) from point p to point pt(n).
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
c The distance from the particle to point b is given by:
c dp = |pt - b| = |(p - b) + v*t + 0.5*a*t**2|
c To find the minimum distances, differentiate dp**2 with respect
c to time, and solve the resulting cubic equation for t = tint,
c with either one or three real roots.
c c0 + c1 * t * c2 * t**2 + c3 * t**3 = 0, where
c c0 = 2 * (p - b) dot v
c c1 = 2 * ((p - b) dot a + v**2)
c c2 = 3 * v dot a
c c4 = a**2
c If a = 0, the particle path is linear, and one root occurs:
c t = -(p - b) dot v / a**2.
c The proximal or intersection points pt, and the velocity are
c found by substituting the roots tint in the equations for pt
c and vt. The proximal distance is the magnitude of pt - b.
c
c See subroutine aptparb for finding dint.
c
c History:
c
c Input: px, py, pz, vx, vy, vz, ax, ay, az, bx, by, bz, tol.
c
c Output: nints, atype, ptx, pty, ptz, vtx, vty, vtz, vtlen, dint, dprox,
c tint.
c
c Calls: aptcubs, aptparb, aptvdis, aptvdot, 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 point", 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".
c
c dint Output The distance along the particle path to the proximal
c point. Size = 3.
c
c dprox Output The distances from the proximal points or intersections
c on the partice path to point "b". Size = 3.
c
c nints Output The number of proximal points or intersections of the
c particle path and point "b".
c
c ptx,y,z Output The x, y and z coordinates of the proximal points or
c intersections of the particle path with point "b".
c Size = 3.
c
c pz,py,pz Input The initial x, y and z coordinates of the particle.
c
c tint Output The times along the particle path to the proximal
c points or intersections with point "b". Size = 3.
c
c tol Input Numerical tolerance limit.
c On Cray computers, recommend 1.e-5 to 1.e-11.
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(3) ! Type of point, "proximal" or "on point".
character*8 atype ! Type of point, "proximal" or "on point".
dimension dint(3) ! Path distance to prox or intersection.
dimension dprox(3) ! Distance from prox point to point "b".
dimension ptx(3) ! Coordinate x of prox or int point.
dimension pty(3) ! Coordinate y of prox or int point.
dimension ptz(3) ! Coordinate z of prox or int point.
dimension tint(3) ! Time to proximity or intersection.
dimension vtlen(3) ! Magnitude of particle velocity.
dimension vtx(3) ! Component x of particle velocity.
dimension vty(3) ! Component y of particle velocity.
dimension vtz(3) ! Component z of particle velocity.
c.... Local variables.
dimension timag(3) ! Imaginary part of complex time.
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptpars finding proximal points or interscctions' /
cbug & ' of particle path and a point. 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 (' Point b:' /
cbug & ' bx,by,bz= ',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
cbugc***DEBUG ends.
c.... Initialize.
nints = 0
do 100 n = 1, 3
dint(n) = -1.e99
ptx(n) = -1.e99
pty(n) = -1.e99
ptz(n) = -1.e99
tint(n) = -1.e99
vtlen(n) = -1.e99
vtx(n) = -1.e99
vty(n) = -1.e99
vtz(n) = -1.e99
100 continue
c.... Find the magnitude of the initial velocity vector "v".
call aptvunb (vx, vy, vz, 1, tol, uvx, uvy, uvz, vlen, nerr)
c.... Find the magnitude of the acceleration vector "a".
call aptvunb (ax, ay, az, 1, tol, uax, uay, uaz, alen, 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 dot product of vector "bp" and the initial velocity "v".
call aptvdot (bpx, bpy, bpz, vx, vy, vz, 1, tol, bpdotv, nerr)
c.... Find the dot product of vector "bp" and the acceleration vector "a".
call aptvdot (bpx, bpy, bpz, ax, ay, az, 1, tol, bpdota, nerr)
c.... Find the dot product of vector "v" and the acceleration vector "a".
call aptvdot (vx, vy, vz, ax, ay, az, 1, tol, vdota, nerr)
c.... Find the times of proximity or intersection.
if ((vlen .eq. 0.0) .and. (alen .eq. 0.0)) then ! Particle is fixed.
nints = 1
tint(1) = 0.0
elseif (alen .eq. 0.0) then ! Path is linear.
nints = 1
tint(1) = -bpdotv / vlen**2
else ! Path is parabolic.
c.... Find the coefficients of the cubic equation for time.
c0 = 2.0 * bpdotv / alen**2
c1 = 2.0 * (bpdota + vlen**2) / alen**2
c2 = 3.0 * vdota / alen**2
tolx = tol
if (tol .eq. 0.0) tolx = 1.e-11
call aptcubs (c0, c1, c2, tolx,
& nints, tint, timag, atype, itrun)
endif ! Tested type of path.
c.... Find the proximal or intersection points on the path.
do 140 n = 1, nints
call aptparb (px, py, pz, vx, vy, vz, ax, ay, az,
& tint(n), tol,
& ptx(n), pty(n), ptz(n),
& vtx(n), vty(n), vtz(n), vtlen(n),
& dint(n), ntype)
call aptvdis (bx, by, bz, ptx(n), pty(n), ptz(n), 1, tol,
& bpx, bpy, bpz, dprox(n), nerr)
if (itrun .eq. 2) then
atype(n) = 'unknown'
elseif (dprox(n) .eq. 0.0) then
atype(n) = 'ON POINT'
else
atype(n) = 'proximal'
endif
140 continue
999 continue
cbugc***DEBUG begins.
cbug 9904 format (/ 'aptpars 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 write ( 3, 9904) nints
cbug do 710 n = 1, nints ! Loop over proximal points.
cbug write ( 3, 9905) atype(n), dprox(n)
cbug write ( 3, 9906) tint(n), dint(n)
cbug write ( 3, 9907) ptx(n), pty(n), ptz(n)
cbug write ( 3, 9908) vtx(n), vty(n), vtz(n), vtlen(n)
cbug 710 continue ! End of loop over proximal points.
cbugc***DEBUG ends.
return
c.... End of subroutine aptpars. (+1 line.)
end
UCRL-WEB-209832