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