subroutine aptqprr (ax, ay, az, bx, by, bz, qc, qx, qy, qz, & qxy, qyz, qzx, qxx, qyy, qzz, nran, tol, & nprox, cx, cy, cz, dist, dx, dy, dz, & snx, sny, snz, costh, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTQPRR c c call aptqprr (ax, ay, az, bx, by, bz, qc, qx, qy, qz, c qxy, qyz, qzx, qxx, qyy, qzz, nran, tol, c nprox, cx, cy, cz, dist, dx, dy, dz, c snx, sny, snz, costh, nerr) c c Version: aptqprr Updated 1994 December 20 12:50. c aptqprr Originated 1994 October 31 13:40. c c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123. c c c Purpose: To try to find the proximal point "c" = (cx, cy, cz) on a c quadric surface, nearest point "a" = (ax, ay, az), with an c initial trial direction "b" = (bx, by, bz). c c Note: use APT subroutine aptwhis to find the proximal point. c Use this subroutine only if aptwhis fails. c c The quadric surface is represented by: c c f(x,y,z) = qc + qx * x + qy * y + qz * z + c qxy * x * y + qyz * y * z + qzx * z * x + c qxx * x**2 + qyy * y**2 + qzz * z**2 = 0, c c The method is to try various directions from point "a", starting c with direction "b", and the directions to any extrema and axial c intercept points on the quadric surface, and the x, y and z c directions, then trying nran directions randomly sampled from c an isotropic distribution, and then trying up to 10 * nran c (or until the result converges, based on tol) directions c randomly sampled from a cosine-power distribution directed c toward the last best direction, with the power geometrically c increasing from 1.0 to 1.e10. c c The final vector from point "a" to point "c" is vector c "d" = (dx, dy, dz) = (cx - ax, cy - ay, cz - az). c c The vector normal to the quadric surface at point "c" is c "sn" = (snx, sny, snz) = grad f: c snx = qx + 2.0 * qxx * cx + qxy * cy + qzx * cz, c sny = qy + qxy * cx + 2.0 * qyy * cy + qyz * cz, c snz = qz + qzx * cx + qyz * cy + 2.0 * qzz * cz, c and the cosine of the angle between vector "d" and the normal c vector "sn" is costh, which is -1 or 1 if the method converges c (nprox = 1). c c Any result less than the estimated error in its calculation, c based on tol, will be truncated to zero. c c Flag nerr indicates any input error. c c Input: ax, ay, az, bx, by, bz, qc, qx, qy, qz, qxy, qyz, qzx, c qxx, qyy, qzz, nran, tol. c c Output: nprox, cx, cy, cz, dist, dx, dy, dz, snx, sny, snz, costh, nerr. c c Calls: aptrkis, aptscat, aptscav, aptvusz c c c Glossary: c c ax,ay,az Input The x, y, z coordinates of point "a". c c bx,by,bz Input The x, y, z components of the first trial vector "b". c c costh Output The cosine of the angle between vector "d" (from point c "a" to point "c") and the vector normal to the c quadric surface at point "c". If nprox = 1, costh c should be either -1.0 or 1.0. c c cx,cy,cz Output The x, y, z coordinates of point "c" on the quadric c surface, found nearest to point "a". If nprox = 1, c the proximal point, to within tolerance limit tol. c c dx,dy,dz Output The x, y, z components of the vector from point "a" to c point "c". c c dist Output The distance from point "a" to point "c". c If nprox = 1, the proximal distance, to within c tolerance limit tol. c c nerr Output Indicates an input error, if not 0. c 1 if vector "b" is too small. c 2 if nran is not positive. c c nprox Output The number of exact proximal points "c" and proximal c distances dist found. c -1 if none, and none of the trial directions from c point "a" intersect the quadric surface. c The coefficients qc, ..., qzz may be bad. c 0 if none, but the surface point "c" in direction "d" c from point "a" produced the smallest distance dist c of any directions tried. c 1 if an exact proximal point "c" was found, within c tolerance limit tol. c c q.. Input Coefficients of the implicit equation of a quadric c surface in xyz space (qc, qx, qy, qz, qxy, qyz, qzx, c qxx, qyy, qzz). c c snx,... Output The x, y, z components of the vector normal to the c quadric surface at point "c". c c tol Input Numerical tolerance limit. c On Cray computers, recommend 1.e-5 to 1.e-11. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccend. c.... Local variables. common /laptqprr/ ccx(4) ! X coordinate of surface point. common /laptqprr/ ccy(4) ! Y coordinate of surface point. common /laptqprr/ ccz(4) ! Z coordinate of surface point. common /laptqprr/ ddist(4) ! Distance to quadric surface. common /laptqprr/ aqp(25) ! Type of surface point. character*8 aqp common /laptqprr/ qpx(25) ! Coordinate x of surface point. common /laptqprr/ qpy(25) ! Coordinate y of surface point. common /laptqprr/ qpz(25) ! Coordinate z of surface point. cbugc***DEBUG begins. cbug 9901 format (/ 'aptqprr finding (nran=',i8',) the distance from', cbug & ' point "a"' / cbug & ' to the quadric surface with the coefficients qc-qzz,' / cbug & ' and initial trial direction "b":' / cbug & ' ax, ay, az =',1p3e22.14 / cbug & ' bx, by, bz =',1p3e22.14 / 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) nran, ax, ay, az, bx, by, bz, cbug & qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz cbugc***DEBUG ends. c.... Initialize. nerr = 0 nprox = 0 ndecr = 0 ! Total number of improved guesses. distmin = 1.e99 ! Minimum distance found. dist = -1.e99 cx = -1.e99 cy = -1.e99 cz = -1.e99 dx = -1.e99 dy = -1.e99 dz = -1.e99 snx = -1.e99 sny = -1.e99 snz = -1.e99 costh = -1.e99 c.... Test for errors. call aptvusz (bx, by, bz, tol, ubx, uby, ubz, blen) if (blen .le. tol) then nerr = 1 nprox = -1 go to 710 endif if (nran .lt. 1) then nerr = 2 nprox = -1 go to 710 endif c=======================================================================******** c.... Find any extrema and axis intercepts on the quadric surface. call aptqupt (qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz, & tol, aqp, qpx, qpy, qpz, nqp, nerr) c... Add specified direction "b", and x, y and z directions. nqp = nqp + 1 qpx(nqp) = ax + bx qpy(nqp) = ay + by qpz(nqp) = az + bz nqp = nqp + 1 qpx(nqp) = ax + 1.0 qpy(nqp) = ay + 0.0 qpz(nqp) = az + 0.0 nqp = nqp + 1 qpx(nqp) = ax + 0.0 qpy(nqp) = ay + 1.0 qpz(nqp) = az + 0.0 nqp = nqp + 1 qpx(nqp) = ax + 0.0 qpy(nqp) = ay + 0.0 qpz(nqp) = az + 1.0 c=======================================================================******** c.... Loop over input vector and directions from point "a" to surface points, c.... then loop over nran random directions. nranp = nran + nqp do 210 nd = 1, nranp ! Loop over random directions. if (nd .le. nqp) then ! Use directions to surface points, and "b". call aptvdis (ax, ay, az, qpx(nd), qpy(nd), qpz(nd), 1, tol, & vx, vy, vz, vlen, nerra) if (vlen .eq. 0.0) go to 210 vx = vx / vlen vy = vy / vlen vz = vz / vlen else ! Use a random direction. call aptscat (1, vx, vy, vz, nerra) endif c.... Find the intersection of the track from point "a" in direction "v". dintmn = -distmin dintmx = distmin call aptrkis (ax, ay, az, vx, vy, vz, & qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz, & dintmn, dintmx, 1, tol, & nint, sx, sy, sz, distest, & ssnx, ssny, ssnz, cosn, nerra) if (nint .eq. 0) go to 210 if (abs (distest) .le. distmin) then if (distest .lt. 0.0) then distest = -distest cosn = -cosn vx = -vx vy = -vy vz = -vz endif ndecr = ndecr + 1 cx = sx cy = sy cz = sz distmin = abs (distest) dist = distest dx = vx * distest dy = vy * distest dz = vz * distest snx = ssnx sny = ssny snz = ssnz costh = cosn cbugc***DEBUG begins. cbug 9902 format (/ 'Proximal x,y,z ?',1p3e20.12 / cbug & 'Proximal dx,dy,dz?',1p3e20.12 / cbug & ' Approx distance?',1pe20.12,' Test',i8 ) cbug 9903 format ( ' Cosine of angle ',1pe20.12,' (from normal)', cbug & 5x,1pe20.12 ) cbug power = 0.0 cbug write ( 3, 9902) cx, cy, cz, dx, dy, dz, dist, ndecr cbug write ( 3, 9903) costh, power cbugc***DEBUG ends. if (distmin .eq. 0.0) then nprox = 1 costh = 1.0 go to 710 endif endif 210 continue ! End of loop over random directions. if (ndecr .eq. 0) then ! No intersections. nprox = -1 go to 710 endif delcos = 1.0 - abs (costh) if (delcos .le. tol) then ! Found an exact proximal point. nprox = 1 go to 710 endif c=======================================================================******** c.... Search in cosine-power distribution around the best direction. numtot = 10 * nran power = 1.0 powmult = 10.0**(1.0 / nran) do 250 nd = 1, numtot ! Loop over cosine-power distribution. call aptvusz (dx, dy, dz, tol, & ugoodx, ugoody, ugoodz, dlen) call aptscav (1, power, ugoodx, ugoody, ugoodz, & vx, vy, vz, nerra) power = power * powmult c.... Find the intersection of the track from point P in direction V. dintmn = -distmin dintmx = distmin call aptrkis (ax, ay, az, vx, vy, vz, & qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz, & dintmn, dintmx, 1, tol, & nint, sx, sy, sz, distest, & ssnx, ssny, ssnz, cosn, nerra) if (nint .eq. 0) go to 250 ! No intersection. if (abs (distest) .le. distmin) then if (distest .lt. 0.0) then distest = -distest cosn = -cosn vx = -vx vy = -vy vz = -vz endif ndecr = ndecr + 1 cx = sx cy = sy cz = sz distmin = abs (distest) dist = distest dx = vx * distest dy = vy * distest dz = vz * distest snx = ssnx sny = ssny snz = ssnz costh = cosn cbugc***DEBUG begins. cbug write ( 3, 9902) cx, cy, cz, dx, dy, dz, dist, ndecr cbug write ( 3, 9903) costh, power cbugc***DEBUG ends. deldist = abs (abs (distest) - distmin) delcos = 1.0 - abs (costh) c.... Test for final convergence. if ((deldist .le. tol * distmin) .and. & (delcos .le. tol)) then nprox = 1 go to 710 endif endif ! Improved guess. 250 continue ! End of loop over cosine-power distr. c=======================================================================******** 710 continue cbugc***DEBUG begins. cbug 9904 format (/ 'aptqprr results: nprox=',i2,' ndecr=',i8,' nerr=',i3 / cbug & ' dist =',1pe22.14 / cbug & ' cx, cy, cz =',1p3e22.14 / cbug & ' dx, dy, dz =',1p3e22.14 / cbug & ' snx,sny,snz=',1p3e22.14 / cbug & ' costh =',1pe22.14 ) cbug write ( 3, 9904) nprox, ndecr, nerr, dist, cx, cy, cz, cbug & dx, dy, dz, snx, sny, snz, costh cbugc***DEBUG ends. return c.... End of subroutine aptqprr. (+1 line.) end UCRL-WEB-209832