subroutine aptparx (px, py, pz, vx, vy, vz, ax, ay, az, & qc, qx, qy, qz, qxy, qyz, qzx, & qxx, qyy, qzz, tol, & nprox, atype, ptx, pty, ptz, & vtx, vty, vtz, vtlen, dprox, tprox, & dside, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTPARX c c call aptparx (px, py, pz, vx, vy, vz, ax, ay, az, c & qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz, tol, c & nprox, atype, ptx, pty, ptz, c & vtx, vty, vtz, vtlen, dprox, tprox, c & dside, nerr) c c Version: aptparx Updated 1997 May 6 15:25. c aptparx Originated 1997 May 6 15:25. 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 quadric surface with the implicit c equation F(x,y,z) = qc + qx*x + qy*y + qz*z + qxy*x*y + c qyz*y*z + qzx*z*x + qxx*x**2 + qyy*y**2 + qzz*z**2 = 0, c to find the number nprox (0, 1, 2 or 3) of extrema of F on c the particle path, c the extremum type ("minimum", "min/max" or "maximum"), c the coordinates pt(n) = (ptx(n), pty(n), ptz(n)), c the velocities vt(n) = (vtx(n), vty(n), vtz(n)), c the velocity magnitudes vtlen(n), c the path lengths to extremum dprox(n), c the times to extremum tprox(n), and c the value dside(n) of F(x,y,z), all for n = 1, nprox. c Flag nerr indicates any input error. c c The x, y and z components of the particle velocity and c coordinates are given by: c vtx = vx + ax*t c vty = vy + ay*t c vtz = vz + az*t c ptx = px + vx*t + 0.5*ax*t c pty = py + vy*t + 0.5*ay*t c ptz = pz + vz*t + 0.5*az*t c c The times of any extrema of the particle path with the c quadric surface, tprox, are obtained by substituting ptx, pty c and ptz for x, y and z in the equation of the quadric surface, c differentiating with respect to time, and solving the resulting c cubic, quadratic or linear equation for the times. c There may be 0, 1, 2 or 3 such extrema. c The path distances dprox are found in aptparb. c c Flag nerr indicates any input error. c c History: c c Input: px, py, pz, vx, vy, vz, qc, qx, qy, qz, qxy, qyz, qzx, c qxx, qyy, qzz, tol. c c Output: nprox, atype, ptx, pty, ptz, vtx, vty, vtz, vtlen, c dprox, tprox, nerr. c c Calls: aptcubs, aptparb, aptqrts, aptquar c c c Glossary: c c atype Output The extrema types, 'minimum' or 'maximum'. c Size = 3. c c ax,ay,az Input The x, y and z components of the constant acceleration c vector "a". c c dprox Output The distance along the particle path to the quadric c surface, if nprox is positive. Size = 3. c c nerr Output Indicates an input error, if not zero. c 1 if equation for extrema times is null, indicating c particle path is a point. If so, nprox = 1, and c extremum is the particle point. c 2 if equation for extrema times is impossible. c c nprox Output The number of extrema of the particle path and c the quadric surface. Range 0 to 3. c c ptx,y,z Output The x, y and z coordinates of the extrema on the c particle path, if nprox is positive. Size = 3. c c pz,py,pz Input The initial x, y and z coordinates of the particle. c c qc-qzz Input The coefficients of the implicit equation for the c quartic surface. c c tprox Output The times along the particle path to the extrema, c if nprox is positive. Size = 3. c c tol Input Numerical tolerance limit. With 64-bit floating point c numbers, recommend 1.e-5 to 1.e-11. c c vtx,y,z Output The x, y and z components of the particle velocity c at points pt, if nprox is positive. Size = 3. c c vx,vy,vz Input The initial x, y and z components of the particle c velocity. ccend. c.... Dimensioned arguments. dimension atype(3) ! Extremum type 'minimum' or 'maximum'. character*8 atype ! Extremum type 'minimum' or 'maximum'. dimension dprox(3) ! Path distance to extremum. dimension dside(3) ! Value of quadric equation. dimension ptx(3) ! Coordinate x of int point. dimension pty(3) ! Coordinate y of int point. dimension ptz(3) ! Coordinate z of int point. dimension tprox(3) ! Time to extremum. dimension timag(3) ! Imaginary component of time. dimension vtlen(3) ! Magnitude of particle velocity. dimension vtx(3) ! Component x of particle velocity. dimension vty(3) ! Component y of particle velocity. dimension vtz(3) ! Component z of particle velocity. cbugc***DEBUG begins. cbug 9901 format (/ 'aptparx finding extrema of quadric surface' / cbug & ' equation on particle path.', 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 9903 format (' Quadric surface with implicit equation', cbug & ' coefficients:' / cbug & ' qc ',1pe22.14 / cbug & ' qx,qy,qz= ',1p3e22.14 / cbug & ' qxy,qyz,qzx=',1p3e22.14 / cbug & ' qxx,qyy,qzz=',1p3e22.14 ) cbug write ( 3, 9901) tol cbug write ( 3, 9902) px, py, pz, vx, vy, vz, ax, ay, az cbug write ( 3, 9903) qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz cbugc***DEBUG ends. c.... Initialize. nprox = 0 do 100 n = 1, 3 atype(n) = 'unknown' dprox(n) = -1.e99 ptx(n) = -1.e99 pty(n) = -1.e99 ptz(n) = -1.e99 tprox(n) = -1.e99 timag(n) = -1.e99 vtlen(n) = -1.e99 vtx(n) = -1.e99 vty(n) = -1.e99 vtz(n) = -1.e99 100 continue c.... Find the coefficients of the times in the implicit quadric equation. c0 = qc + & qx *px + qy *py + qz *pz + & qxy*px*py + qyz*py*pz + qzx*pz*px + & qxx*px*px + qyy*py*py + qzz*pz*pz e0 = abs (qc) + & abs (qx *px) + abs (qy *py) + abs (qz *pz) + & abs (qxy*px*py) + abs (qyz*py*pz) + abs (qzx*pz*px) + & abs (qxx*px*px) + abs (qyy*py*py) + abs (qzz*pz*pz) if (abs (c0) .le. tol * e0) c0 = 0.0 c1 = qx*vx + qy*vy + qz*vz + & qxy*(px*vy + py*vx) + & qyz*(py*vz + pz*vy) + & qzx*(pz*vx + px*vz) + & 2.0*(qxx*px*vx + qyy*py*vy + qzz*pz*vz) e1 = abs (qx*vx) + abs (qy*vy) + abs (qz*vz) + & abs (qxy)*(abs (px*vy) + abs (py*vx)) + & abs (qyz)*(abs (py*vz) + abs (pz*vy)) + & abs (qzx)*(abs (pz*vx) + abs (px*vz)) + & 2.0*(abs (qxx*px*vx) + abs (qyy*py*vy) + abs (qzz*pz*vz)) if (abs (c1) .le. tol * e1) c1 = 0.0 c2 = 0.5* (qx*ax + qy*ay + qz*az) + & 0.5*qxy*(px*ay + 2.0*vx*vy + py*ax) + & 0.5*qyz*(py*az + 2.0*vy*vz + pz*ay) + & 0.5*qzx*(pz*ax + 2.0*vz*vx + px*az) + & qxx*(px*ax + vx*vx) + & qyy*(py*ay + vy*vy) + & qzz*(pz*az + vz*vz) e2 = 0.5* ((abs (qx*ax) + abs (qy*ay) + abs (qz*az)) + & abs (qxy)*(abs (px*ay) + 2.0*abs (vx*vy) + abs (py*ax)) + & abs (qyz)*(abs (py*az) + 2.0*abs (vy*vz) + abs (pz*ay)) + & abs (qzx)*(abs (pz*ax) + 2.0*abs (vz*vx) + abs (px*az))) + & abs (qxx)*(abs (px*ax) + abs (vx*vx)) + & abs (qyy)*(abs (py*ay) + abs (vy*vy)) + & abs (qzz)*(abs (pz*az) + abs (vz*vz)) if (abs (c2) .le. tol * e2) c2 = 0.0 c3 = 0.5*qxy*(vx*ay + vy*ax) + & 0.5*qyz*(vy*az + vz*ay) + & 0.5*qzx*(vz*ax + vx*az) + & qxx*vx*ax + qyy*vy*ay + qzz*vz*az e3 = 0.5*abs (qxy)*(abs (vx*ay) + abs (vy*ax)) + & 0.5*abs (qyz)*(abs (vy*az) + abs (vz*ay)) + & 0.5*abs (qzx)*(abs (vz*ax) + abs (vx*az)) + & abs (qxx*vx*ax) + abs (qyy*vy*ay) + abs (qzz*vz*az) if (abs (c3) .le. tol * e3) c3 = 0.0 c4 = 0.25*(qxy*ax*ay + qyz*ay*az + qzx*az*ax + & qxx*ax*ax + qyy*ay*ay + qzz*az*az) e4 = 0.25*(abs (qxy*ax*ay) + abs (qyz*ay*az) + abs (qzx*az*ax) + & abs (qxx*ax*ax) + abs (qyy*ay*ay) + abs (qzz*az*az)) if (abs (c4) .le. tol * e4) c4 = 0.0 c.... Find any extrema of the particle path and the quadric. d0 = c1 d1 = 2.0 * c2 d2 = 3.0 * c3 d3 = 4.0 * c4 cbugcbugc***DEBUG begins. cbugcbug 9904 format ('Coefficients of the equation for the', cbugcbug & ' extremum times:' / cbugcbug & ' set d0 =',1pe20.12 / cbugcbug & ' set d1 =',1pe20.12 / cbugcbug & ' set d2 =',1pe20.12 / cbugcbug & ' set d3 =',1pe20.12 ) cbugcbug cbugcbug write ( 3, 9904) d0, d1, d2, d3 cbugcbugc***DEBUG ends. if (d3 .ne. 0.0) then ! Solve cubic equation. cbugcbugc***DEBUG begins. cbugcbug 9906 format ('aptparx DEBUG. Solving cubic.') cbugcbug write ( 3, 9906) cbugcbugc***DEBUG ends d0 = d0 / d3 d1 = d1 / d3 d2 = d2 / d3 call aptcubs (d0, d1, d2, tol, & nprox, tprox, timag, atype, itrun) elseif (d2 .ne. 0.0) then ! Solve quadratic equation. cbugcbugc***DEBUG begins. cbugcbug 9907 format ('aptparx DEBUG. Solving quadratic.') cbugcbug write ( 3, 9907) cbugcbugc***DEBUG ends call aptqrts (0, d2, d1, d0, qq, tol, & nprox, tprox(1), tprox(2), itrun) elseif (d1 .ne. 0.0) then ! Solve linear equation. cbugcbugc***DEBUG begins. cbugcbug 9908 format ('aptparx DEBUG. Solving linear.') cbugcbug write ( 3, 9908) cbugcbugc***DEBUG ends nprox = 1 tprox(1) = -d0 / d1 elseif (d0 .ne. 0.0) then ! Impossible equation. nerr = 2 cbugcbugc***DEBUG begins. cbugcbug 9909 format ('aptparx DEBUG. Impossible equation.') cbugcbug write ( 3, 9909) cbugcbugc***DEBUG ends else ! Null equation. Path is a point. nerr = 1 nprox = 1 tprox(1) = 0.0 cbugcbugc***DEBUG begins. cbugcbug 9910 format ('aptparx DEBUG. Null equation.') cbugcbug write ( 3, 9910) cbugcbugc***DEBUG ends endif ! Tested order of equation for time. c.... Find path, speed, velocity and position at extrema. if (nprox .gt. 0) then do 850 n = 1, nprox ! Loop over times of extrema. curve = d1 + 2.0 * d2 * tprox(n) + 3.0 * d3 * tprox(n)**2 curverr = tol * (abs (d1) + 2.0 * abs (d2 * tprox(n)) + & 3.0 * abs (d3) * tprox(n)**2) if (abs (curve) .le. curverr) curve = 0.0 if (curve .lt. 0.0) then atype(n) ='maximum' elseif (curve .eq. 0.0) then atype(n) = 'min/max' else atype(n) ='minimum' endif call aptparb (px, py, pz, vx, vy, vz, ax, ay, az, & tprox(n), tol, & ptx(n), pty(n), ptz(n), & vtx(n), vty(n), vtz(n), & vtlen(n), dprox(n), ntype) ta = tprox(n) dside(n) = c0 + ta * (c1 + ta * (c2 + ta * (c3 + ta * c4))) ta = abs (tprox(n)) eside = tol * (abs (c0) + ta * (abs (c1) + ta * & (abs (c2) + ta * (abs (c3) + ta * abs (c4))))) if (abs (dside(n)) .lt. eside) dside(n) = 0.0 850 continue ! End of loop over times of extrema. endif cbugc***DEBUG begins. cbug 9912 format (/ 'aptparx results: nerr=',i2,' nprox=',i2 ) cbug write ( 3, 9912) nerr, nprox cbug if (nprox .gt. 0) then ! One or more extrema. cbug 9913 format (/ 'Extremum type ',a8 / cbug & 'Time, path, speed ',1p3e20.12 / cbug & 'Velocity (xyz) ',1p3e20.12 / cbug & 'Position (xyz) ',1p3e20.12 / cbug & 'dside ',1pe20.12 ) cbug write ( 3, 9913) (atype(n), tprox(n), dprox(n), vtlen(n), cbug & vtx(n), vty(n), vtz(n), cbug & ptx(n), pty(n), ptz(n), dside(n), cbug & n = 1, nprox) cbug endif ! Tested number of extrema. cbugc***DEBUG ends. return c.... End of subroutine aptparx. (+1 line.) end UCRL-WEB-209832