subroutine aptesis (ax, ay, az, np, tol, qc, qx, qy, qz, & qxy, qyz, qzx, qxx, qyy, qzz, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTESIS c c call aptesis (ax, ay, az, np, tol, qc, qx, qy, qz, c & qxy, qyz, qzx, qxx, qyy, qzz, nerr) c c Version: aptesis Updated 1990 December 13 10:30. c aptesis Originated 1990 December 13 10:30. c c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123. c c c Purpose: For each of np sets of input data, to find the equation of the c ellipsoid centered at the origin, with semiaxes of ax, ay c and az on the x, y and z axes, respectively. c c The equation for the ellipsoid is as follows: c c f(u,v,w) = qc + qx * u + qy * v + qz * w + c qxy * u * v + qyz * v * w + qzx * w * u + c qxx * u**2 + qyy * v**2 + qzz * w**2 = 0 c c where qx = qy = qz = qxy = qyz = qzx = 0, c qc = - ax**2 * ay**2 * az**2, qxx = ay**2 * az**2, c qyy = az**2 * ax**2, qzz = ax**2 * ay**2. c c General case: c c If ax, ay and az are all non-zero, the surface will be c an ellipsoid with semiaxes ax, ay, az on the x, y and c z axes, respectively. c c Indeterminate cases: c c If none, or only one of ax, ay or az is nonzero, as determined c by tol, the surface is indeterminate (or a point at the origin), c and all coefficients will be zero. c c Other singular cases: c c If ax = ay = az, the surface will be a sphere of radius ax. c c If any two of ax, ay or az are equal, the surface will be an c ellipsoid symmetric around the axis of the non-equal semiaxis. c c If just one of the values ax, ay or az is zero, as determined c by tol, the surface will be an x, y or z plane through the c origin, respectively. c c Note: Use apttris and aptrois to translate and rotate the ellipsoid c into any desired position and orientation. c Use aptelis to find the equation of an ellipsoid symmetric c around a major axis and passing through two given points. c c Input: ax, ay, az, np, tol. c c Output: qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz, nerr. c c Glossary: c c ax,ay,az Input Semiaxes of the ellipsoid. Aligned in the directions c of the x, y and z axis, respectively. c c nerr Output Indicates an input error, if not zero: c 1 if np is not positive. c c np Input Size of arrays ax, ay, az, qc, qx, qy, qz, c qxy, qyz, qzx, qxx, qyy, qzz. c c qc Output Constant term in the equation of the ellipsoid. c Size np. qc = -ax**2 * ay**2 * az**2. c c qx,qy,qz Output Coefficients of u, v, w in the equation of the c ellipsoid. Size np. qx = qy = qz = 0. c c q.. Output Coefficients of second-order terms in the equation of c the ellipsoid: qxy, qyz, qzx, qxx, qyy, qzz. c Coefficients of u*v, v*w, w*u, u**2, v**2, w**2, c respectively. Size np. qxy = qyz = qzx = 0. c qxx = ay**2 * az**2, qyy = az**2 * ax**2, c qzz = ax**2 * ay**2. c c tol Input Numerical tolerance limit. c On Cray computers, recommend 1.e-5 to 1.e-11. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccend. c.... Dimensioned arguments. c---- Semiaxis along the x axis. dimension ax (1) c---- Semiaxis along the y axis. dimension ay (1) c---- Semiaxis along the z axis. dimension az (1) c---- Coefficient of term 1.0. dimension qc (1) c---- Coefficient of term u. dimension qx (1) c---- Coefficient of term u**2. dimension qxx (1) c---- Coefficient of term u * v. dimension qxy (1) c---- Coefficient of term v. dimension qy (1) c---- Coefficient of term v**2. dimension qyy (1) c---- Coefficient of term v * w. dimension qyz (1) c---- Coefficient of term w. dimension qz (1) c---- Coefficient of term w * u. dimension qzx (1) c---- Coefficient of term w**2. dimension qzz (1) c.... Local variables. c---- Estimated error in ax, ay, az. common /laptesis/ aerr c---- Index in external array. common /laptesis/ n cbugc***DEBUG begins. cbug 9901 format (/ 'aptesis finding the equation of a ellipsoid.', cbug & ' tol=',1pe13.5) cbug 9902 format (i3,' ax,ay,az=',1p3e22.14) cbug write ( 3, 9901) tol cbug write ( 3, 9902) (n, ax(n), ay(n), az(n), n = 1, np) cbugc***DEBUG ends. c.... Initialize. nerr = 0 c.... Test for input errors. if (np .le. 0) then nerr = 1 cbugc***DEBUG begins. cbug 9903 format (/ "aptesis fatal. Non-positive np =",i5) cbug write ( 3, 9903) np cbugc***DEBUG ends. go to 210 endif c.... Find the coefficients of the equation of the ellipsoid. c---- Loop over data. do 120 n = 1, np qc(n) = -ax(n)**2 * ay(n)**2 * az(n)**2 qx(n) = 0.0 qy(n) = 0.0 qz(n) = 0.0 qxy(n) = 0.0 qyz(n) = 0.0 qzx(n) = 0.0 qxx(n) = ay(n)**2 * az(n)**2 qyy(n) = az(n)**2 * ax(n)**2 qzz(n) = ax(n)**2 * ay(n)**2 c---- End of loop over data. 120 continue c.... See if the results should be tested for truncation error. c---- Test for truncation to zero. if (tol .gt. 0.0) then c---- Loop over data. do 130 n = 1, np aerr = tol * (abs (ax(n)) + abs (ay(n)) + abs (az(n))) if (abs (ax(n)) .lt. aerr) then qc(n) = 0.0 qyy(n) = 0.0 qzz(n) = 0.0 endif if (abs (ay(n)) .lt. aerr) then qc(n) = 0.0 qxx(n) = 0.0 qzz(n) = 0.0 endif if (abs (az(n)) .lt. aerr) then qc(n) = 0.0 qxx(n) = 0.0 qyy(n) = 0.0 endif c---- End of loop over data. 130 continue c---- Tested tol. endif cbugc***DEBUG begins. cbug 9904 format (/ 'aptesis results:' / cbug & (i3,' qc= ',1pe22.14 / cbug & ' qx,qy,qz= ',1p3e22.14 / cbug & ' qxy,qyz,qzx=',1p3e22.14 / cbug & ' qxx,qyy,qzz=',1p3e22.14)) cbug write ( 3, 9904) (n, qc(n), qx(n), qy(n), qz(n), cbug & qxy(n), qyz(n), qzx(n), qxx(n), qyy(n), qzz(n), cbug & n = 1, np) cbugc***DEBUG ends. 210 return c.... End of subroutine aptesis. (+1 line.) end UCRL-WEB-209832