subroutine aptparb (px, py, pz, vx, vy, vz, ax, ay, az, t, tol,
& ptx, pty, ptz, vtx, vty, vtz, vtlen,
& dpath, ntype)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTPARB
c
c call aptparb (px, py, pz, vx, vy, vz, ax, ay, az, t, tol,
c & ptx, pty, ptz, vtx, vty, vtz, vtlen, dpath, ntype)
c
c Version: aptparb Updated 1997 April 11 15:00.
c aptparb Originated 1996 December 19 18:00.
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) at time zero, and constant
c acceleration a = (ax, ay, az), to find at time t the coordinates
c pt = (ptx, pty, ptz), the velocity vt = (vtx, vty, vtz) and
c its magnitude vtlen, and the path length dpath.
c Flag ntype indicates the path is linear (1) or parabolic (2).
c
c The particle velocity and coordinates at time t are given by:
c vt = v + a * t,
c pt = p + v * t + 0.5 * a * t**2.
c
c The path length is given by the integral from time 0 to t of
c the absolute value of vt:
c dpath = integral {|v + a * t| dt}, from time = 0 to t,
c where
c |v + a * t|**2 = v dot v + 2 * v dot a * t + a dot a * t**2,
c where "dot" indicates the vector dot product.
c
c History:
c
c Input: px, py, pz, vx, vy, vz, ax, ay, az, t, tol.
c
c Output: ptx, pty, ptz, vtx, vty, vtz, vtlen, dpath.
c
c Calls: aptvdot, 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 dpath Output The distance along the particle path from time 0 to t.
c
c ntype Output If 1, the path is linear, but may have a vertex.
c If 2, the path is parabolic, in the plane containing
c point p and vectors v and a.
c
c ptx,... Output The x, y and z coordinates of the particle at time t.
c
c pz,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 vtlen Output The magnitude of the velocity at time t.
c
c vtx,... Output The x, y and z components of the particle velocity at
c time t.
c
c vx,vy,vz Input The initial x, y and z components of the particle
c velocity.
ccend.
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptparb finding particle path at time',1pe22.14 /
cbug & ' 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 write ( 3, 9901) t, tol
cbug write ( 3, 9902) px, py, pz, vx, vy, vz, ax, ay, az
cbugc***DEBUG ends.
c.... Initialize.
ntype = 1
c.... Find the magnitude of the initial velocity vector.
call aptvunb (vx, vy, vz, 1, tol, ux, uy, uz, vlen, nerr)
c.... Find the magnitude of the acceleration vector.
call aptvunb (ax, ay, az, 1, tol, ux, uy, uz, alen, nerr)
c.... Find the velocity at time t.
call aptvsum (0, 1.0, vx, vy, vz, t, ax, ay, az, 1, tol,
& vtx, vty, vtz, vtlen, nerr)
c.... Find the average velocity.
tmid = 0.5 * t
call aptvsum (0, 1.0, vx, vy, vz, tmid, ax, ay, az, 1, tol,
& vmx, vmy, vmz, vmlen, nerr)
c.... Find the new position.
call aptvsum (0, 1.0, px, py, pz, t, vmx, vmy, vmz, 1, tol,
& ptx, pty, ptz, ptlen, nerra)
c.... Find the arc length from time 0 to t.
if (alen .eq. 0.0) then
dpath = vmlen * t
elseif (vlen .eq. 0.0) then
dpath = vmlen * t
else
call aptvdot (vx, vy, vz, ax, ay, az, 1, tol, vdota, nerr)
tver = -vdota / alen**2 ! Time to reach vertex of curve.
qq = vlen**2 + vdota * tver
qqerr = tol * (vlen**2 - vdota * tver)
if (abs (qq) .le. qqerr) qq = 0.0
if (qq .eq. 0.0) then
dpath = 0.5 * (vlen * tver + vtlen * (t - tver))
epath = 0.5 * tol * (vlen * abs (tver) +
& vtlen * (abs (t) + abs (tver)))
if (abs (dpath) .le. epath) dpath = 0.0
else
ntype = 2
fl1 = log (vlen - alen * tver)
fl2 = log (vtlen + alen * (t - tver))
dpath = 0.5 * (vlen * tver + vtlen * (t - tver)) +
& 0.5 * qq * (fl2 - fl1) / alen
epath = 0.5 * tol * (vlen * abs (tver) +
& vtlen * (abs (t) + abs (tver)) +
& 0.5 * abs (qq) * (abs (fl2) + abs (fl1))) / alen
if (abs (dpath) .le. epath) dpath = 0.0
endif
endif
210 continue
cbugc***DEBUG begins.
cbug 9903 format (/ 'aptparb results: ntype=',i2 /
cbug & ' pt(x,y,z)= ',1p3e22.14 /
cbug & ' vt(x,y,z)= ',1p3e22.14 /
cbug & ' vtlen= ',1pe22.14 /
cbug & ' dpath= ',1pe22.14 )
cbug write ( 3, 9903) ntype, ptx, pty, ptz, vtx, vty, vtz, vtlen,
cbug & dpath
cbug 140 continue
cbugc***DEBUG ends.
return
c.... End of subroutine aptparb. (+1 line.)
end
UCRL-WEB-209832