subroutine aptparh (px, py, pz, vlen, ax, ay, az, bx, by, bz, tol,
& nints, vx, vy, vz, vtx, vty, vtz, vtlen,
& dint, tint)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTPARH
c
c call aptparh (px, py, pz, vlen, ax, ay, az, bx, by, bz, tol,
c & nints, vx, vy, vz, vtx, vty, vtz, vtlen,
c & dint, tint)
c
c Version: aptparh Updated 1997 April 21 17:30.
c aptparh Originated 1997 April 21 17:30.
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 speed vlen, and constant acceleration a = (ax, ay, az),
c and for the point b = (bx, by, bz), to find the number of
c intersections nints (0, 1 or 2) of the particle path with point
c "b" that are possible by finding an initial velocity vector
c v(n) with magnitude vlen(n), and for each such intersection,
c to find the velocity vt(n) = (vtx(n), vty(n), vtz(n)), its
c magnitude vtlen(n), the path length to the intersection dint(n),
c and the time to the intersection tint(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 Substituting point "b" for point "pt", solving for v, squaring,
c and rearranging:
c d**2 - ((d dot a) + vlen**2) * t**2 + 0.25 * a**2 * t**4 = 0
c
c This equation may have 0, 1 or 2 real roots tint(n).
c If tint(n) = 0, points "p" and "b" are the same, any vector "v"
c will do, and a vector v(n) = (vlen, 0, 0) will be returned.
c If tint(n) is not zero, vector v(n) is
c v(n) = (d - 0.5 * a * tint(n)**2) / tint(n)
c and the particle velocity at the intersection, vector vt(n) is:
c vt(n) = v(n) + a * tint(n)
c If no real roots exist, point "b" can not be reached with the
c specified value of vlen.
c
c See subroutine aptparb for finding dint.
c
c History:
c
c Input: px, py, pz, vlen, ax, ay, az, bx, by, bz, tol.
c
c Output: nints, vx, vy, vz, vmlen, vtx, vty, vtz, vtlen, dint, tint.
c
c Calls: aptparb, aptqrts, aptvdis, aptvsum, 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 bx,by,bz Input The x, y and z coordinates of point "b".
c
c dint Output The distance along the particle path from point "p" to
c point "b". Size = 2.
c
c nints Output The number of intersections (0, 1 or 2) of the particle
c path and point "b".
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
c intersections with point "b". Size = 2.
c
c tol Input Numerical tolerance limit.
c On Cray computers, recommend 1.e-5 to 1.e-11.
c
c vlen Input The magnitude of the initial particle velocity.
c
c vtlen Output The magnitudes of the particle velocities vt at the
c intersection times. Size = 2.
c
c vtx,y,z Output The x, y and z components of the particle velocity
c at point "b". Size = 2.
c
c vx,vy,vz Output The x, y and z components of the initial particle
c velocity required to intersect point "b". Size = 2.
c
ccend.
c.... Dimensioned arguments.
dimension dint(2) ! Path distance to intersection.
dimension tint(2) ! Time to intersection.
dimension vtlen(2) ! Magnitude of final particle velocity.
dimension vtx(2) ! Component x of final particle velocity.
dimension vty(2) ! Component y of final particle velocity.
dimension vtz(2) ! Component z of final particle velocity.
dimension vx(2) ! Component x of init particle velocity.
dimension vy(2) ! Component y of init particle velocity.
dimension vz(2) ! Component z of init particle velocity.
c.... Local variables.
dimension tsq(2) ! Square of intersection time.
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptparh finding intersections of particle path',
cbug & ' and a point. tol= ',1pe22.14)
cbug 9902 format (' Particle at p with velocity v, acceleration a:' /
cbug & ' px,py,pz= ',1p3e22.14 /
cbug & ' speed = ',1pe22.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, vlen, ax, ay, az
cbug write ( 3, 9903) bx, by, bz
cbugc***DEBUG ends.
c.... Initialize.
nints = 0
do 100 n = 1, 2
dint(n) = -1.e99
tint(n) = -1.e99
vtlen(n) = -1.e99
vtx(n) = -1.e99
vty(n) = -1.e99
vtz(n) = -1.e99
vx(n) = -1.e99
vy(n) = -1.e99
vz(n) = -1.e99
100 continue
c.... Find the vector pb from point "p" to point "b".
call aptvdis (px, py, pz, bx, by, bz, 1, tol,
& pbx, pby, pbz, pblen, 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 dot product of vector "pb" and the acceleration vector "a".
call aptvdot (pbx, pby, pbz, ax, ay, az, 1, tol, pbdota, nerr)
c.... Find the times of intersection.
c0 = pblen**2
c1 = -(pbdota + vlen**2)
c2 = 0.25 * alen**2
cbugcbugc***DEBUG begins.
cbugcbug 9801 format ( 'c0, c1, c2= ',1p3e22.14 )
cbugcbug write ( 3, 9801) c0, c1, c2
cbugcbugc***DEBUG ends.
call aptqrts (0, c2, c1, c0, qq, tol,
& nroots, tsq(1), tsq(2), itrun)
cbugcbugc***DEBUG begins.
cbugcbug 9802 format ( 'nroots=',i2,' tsq1,2=',1p2e22.14 )
cbugcbug write ( 3, 9802) nroots, tsq(1), tsq(2)
cbugcbugc***DEBUG ends.
if (nroots .eq. 0) then ! Test number of roots of t**2.
nints = 0
elseif (nroots .eq. 1) then
if (tsq(1) .lt. 0.0) then
nints = 0
endif
else
if (tsq(1) .lt. 0.0) then
if (tsq(2) .lt. 0.0) then
nints = 0
else
nints = 1
tint(1) = sqrt (tsq(2))
endif
elseif (tsq(2) .lt. 0.0) then
nints = 0
else
nints = 2
tint(1) = sqrt (tsq(1))
tint(2) = sqrt (tsq(2))
endif
endif ! Tested number of roots of t**2.
cbugcbugc***DEBUG begins.
cbugcbug 9803 format ( 'nints=',i2,' tint1,2=',1p2e22.14 )
cbugcbug write ( 3, 9803) nroots, tint(1), tint(2)
cbugcbugc***DEBUG ends.
if (nints .eq. 0) go to 999
c.... Find the intersection velocities and path lengths.
do 140 n = 1, nints
c.... Find the initial velocities v = (pb - 0.5 * a * t**2) / t.
fpb = 1.0 / tint(n)
fa = -0.5 * tint(n)
call aptvsum (0, fpb, pbx, pby, pbz, fa, ax, ay, az, 1, tol,
& vx(n), vy(n), vz(n), vmlenx, nerr)
call aptparb (px, py, pz, vx(n), vy(n), vz(n), ax, ay, az,
& tint(n), tol,
& ptx, pty, ptz,
& vtx(n), vty(n), vtz(n), vtlen(n),
& dint(n), ntype)
140 continue
999 continue
cbugc***DEBUG begins.
cbug 9904 format (/ 'aptparh results: nints=',i2 )
cbug 9905 format (/ 'Init vel(xyz) ',1p3e22.14 /
cbug & ' speed ',1pe22.14 )
cbug 9906 format ( 'Time, path len',1p2e22.14)
cbug 9908 format ( 'Int vel(xyz) ',1p3e22.14 /
cbug & ' speed ',1pe22.14 )
cbug write ( 3, 9904) nints
cbug do 710 n = 1, nints ! Loop over intersection points.
cbug write ( 3, 9905) vx(n), vy(n), vz(n), vlen
cbug write ( 3, 9906) tint(n), dint(n)
cbug write ( 3, 9908) vtx(n), vty(n), vtz(n), vtlen(n)
cbug 710 continue ! End of loop over intersection points.
cbugc***DEBUG ends.
return
c.... End of subroutine aptparh. (+1 line.)
end
UCRL-WEB-209832