subroutine aptparp (px, py, pz, vx, vy, vz, ax, ay, az, & bx, by, bz, cx, cy, cz, tol, & nints, ptx, pty, ptz, vtx, vty, vtz, vtlen, & dint, tint, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTPARP c c call aptparp (px, py, pz, vx, vy, vz, ax, ay, az, c & bx, by, bz, cx, cy, cz, tol, c & nints, ptx, pty, ptz, vtx, vty, vtz, vtlen, c & dint, tint, nerr) c c Version: aptparp Updated 1997 April 16 16:20. c aptparp 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) and constant acceleration c a = (ax, ay, az), and for the plane through point c b = (bx, by, bz) with normal vector c = (cx, cy, cz), to find c the number nints (0, 1 or 2) of intersections of the particle c path with the plane, and for each such intersection point, c the coordinates pt(n) = (ptx(n), pty(n), ptz(n)), n = 1, nints, c the velocities vt(n) = (vtx(n), vty(n), vtz(n)), n = 1, nints, c the velocity magnitudes vtlen(n), n = 1, nints, c the path lengths to intersection dint(n), n = 1, nints, and c the times to intersection tint(n), n = 1, nints. c If no such intersections occur (nints = 0), to find the proximal c points pt(1) = (ptx(1), pty(1), ptz(1)) on the particle path and c pt(2) = (ptx(2), pty(2), ptz(2)) on the plane, the velocity c vt(1) = (vtx(1), vty(1), vtz(1)) and its magnitude vtlen(1) at c the proximal point, the path length to the proximal point c dint(1), the distance between the proximal points dint(2), c and the time tint(1) to reach the proximal point. c Flag nerr indicates any input error. 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 The distance from the particle to the plane is given by: c dp = (pt - b) dot c / |c|, or c dp = [(p - b) dot c + (v dot c)*t + 0.5*(a dot c)*t**2] / |c|, c where "dot" indicates the vector dot product, and |c| is the c absolute magnitude of vector c. c c If any intersections between the particle path and the plane c occur, they may be found by setting dp = 0 and solving for c t = tint in the equation above, then substituting t = tint c in the equation for pt to find d. c There may be zero, one or two such solutions. c If no solutions for t exist, then either the path is parallel c to the plane, with (v dot c) = 0 and (a dot c) = 0, or the c path is at its nearest point to the plane when is it parallel c to the plane: c vt dot c = 0, or c v dot c + (a dot c) tint = 0, or c tint = -(v dot c) / (a dot c), c which may be used to find d = pt at the proximal point. c c See subroutine aptparb for finding dint. c c History: c c Input: px, py, pz, vx, vy, vz, ax, ay, az, bx, by, bz, cx, cy, cz, tol. c c Output: nints, ptx, pty, ptz, vtx, vty, vtz, vtlen, dint, tint, nerr. c c Calls: aptparb, aptptpl, aptqrts, aptvdis, 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 bx,by,bz Input The x, y and z coordinates of point "b" in the plane. c c cx,cy,cz Input The x, y and z components of the vector "c" normal to c the plane. c c dint Output The distance along the particle path to the proximal c point to the plane, and from the plane to the c particle path (nints = 0), or to the intersection(s) c with the plane (nints = 1 or 2). Size = 2. c c nerr Output Indicates an input error, if not zero. c 1 if normal vector c has zero magnitude. c c nints Output The number of intersections of the particle path and c the plane. 0 if no intersections occur. c -1 if the intersection is a tangent point. c c ptx,y,z Output The x, y and z coordinates of the proximal points on c the particle path and in the plane (nints = 0), c or the point(s) of intersection of the particle path c and the plane (nints = 1 or 2). Size 2. c c pz,py,pz Input The initial x, y and z coordinates of the particle. c c tint Output The time along the particle path to the proximal c point to the plane (nints = 0), or to the c intersection(s) with the plane (nints = 1 or 2). c Size = 2. c c tol Input Numerical tolerance limit. c On Cray computers, recommend 1.e-5 to 1.e-11. c c vtlen Output The magnitudes of velocities vt. Size 2. c c vtx,y,z Output The x, y and z components of the particle velocity c at points pt. Size = 2. c c vx,vy,vz Input The initial x, y and z components of the particle c velocity. ccend. c.... Dimensioned arguments. dimension ptx(2) ! Coordinate x of prox or int point. dimension pty(2) ! Coordinate y of prox or int point. dimension ptz(2) ! Coordinate z of prox or int point. dimension vtlen(2) ! Magnitude of particle velocity. dimension vtx(2) ! Component x of particle velocity. dimension vty(2) ! Component y of particle velocity. dimension vtz(2) ! Component z of particle velocity. dimension dint(2) ! Path distance to prox or intersection. dimension tint(2) ! Time to proximity or intersection. cbugc***DEBUG begins. cbug 9901 format (/ 'aptparp finding intersection or proximal point' / cbug & ' of particle path and plane. 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 9903 format (' Plane through point b with normal vector c:' / cbug & ' bx,by,bz= ',1p3e22.14 / cbug & ' cx,cy,cz= ',1p3e22.14 ) cbug write ( 3, 9901) tol cbug write ( 3, 9902) px, py, pz, vx, vy, vz, ax, ay, az cbug write ( 3, 9903) bx, by, bz, cx, cy, cz cbugc***DEBUG ends. c.... Initialize. nints = 0 do 100 n = 1, 2 dint(n) = -1.e99 ptx(n) = -1.e99 pty(n) = -1.e99 ptz(n) = -1.e99 tint(n) = -1.e99 vtlen(n) = -1.e99 vtx(n) = -1.e99 vty(n) = -1.e99 vtz(n) = -1.e99 100 continue c.... Find the magnitude of the normal vector "c". Test for error. call aptvunb (cx, cy, cz, 1, tol, unx, uny, unz, clen, nerr) if (clen .eq. 0.0) then nerr = 1 go to 210 endif c.... Find the magnitude of the initial velocity vector "v". call aptvunb (vx, vy, vz, 1, tol, uvx, uvy, uvz, vlen, 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 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 dot product of vector "pb" and the unit normal vector "un". call aptvdot (pbx, pby, pbz, unx, uny, unz, 1, tol, pbdotu, nerr) c.... Find the dot product of vector "v" and the unit normal vector "un". call aptvdot (vx, vy, vz, unx, uny, unz, 1, tol, vdotu, nerr) c.... Find the dot product of vector "a" and the unit normal vector "un". call aptvdot (ax, ay, az, unx, uny, unz, 1, tol, adotu, nerr) c.... Find any intersections, or if none, find the proximal points. if (adotu .eq. 0.0) then ! Acceleration is parallel to plane. if (vdotu .eq. 0.0) then ! Path always parallel to plane. nints = 0 tint(1) = 0.0 vtx(1) = vx vty(1) = vy vtz(1) = vz vtlen(1) = vlen dint(1) = 0.0 ptx(1) = px pty(1) = py ptz(1) = pz call aptptpl (px, py, pz, bx, by, bz, cx, cy, cz, 1, tol, & dint(2), ptx(2), pty(2), ptz(2), itrun, nerr) else ! Path intersects plane once. nints = 1 tint(1) = pbdotu / vdotu call aptparb (px, py, pz, vx, vy, vz, ax, ay, az, tint(1), & tol, ptx(1), pty(1), ptz(1), & vtx(1), vty(1), vtz(1), vtlen(1), & dint(1), ntype) endif ! Tested direction of initial velocity. else ! Acceleration is not parallel to plane. if (vdotu .eq. 0.0) then ! Path initially parallel to plane. arg = 2.0 * pbdotu / adotu if (arg .lt. 0.0) then ! No intersection. nints = 0 tint(1) = 0.0 vtx(1) = vx vty(1) = vy vtz(1) = vz vtlen(1) = vlen dint(1) = 0.0 dint(2) = -pbdotu ptx(1) = px pty(1) = py ptz(1) = pz call aptptpl (px, py, pz, bx, by, bz, cx, cy, cz, 1, tol, & dint(2), ptx(2), pty(2), ptz(2), itrun, nerr) else ! Two intersections. nints = 2 tint(1) = -sqrt (arg) tint(2) = sqrt (arg) call aptparb (px, py, pz, vx, vy, vz, ax, ay, az, tint(1), & tol, ptx(1), pty(1), ptz(1), & vtx(1), vty(1), vtz(1), vtlen(1), & dint(1), ntype) call aptparb (px, py, pz, vx, vy, vz, ax, ay, az, tint(2), & tol, ptx(2), pty(2), ptz(2), & vtx(2), vty(2), vtz(2), vtlen(2), & dint(2), ntype) endif ! Tested number of intersections. else ! Path is parabolic, not at vertex. aa = 0.5 * adotu bb = vdotu cc = -pbdotu call aptqrts (0, aa, bb, cc, qq, tol, & nroots, troot1, troot2, itrun) if (nroots .eq. 0) then ! No intersections. nints = 0 tint(1) = -vdotu / adotu call aptparb (px, py, pz, vx, vy, vz, ax, ay, az, tint(1), & tol, ptx(1), pty(1), ptz(1), & vtx(1), vty(1), vtz(1), vtlen(1), & dint(1), ntype) call aptptpl (ptx(1), pty(1), ptz(1), & bx, by, bz, cx, cy, cz, 1, tol, & dint(2), ptx(2), pty(2), ptz(2), itrun, nerr) elseif (nroots .eq. 1) then nints = 1 tint(1) = troot1 call aptparb (px, py, pz, vx, vy, vz, ax, ay, az, tint(1), & tol, ptx(1), pty(1), ptz(1), & vtx(1), vty(1), vtz(1), vtlen(1), & dint(1), ntype) c.... See if the intersection is a tangent point. call aptvdot (vtx(1), vty(1), vtz(1), unx, uny, unz, & 1, tol, vtdotu, nerr) if (vtdotu .eq. 0.0) nints = -1 else nints = 2 tint(1) = troot1 tint(2) = troot2 call aptparb (px, py, pz, vx, vy, vz, ax, ay, az, tint(1), & tol, ptx(1), pty(1), ptz(1), & vtx(1), vty(1), vtz(1), vtlen(1), & dint(1), ntype) call aptparb (px, py, pz, vx, vy, vz, ax, ay, az, tint(2), & tol, ptx(2), pty(2), ptz(2), & vtx(2), vty(2), vtz(2), vtlen(2), & dint(2), ntype) endif ! Tested number of roots. endif ! Tested direction of initial velocity. endif ! Tested direction of acceleration. 210 continue cbugc***DEBUG begins. cbug 9904 format (/ 'aptparp results: nerr=',i2,' nints=',i2 / cbug & ' path1,tint1=',1p2e22.14 / cbug & ' d1(x,y,z)= ',1p3e22.14 / cbug & ' path2,tint2=',1p2e22.14 / cbug & ' d2(x,y,z)= ',1p3e22.14 / cbug & ' v1(x,y,z)= ',1p3e22.14 / cbug & ' vtlen1= ',1p1e22.14 / cbug & ' v2(x,y,z)= ',1p3e22.14 / cbug & ' vtlen2= ',1p1e22.14 ) cbug write ( 3, 9904) nerr, nints, cbug & dint(1), tint(1), ptx(1), pty(1), ptz(1), cbug & dint(2), tint(2), ptx(2), pty(2), ptz(2), cbug & vtx(1), vty(1), vtz(1), vtlen(1), cbug & vtx(2), vty(2), vtz(2), vtlen(2) cbugc***DEBUG ends. return c.... End of subroutine aptparp. (+1 line.) end UCRL-WEB-209832