subroutine aptpart (ptx, pty, ptz, time, np, tol,
& px, py, pz, vx, vy, vz, vlenv,
& ax, ay, az, vlena, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTPART
c
c call aptpart (ptx, pty, ptz, time, np, tol,
c px, py, pz, vx, vy, vz, vlenv,
c ax, ay, az, vlena, nerr)
c
c Version: aptpart Updated 1997 April 15 15:30.
c aptpart Originated 1997 April 14 17:30.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: For the particle with the specified coordinates
c pt(n) = (ptx(n), pty(n), ptz(n)) at times time(n), n = 1, np
c (np = 1, 2 or 3), on a uniformly accelerated parabolic
c trajectory, find the initial coordinates p = (px, py, pz)
c and velocity v = (vx, vy, vz) (zero if np = 1) at time zero,
c and the constant acceleration a = (ax, ay, az) (zero if
c np = 1 or 2), and the magnitudes of the velocity and
c acceleration vectors.
c Flag nerr indicates an input error, if nonzero.
c
c The particle coordinates are given by:
c pt(n) = p + v * t(n) + 0.5 * a * t(n)**2.
c These equations may be solved for the 1, 2 or 3 input
c points and times, to find p, v and a.
c
c
c History:
c
c Input: ptx, pty, ptz, time, np, tol.
c
c Output: px, py, pz, vx, vy, vz, vlen, ax, ay, az, vlena, nerr
c
c Calls: aptvunb
c
c
c Glossary:
c
c ax,ay,az Output The x, y and z components of the constant acceleration
c vector a (zero if np = 1 or 2).
c
c nerr Output Indicates an input error, if positive.
c 1 if np < 1 or np > 3.
c 2 if any two times are the same.
c
c np Input The number of points and times specified on the
c parabolic trajectory of the uniformly accelerated
c particle. Must be 2 or 3. If 2, the uniform
c acceleration is zero, and the velocity is constant.
c
c ptx,... Output The x, y and z coordinates of the particle at times
c time(n), n = 1, np. Size 2 or 3.
c
c pz,py,pz Output The x, y and z coordinates of the particle at time
c zero.
c
c time Input The times corresponding to points pt(n), n = 1, np.
c Size 2 or 3.
c
c tol Input Numerical tolerance limit.
c On Cray computers, recommend 1.e-5 to 1.e-11.
c
c vlena Output The magnitude of the constant acceleration vector.
c
c vlenv Output The magnitude of the initial velocity vector.
c
c vx,vy,vz Output The x, y and z components of the initial particle
c velocity (zero if np = 1).
ccend.
c.... Dimensioned arguments.
dimension ptx(3) ! Coordinate x of a particle position.
dimension pty(3) ! Coordinate y of a particle position.
dimension ptz(3) ! Coordinate z of a particle position.
dimension time(3) ! Time of a particle position.
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptpart finding initial particle position',
cbug & ' and velocity' /
cbug & ' and constant acceleration. tol=',1pe22.14 /
cbug & ' Particle positions and times:' )
cbug 9902 format (' ptx,pty,ptz=',1p3e22.14 /
cbug & ' t= ',1pe22.14 )
cbug write ( 3, 9901) tol
cbug write ( 3, 9902) (ptx(n), pty(n), ptz(n), time(n), n = 1, np)
cbugc***DEBUG ends.
c.... Initialize.
nerr = 0
px = -1.e99
py = -1.e99
pz = -1.e99
vx = -1.e99
vy = -1.e99
vz = -1.e99
vlenv = -1.e99
ax = -1.e99
ay = -1.e99
az = -1.e99
vlenz = -1.e99
if ((np .lt. 1) .or. (np .gt. 3)) then
nerr = 1
go to 999
endif
if (np .eq. 1) then ! One point. No velocity or acceleration.
px = ptx(1)
py = pty(1)
pz = ptz(1)
vx = 0.0
vy = 0.0
vy = 0.0
vlenv = 0.0
ax = 0.0
ay = 0.0
ay = 0.0
vlena = 0.0
elseif (np .eq. 2) then ! Two points. No acceleration.
x1 = ptx(1)
y1 = pty(1)
z1 = ptz(1)
t1 = time(1)
x2 = ptx(2)
y2 = pty(2)
z2 = ptz(2)
t2 = time(2)
dt12 = t2 - t1
e12 = tol * (abs (t1) + abs (t2))
if (abs (dt12) .le. e12) dt12 = 0.0
if (dt12 .eq. 0.0) then
nerr = 2
go to 999
endif
if (nerr .ne. 0) go to 999
c.... Find the position at time zero.
px = (t2 * x1 - t1 * x2) / dt12
ex = tol * (abs (t2 * x1) + abs (t1 * x2)) / abs (dt12)
if (abs (px) .le. ex) px = 0.0
py = (t2 * y1 - t1 * y2) / dt12
ey = tol * (abs (t2 * y1) + abs (t1 * y2)) / abs (dt12)
if (abs (py) .le. ey) py = 0.0
pz = (t2 * z1 - t1 * z2) / dt12
ez = tol * (abs (t2 * z1) + abs (t1 * z2)) / abs (dt12)
if (abs (pz) .le. ez) pz = 0.0
c.... Find the velocity at time zero.
vx = (x2 - x1) / dt12
ex = tol * (abs (x2) + abs (x1)) / dt12
if (abs (vx) .le. ex) vx = 0.0
vy = (y2 - y1) / dt12
ey = tol * (abs (y2) + abs (y1)) / dt12
if (abs (vy) .le. ey) vy = 0.0
vz = (z2 - z1) / dt12
ez = tol * (abs (z2) + abs (z1)) / dt12
if (abs (vz) .le. ez) vz = 0.0
call aptvunb (vx, vy, vz, 1, tol,
& ux, uy, uz, vlenv, nerra)
c.... Find the acceleration.
ax = 0.0
ay = 0.0
az = 0.0
vlena = 0.0
elseif (np .eq. 3) then ! Three points.
x1 = ptx(1)
y1 = pty(1)
z1 = ptz(1)
t1 = time(1)
x2 = ptx(2)
y2 = pty(2)
z2 = ptz(2)
t2 = time(2)
x3 = ptx(3)
y3 = pty(3)
z3 = ptz(3)
t3 = time(3)
dt12 = t2 - t1
e12 = tol * (abs (t1) + abs (t2))
if (abs (dt12) .le. e12) dt12 = 0.0
dt23 = t3 - t2
e23 = tol * (abs (t2) + abs (t3))
if (abs (dt23) .le. e23) dt23 = 0.0
dt31 = t1 - t3
e31 = tol * (abs (t3) + abs (t1))
if (abs (dt31) .le. e31) dt31 = 0.0
if ((dt12 .eq. 0.0) .or.
& (dt23 .eq. 0.0) .or.
& (dt31 .eq. 0.0)) then
nerr = 2
go to 999
endif
c.... Find the particle position at time zero.
c1 = -t2 * t3 / (dt31 * dt12)
c2 = -t3 * t1 / (dt12 * dt23)
c3 = -t1 * t2 / (dt23 * dt31)
px = (c1 * x1 + c2 * x2 + c3 * x3)
ex = tol * (abs (c1 * x1) + abs (c2 * x2) + abs (c3 * x3))
if (abs (px) .le. ex) px = 0.0
py = (c1 * y1 + c2 * y2 + c3 * y3)
ey = tol * (abs (c1 * y1) + abs (c2 * y2) + abs (c3 * y3))
if (abs (py) .le. ey) py = 0.0
pz = (c1 * z1 + c2 * z2 + c3 * z3)
ez = tol * (abs (c1 * z1) + abs (c2 * z2) + abs (c3 * z3))
if (abs (pz) .le. ez) pz = 0.0
c.... Find the velocity at time zero.
c1 = (t2 + t3) / (dt31 * dt12)
e1 = tol * (abs (t2) + abs (t3)) / (dt12 * dt31)
if (abs (c1) .le. e1) c1 = 0.0
c2 = (t3 + t1) / (dt12 * dt23)
e2 = tol * (abs (t3) + abs (t1)) / (dt23 * dt12)
if (abs (c2) .le. e2) c2 = 0.0
c3 = (t1 + t2) / (dt23 * dt31)
e3 = tol * (abs (t1) + abs (t2)) / (dt31 * dt23)
if (abs (c3) .le. e3) c3 = 0.0
vx = c1 * x1 + c2 * x2 + c3 * x3
ex = tol * (abs (c1 * x1) + abs (c2 * x2) + abs (c3 * x3))
if (abs (vx) .le. ex) vx = 0.0
vy = c1 * y1 + c2 * y2 + c3 * y3
ey = tol * (abs (c1 * y1) + abs (c2 * y2) + abs (c3 * y3))
if (abs (vy) .le. ey) vy = 0.0
vz = c1 * z1 + c2 * z2 + c3 * z3
ez = tol * (abs (c1 * z1) + abs (c2 * z2) + abs (c3 * z3))
if (abs (vz) .le. ez) vz = 0.0
call aptvunb (vx, vy, vz, 1, tol,
& ux, uy, uz, vlenv, nerra)
c.... Find the constant acceleration.
c1 = -2.0 / (dt31 * dt12)
c2 = -2.0 / (dt12 * dt23)
c3 = -2.0 / (dt23 * dt31)
ax = c1 * x1 + c2 * x2 + c3 * x3
ex = tol * (abs (c1 * x1) + abs (c2 * x2) + abs (c3 * x3))
if (abs (ax) .le. ex) ax = 0.0
ay = c1 * y1 + c2 * y2 + c3 * y3
ey = tol * (abs (c1 * y1) + abs (c2 * y2) + abs (c3 * y3))
if (abs (ay) .le. ey) ay = 0.0
az = c1 * z1 + c2 * z2 + c3 * z3
ez = tol * (abs (c1 * z1) + abs (c2 * z2) + abs (c3 * z3))
if (abs (az) .le. ez) az = 0.0
call aptvunb (ax, ay, az, 1, tol,
& ux, uy, uz, vlena, nerra)
endif ! Tested number of points.
999 continue
cbugc***DEBUG begins.
cbug 9903 format (/ 'aptpart results: nerr=',i2 /
cbug & ' p(x,y,z)= ',1p3e22.14 /
cbug & ' v(x,y,z)= ',1p3e22.14 /
cbug & ' vlenv= ',1pe22.14 /
cbug & ' a(x,y,z)= ',1p3e22.14 /
cbug & ' vlena= ',1pe22.14 )
cbug write ( 3, 9903) nerr, px, py, pz, vx, vy, vz, vlenv,
cbug & ax, ay, az, vlena
cbug 140 continue
cbugc***DEBUG ends.
return
c.... End of subroutine aptpart. (+1 line.)
end
UCRL-WEB-209832