subroutine aptacct (px, py, pz, vx, vy, vz, cx, cy, cz, np, t,
& tol, ptx, pty, ptz, vtx, vty, vtz, vtlen,
& atx, aty, atz, atlen, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTACCT
c
c call aptacct (px, py, pz, vx, vy, vz, cx, cy, cz, np, t,
c & tol, ptx, pty, ptz, vtx, vty, vtz, vtlen,
c & atx, aty, atz, atlen, nerr)
c
c Version: aptacct Updated 1999 February 16 14:30.
c aptacct Originated 1999 February 16 14: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) and
c initial velocity vector v = (vx, vy, vz) at time zero, and with
c time-dependent acceleration vector a = (ax, ay, az), where
c ax = cx(0) + cx(1)*t + ... + cx(n)*t**n + ... + cx(np)*t**np
c ay = cy(0) + cy(1)*t + ... + cy(n)*t**n + ... + cy(np)*t**np
c az = cz(0) + cz(1)*t + ... + cz(n)*t**n + ... + cz(np)*t**np,
c to find at time t the coordinates pt = (ptx, pty, ptz), the
c velocity vector vt = (vtx, vty, vtz), its magnitude vtlen, and
c the acceleration vector at = (atx, aty, atz) and its magnitude
c atlen. The expression for ax may be the initial terms of a
c Taylor series expansion of ax at time zero, with cx(n) being the
c n'th derivative of ax with respect to time, divided by n
c factorial, at time zero. Similarly for ay and az.
c
c The particle velocity and coordinates are given by:
c vt = v + I {a, 0, t},
c pt = p + I {vt, 0, t},
c where I {f, 0, t} indicates the integral of the time-dependent
c function f over time, from 0 to t.
c
c A non-zero value of nerr indicates an input error.
c
c History:
c
c Input: px, py, pz, vx, vy, vz, cx, cy, cz, nt, t, tol.
c
c Output: ptx, pty, ptz, vtx, vty, vtz, vtlen, atx, aty, atz, atlen, nerr.
c
c Calls: aptvunb
c
c Glossary:
c
c atx,... Output The x, y and z components of the acceleration at
c time t.
c
c cx Input The coefficients of the time-dependent x component of
c the acceleration vector a, cx(n), n = 0, np.
c ax = cx(0)+cx(1)*t+...+cx(n)*t**n+...+cx(np)*t**np.
c
c cy Input The coefficients of the time-dependent y component of
c the acceleration vector a, cy(n), n = 1, np.
c ay = cy(0)+cy(1)*t+...+cy(n)*t**n+...+cy(np)*t**np.
c
c cz Input The coefficients of the time-dependent z component of
c the acceleration vector a, cz(n), n = 1, np.
c az = cz(0)+cz(1)*t+...+cz(n)*t**n+...+cz(np)*t**np.
c
c nerr Output Indicates an input error, if not zero:
c 1 if np is negative.
c
c np Input The highest power of t in the expressions for ax,
c ay and az.
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 computers with 64-bit floating point numbers,
c 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.
c.... Dimensioned arguments.
dimension cx(0:1)
dimension cy(0:1)
dimension cz(0:1)
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptacct 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 9903 format (' cx= ',1p3e22.14 /
cbug & (' ',1p3e22.14 ))
cbug 9904 format (' cy= ',1p3e22.14 /
cbug & (' ',1p3e22.14 ))
cbug 9905 format (' cz= ',1p3e22.14 /
cbug & (' ',1p3e22.14 ))
cbug write ( 3, 9901) t, tol
cbug write ( 3, 9902) px, py, pz, vx, vy, vz
cbug write ( 3, 9903) (cx(n), n = 0, np)
cbug write ( 3, 9904) (cy(n), n = 0, np)
cbug write ( 3, 9905) (cz(n), n = 0, np)
cbugc***DEBUG ends.
c.... Initialize.
nerr = 0
c.... Test for input errors.
if (np .lt. 0) then
nerr = 1
go to 210
endif
c.... Find the magnitude of the initial acceleration vector.
call aptvunb (cx(0), cy(0), cz(0), 1, tol,
& ux, uy, uz, alen, nerr)
c.... Find the magnitude of the initial velocity vector.
call aptvunb (vx, vy, vz, 1, tol, ux, uy, uz, vlen, nerr)
c.... Test for zero time interval.
if (t .eq. 0.0) then
atx = cx(0)
aty = cy(0)
atz = cz(0)
atlen = alen
vtx = vx
vty = vy
vtz = vz
vtlen = vlen
ptx = px
pty = py
ptz = pz
go to 210
endif
c.... Find the velocity and position at time t.
if (np .ge. 1) then
atx = cx(0)
aty = cy(0)
atz = cz(0)
vtx = vx + cx(0) * t
vty = vy + cy(0) * t
vtz = vz + cz(0) * t
ptx = px + vx * t + 0.5 * cx(0) * t**2
pty = py + vy * t + 0.5 * cy(0) * t**2
ptz = pz + vz * t + 0.5 * cz(0) * t**2
ta = 1.0
do 120 n = 1, np ! Loop over coefficients.
ta = ta * t
tv = ta * t
tp = tv * t
fv = n + 1
fp = (n + 1) * (n + 2)
atx = atx + cx(n) * ta
aty = aty + cy(n) * ta
atz = atz + cz(n) * ta
vtx = vtx + cx(n) * tv / fv
vty = vty + cy(n) * tv / fv
vtz = vtz + cz(n) * tv / fv
ptx = ptx + cx(n) * tp / fp
pty = pty + cy(n) * tp / fp
ptz = ptz + cz(n) * tp / fp
120 continue
endif ! End of loop over coefficients.
c.... Truncate results near zero.
if (tol .gt. 0.0) then ! Truncate results.
atxerr = abs (cx(0))
atyerr = abs (cy(0))
atzerr = abs (cz(0))
vtxerr = abs (vx) + 2.0 * abs (cx(0) * t)
vtyerr = abs (vy) + 2.0 * abs (cy(0) * t)
vtzerr = abs (vz) + 2.0 * abs (cz(0) * t)
ptxerr = abs (px) + 1.5 * abs (vx * t) * abs (cx(0)) * t**2
ptyerr = abs (py) + 1.5 * abs (vy * t) * abs (cy(0)) * t**2
ptzerr = abs (pz) + 1.5 * abs (vz * t) * abs (cz(0)) * t**2
if (np .ge. 1) then
ta = 1.0
tv = t
do 130 n = 1, np ! Loop over coefficients.
ta = ta * t
tv = ta * t
tp = tv * t
fv = n + 1
fp = (n + 1) * (n + 2)
atxerr = atxerr + fv * abs (cx(n) * ta)
atyerr = atyerr + fv * abs (cy(n) * ta)
atzerr = atzerr + fv * abs (cz(n) * ta)
vtxerr = vtxerr + (fv + 1.0) * abs (cx(n) * tv) / fv
vtyerr = vtyerr + (fv + 1.0) * abs (cy(n) * tv) / fv
vtzerr = vtzerr + (fv + 1.0) * abs (cz(n) * tv) / fv
ptxerr = ptxerr + (fv + 2.0) * abs (cx(n) * tp) / fp
ptyerr = ptyerr + (fv + 2.0) * abs (cy(n) * tp) / fp
ptzerr = ptzerr + (fv + 2.0) * abs (cz(n) * tp) / fp
130 continue
endif ! End of loop over coefficients.
if (abs (atx) .le. tol * atxerr) atx = 0.0
if (abs (aty) .le. tol * atyerr) aty = 0.0
if (abs (atz) .le. tol * atzerr) atz = 0.0
if (abs (vtx) .le. tol * vtxerr) vtx = 0.0
if (abs (vty) .le. tol * vtyerr) vty = 0.0
if (abs (vtz) .le. tol * vtzerr) vtz = 0.0
if (abs (ptx) .le. tol * ptxerr) ptx = 0.0
if (abs (pty) .le. tol * ptyerr) pty = 0.0
if (abs (ptz) .le. tol * ptzerr) ptz = 0.0
endif ! Truncate results.
c.... Find the magnitude of the final acceleration vector.
call aptvunb (atx, aty, atz, 1, tol, ux, uy, uz, atlen, nerr)
c.... Find the magnitude of the final velocity vector.
call aptvunb (vtx, vty, vtz, 1, tol, ux, uy, uz, vtlen, nerr)
210 continue
cbugc***DEBUG begins.
cbug 9906 format (/ 'aptacct results: nerr=',i2 /
cbug & ' pt(x,y,z)= ',1p3e22.14 /
cbug & ' vt(x,y,z)= ',1p3e22.14 /
cbug & ' vtlen= ',1pe22.14 /
cbug & ' at(x,y,z)= ',1p3e22.14 /
cbug & ' atlen= ',1pe22.14 )
cbug write ( 3, 9906) nerr, ptx, pty, ptz, vtx, vty, vtz, vtlen,
cbug & atx, aty, atz, atlen
cbug 140 continue
cbugc***DEBUG ends.
return
c.... End of subroutine aptacct. (+1 line.)
end
UCRL-WEB-209832