subroutine aptqexv (dx, dy, dz, qc, qx, qy, qz, qxy, qyz, qzx, & qxx, qyy, qzz, tol, & nextr, ex, ey, ez, adir, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTQEXV c c call aptqexv (dx, dy, dz, qc, qx, qy, qz, qxy, qyz, qzx, c & qxx, qyy, qzz, tol, c & nextr, ex, ey, ez, adir, nerr) c c Version: aptqexv Updated 1994 October 4 16:50. c aptqexv Originated 1994 October 4 16:50. c c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123. c c c Purpose: For the quadric surface specified by the coordinates qc to qzz, c to find any extremities e = (ex, ey, ez) in the direction of the c vector d = (dx, dy, dz), and the type adir of extremities. c The possible types of extrema are minima and maxima (points, c lines or planes), saddle points, constant planes, intersection c lines, and isolated points. c c The implicit equation of the quadric curve is: c c F(x,y,z) = qc + qx * x + qy * y + c qxy * x * y + qyz * y * z + qzx * z * x + c qxx * x**2 + qyy * y**2 + qzz * z**2 = 0 . c c Results may be unpredictable if the coefficients of the quadric c surface differ by many orders of magnitude. c c Method: See APT subroutine aptextr. c c Input: dx, dy, dz, qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz, tol. c c Output: nextr, ex, ey, ez, adir, nerr. c c Calls: aptmopv, aptqext, aptrois, aptrotv c c c Glossary: c c adir Output An ASCII 8-character word, the type of extreme point: c "Minimum", "Maximum", "Plane", "IntPlane", c "Saddle", "Point". c If none, return "unknown". Size = 2. c c dx,dy,dz Input Components x, y and z of the direction in which any c extreme points are to be found. c c ex,ey,ez Output The coordinates of any extreme points. Size = 2. c c nerr Output Indicates an input error, if not zero: c 1 if all coefficients, except possible qc, are zero, c or the quadric surface type is 17, 19, 20, 23 or 24. c 2 if the magnitude of vector "d" is too small. c c nextr Output Number of extreme or intersection points found. c c qc Input Constant term in the equation of the quadric curve. c c qx,qy,qz Input Coefficients of x, y and z, respectively, in the c equation of the quadric curve. c c qxx Input Coefficient of x**2 in the equation of the quadric c curve. c c qxy Input Coefficient of x * y in the equation of the quadric c curve. c c qyy Input Coefficient of y**2 in the equation of the quadric c curve. c c qyz Input Coefficient of y * z in the equation of the quadric c curve. c c qzx Input Coefficient of z * x in the equation of the quadric c curve. c c qzz Input Coefficient of z**2 in the equation of the quadric c curve. c c tol Input Numerical tolerance limit. c On Cray computers, recommend 1.e-11. c On 32-bit computers, recommend 1.e-6. c c History: c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccend. c.... Dimensioned arguments. dimension adir(2) ! Type of extreme point. character*8 adir ! Type of extreme point. dimension ex(2) ! Coordinate x of an extreme point. dimension ey(2) ! Coordinate y of an extreme point. dimension ez(2) ! Coordinate z of an extreme point. c.... Local variables. common /laptqexv/ abdir(6) ! Type of extreme point. character*8 abdir ! Type of extreme point. common /laptqexv/ rotm(3,3) ! Rotation matrix, vector "d" to z axis. common /laptqexv/ bx(6) ! Coordinate x of an extreme point. common /laptqexv/ by(6) ! Coordinate y of an extreme point. common /laptqexv/ bz(6) ! Coordinate z of an extreme point. cbugc***DEBUG begins. cbug 9901 format (/ 'aptqexv finding extrema of quadric in "d" dir.' / cbug & ' tol=',1pe13.5 / cbug & ' dx,dy,dz= ',1p3e22.14 ) cbug 9902 format ( 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, dx, dy, dz cbug write ( 3, 9902) qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz cbugc***DEBUG ends. c.... Initialize. nerr = 0 do 110 n = 1, 2 ex(n) = -1.e99 ey(n) = -1.e99 ez(n) = -1.e99 adir(n) = 'unknown' 110 continue c.... Find any extrema in the quadric in the direction (dx, dy, dz). nextr = 0 c.... Find the rotation operator to rotate the axis vector (dx, dy, dz) to the c.... z axis. call aptrotv (dx, dy, dz, 0., 0., 1., tol, rotm, nerra) cbugc***DEBUG begins. cbug 9990 format (' op11,op12,op13 ',1p3e20.12 / cbug & ' op21,op22,op23 ',1p3e20.12 / cbug & ' op31,op32,op33 ',1p3e20.12) cbug write ( 3, 9990) ((rotm(i,j), j = 1, 3), i = 1, 3) cbugc***DEBUG ends. if (nerra .ne. 0) then nerr = 2 go to 999 endif c.... Rotate the quadric surface. ac = qc ax = qx ay = qy az = qz axy = qxy ayz = qyz azx = qzx axx = qxx ayy = qyy azz = qzz cbugc***DEBUG begins. cbug 9960 format (' ac ',1pe20.12 ) cbug 9980 format (' ax, ay, az ',1p3e20.12 ) cbug 9985 format (' axy, ayz, azx ',1p3e20.12 ) cbug 9991 format (' axx, ayy, azz ',1p3e20.12 ) cbug write ( 3, 9960) ac cbug write ( 3, 9980) ax, ay, az cbug write ( 3, 9985) axy, ayz, azx cbug write ( 3, 9991) axx, ayy, azz cbugc***DEBUG ends. call aptrois (rotm, 0, 0., 0., 0., ac, ax, ay, az, & axy, ayz, azx, axx, ayy, azz, tol, nerra) cbugc***DEBUG begins. cbug write ( 3, 9960) ac cbug write ( 3, 9980) ax, ay, az cbug write ( 3, 9985) axy, ayz, azx cbug write ( 3, 9991) axx, ayy, azz cbugc***DEBUG ends. c.... Find the extrema in the z direction. call aptqext (ac, ax, ay, az, axy, ayz, azx, axx, ayy, azz, & tol, ntyps, cenx, ceny, cenz, nmm, bx, by, bz, & abdir, nerra) cbugc***DEBUG begins. cbug 9922 format ('nmm=',i2) cbug 9923 format ('n=',i2,' abdir="',a8,'"') cbug write ( 3, 9922) nmm cbug if (nmm .gt. 0) then cbug write ( 3, 9923) (n, abdir(n), n = 1, nmm) cbug endif cbugc***DEBUG ends. if (nerra .ne. 0) then nerr = 1 go to 999 endif c.... Rotate any extrema back to the original direction. if (nmm .gt. 0) then call aptmopv (rotm, 1, 0., 0., 0., bx, by, bz, nmm, tol, nerra) endif c.... Find any 'Min z' or 'Max z', relabel 'Minimum' or 'Maximum'. c.... Find any 'Const z' or 'Inter z', relabel 'Plane' or 'IntPlane'. c.... Find any 'Saddle z' or 'Point z', relabel 'Saddle' or 'Point'. do 150 n = 1, nmm ! Loop over extreme points. if (abdir(n) .eq. 'Min z') then nextr = nextr + 1 adir(nextr) = 'Minimum' ex(nextr) = bx(n) ey(nextr) = by(n) ez(nextr) = bz(n) elseif (abdir(n) .eq. 'Max z') then nextr = nextr + 1 adir(nextr) = 'Maximum' ex(nextr) = bx(n) ey(nextr) = by(n) ez(nextr) = bz(n) elseif (abdir(n) .eq. 'Const z') then nextr = nextr + 1 adir(nextr) = 'Plane' ex(nextr) = bx(n) ey(nextr) = by(n) ez(nextr) = bz(n) elseif (abdir(n) .eq. 'Inter z') then nextr = nextr + 1 adir(nextr) = 'IntPlane' ex(nextr) = bx(n) ey(nextr) = by(n) ez(nextr) = bz(n) elseif (abdir(n) .eq. 'Saddle z') then nextr = nextr + 1 adir(nextr) = 'Saddle' ex(nextr) = bx(n) ey(nextr) = by(n) ez(nextr) = bz(n) elseif (abdir(n) .eq. 'Point z') then nextr = nextr + 1 adir(nextr) = 'Point' ex(nextr) = bx(n) ey(nextr) = by(n) ez(nextr) = bz(n) endif 150 continue ! End of loop over extreme points. 999 continue cbugc***DEBUG begins. cbug 9910 format (/ 'aptqexv results: nextr=',i2,' nerr=',i2) cbug 9912 format (a8,' xyz=',1p3e22.14) cbug write ( 3, 9910) nextr, nerr cbug do 991 n = 1, 2 cbug write ( 3, 9912) adir(n), ex(n), ey(n), ez(n) cbug 991 continue cbugc***DEBUG ends. return c.... End of subroutine aptqexv. (+1 line) end UCRL-WEB-209832