subroutine aptcnis (aw, ar, bw, br, np, tol, qc, qu, qv, qw, & quv, qvw, qwu, quu, qvv, qww, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTCNIS c c call aptcnis (aw, ar, bw, br, np, tol, qc, qu, qv, qw, c & quv, qvw, qwu, quu, qvv, qww, nerr) c c Version: aptcnis Updated 1990 December 12 14:30. c aptcnis Originated 1990 December 7 11:50. 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 c a cone symmetric around the w axis, that passes through the c two distinct points a = (aw, ar) and b = (bw, br), c where ar = sqrt (au**2 + av**2) and br = sqrt (bu**2 + bv**2). c Points "a" and "b" are on the same side of the vertex of the c cone if ar and br have the same sign, otherwise are on opposite c sides. The coordinates (u,v,w) are along the three major c orthogonal axes, and may represent (x,y,z), (y,z,x) or (z,x,y). c c The equation for the conical surface is as follows: c c f(u,v,w) = qc + qu * u + qv * v + qw * w + c quv * u * v + qvw * v * w + qwu * w * u + c quu * u**2 + qvv * v**2 + qww * w**2 = 0 c c where qu = qv = quv = qvw = qwu = 0, quu = qvv, and c qww has the opposite sign of quu. c c Degenerate case: c c If points "a" and "b" coincide, within the limits determined c by tol, all coefficients will be zero. c c Other special cases: c c If aw = bw, the surface will be a w plane. c c If ar = br, the surface will be a right circular cylinder. c c If ar = br = 0, the surface will be a line on the w axis. c c Note: Use apttris and aptrois to translate and rotate the cone c into any desired position and orientation. c Use aptcois to find the equation of a cone on a given axis, c with a given vertex angle and vertex axial coordinate. c Use aptcyis to find the equation of a cylinder symmetric c Around a major axis. c Use aptelis to find the equation of an ellipsoid symmetric c around a major axis. c c Input: aw, ar, bw, br, np, tol. c c Output: qc, qu, qv, qw, quv, qvw, qwu, quu, qvv, qww, nerr. c c Calls: aptvdic, aptvaxc c c c Glossary: c c ar, aw Input Coordinates r and w of first point in conical surface. c Size np. c If points "a" and "b" are on the same side of the c vertex of the cone, the signs of ar and br must be c the same. Otherwise, the signs must be different. c c br, bw Input Coordinates r and w of second point in conical surface. c Size np. c If points "a" and "b" are on the same side of the c vertex of the cone, the signs of ar and br must be c the same. Otherwise, the signs must be different. 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 aw, ar, bw, br, qc, qu, qv, qw, c quv, qvw, qwu, quu, qvv, qww. c c qc Output Constant term in quadric cone equation. Size np. c Zero only if the vertex of the cone is at the origin. c qc = -(aw * br - ar * bw)**2. c c qu,qv,qw Output Coefficients of u, v, w in quadric cone equation. c Size np. c qu = qv = 0. c qw = 2.0 * (aw * br - ar * bw) * (br - ar). c c q.. Output Coefficients of second-order terms in quadric c cone equation: quv, qvw, qwu, quu, qvv, qww. c Coefficients of u*v, v*w, w*u, u**2, v**2, w**2, c respectively. Size np. c quv = qvw = qwu = 0. c quu = qvv = (bw - aw)**2. c qww = -(br - ar)**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---- Coordinate r of point "a". dimension ar (1) c---- Coordinate w of point "a". dimension aw (1) c---- Coordinate r of point "b". dimension br (1) c---- Coordinate w of point "b". dimension bw (1) c---- Coefficient of term 1.0. dimension qc (1) c---- Coefficient of term u. dimension qu (1) c---- Coefficient of term u**2. dimension quu (1) c---- Coefficient of term u * v. dimension quv (1) c---- Coefficient of term v. dimension qv (1) c---- Coefficient of term v**2. dimension qvv (1) c---- Coefficient of term v * w. dimension qvw (1) c---- Coefficient of term w. dimension qw (1) c---- Coefficient of term w * u. dimension qwu (1) c---- Coefficient of term w**2. dimension qww (1) c.... Local variables. c---- Distance from point "a" to point "b". common /laptcnis/ abd (64) c---- Cross product of vectors "a" and "b". common /laptcnis/ abm (64) c---- Difference br - ar. common /laptcnis/ abr (64) c---- Difference bw - aw. common /laptcnis/ abw (64) c---- Index in local array. common /laptcnis/ n c---- First index of subset of data. common /laptcnis/ n1 c---- Last index of subset of data. common /laptcnis/ n2 c---- Index in external array. common /laptcnis/ nn c---- Size of current subset of data. common /laptcnis/ ns cbugc***DEBUG begins. cbug 9901 format (/ 'aptcnis finding the quadric coefficients for a cone.', cbug & ' tol=',1pe13.5 / cbug & ' Cone passes through the points:') cbug 9902 format ( cbug & i3,' aw,ar=',1p2e22.14 / cbug & ' bw,br=',1p2e22.14) cbug write ( 3, 9901) tol cbug write ( 3, 9902) (n, aw(n), ar(n), bw(n), br(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 (/ "aptcnis fatal. Non-positive np =",i5) cbug write ( 3, 9903) np cbugc***DEBUG ends. go to 210 endif c.... Set up the indices of the first subset of data. n1 = 1 n2 = min (np, 64) 110 ns = n2 - n1 + 1 c.... Find the 2-D vector from point "a" to point "b". call aptvdic (aw(n1), ar(n1), bw(n1), br(n1), ns, tol, & abw, abr, abd, nerr) c.... Find the cross product of the 2-D vectors "a" and "b". call aptvaxc (aw(n1), ar(n1), bw(n1), br(n1), ns, tol, abm, nerr) c.... Find the coefficients of the equation of the cone. c---- Loop over subset of data. do 120 n = 1, ns nn = n + n1 - 1 qc(nn) = -abm(n)**2 qu(nn) = 0.0 qv(nn) = 0.0 qw(nn) = 2.0 * abm(n) * abr(n) quv(nn) = 0.0 qvw(nn) = 0.0 qwu(nn) = 0.0 quu(nn) = abw(n)**2 qvv(nn) = quu(nn) qww(nn) = -abr(n)**2 c---- End of loop over subset of data. 120 continue c.... See if all data subsets are done. c---- Do another subset of data. if (n2 .lt. np) then n1 = n2 + 1 n2 = min (np, n1 + 63) go to 110 endif cbugc***DEBUG begins. cbug 9904 format (/ 'aptcnis results:' / cbug & (i3,' qc= ',1pe22.14 / cbug & ' qu,qv,qw= ',1p3e22.14 / cbug & ' quv,qvw,qwu=',1p3e22.14 / cbug & ' quu,qvv,qww=',1p3e22.14)) cbug write ( 3, 9904) (n, qc(n), qu(n), qv(n), qw(n), cbug & quv(n), qvw(n), qwu(n), quu(n), qvv(n), qww(n), cbug & n = 1, np) cbugc***DEBUG ends. 210 return c.... End of subroutine aptcnis. (+1 line.) end UCRL-WEB-209832