subroutine aptcois (kang, theta, aw, np, tol, qc, qu, qv, qw, & quv, qvw, qwu, quu, qvv, qww, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTCOIS c c call aptcois (kang, theta, aw, np, tol, qc, qu, qv, qw, c & quv, qvw, qwu, quu, qvv, qww, nerr) c c Version: aptcois Updated 1990 December 12 14:30. c aptcois Originated 1990 December 7 11:40. 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, with the vertex at w = aw, c and the angle theta between the w axis and the conical surface. c The coordinates (u,v,w) are along the three major orthogonal c 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 Special cases: c c The result will be a line if theta is zero plus or minus c any integer multiple of 180 degrees (pi radians). c c The result will be a plane if theta is 90 degrees c (pi/2 radians) plus or minus any integer multiple of c 180 degrees (pi radians). c c Note: Use apttris and aptrois to translate and rotate the cone c into any desired position and orientation. c Use aptcnis to find the equation of a cone on a given axis, c and passing through two given points. 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: kang, theta, aw, np, tol. c c Output: qc, qu, qv, qw, quv, qvw, qwu, quu, qvv, qww, nerr. c c Glossary: c c aw Input Coordinate w of the vertex of the cone on the w axis. c Size np. c c kang Input Indicates theta units are degrees (0) or radians (1). c c nerr Output Indicates an input error, if not zero: c 1 if np is not positive. c 2 if kang is not 0 or 1. c c np Input Size of arrays theta, aw, 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 cone passes through the origin. c qc = -aw**2 * (sin (theta))**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 * (sin (theta))**2. 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 = (cos (theta))**2. c qww = -(sin (theta))**2. c c theta Input Angle from w axis to conical surface. Units are c degrees (kang = 0) or radians (kang = 1). Size np. 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. Size np. c---- Coordinate w of vertex of cone. dimension aw (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---- Angle from w axis to cone. dimension theta (1) c.... Local variables. c---- Vertex angle in radians. common /laptcois/ angle c---- Square of cos (theta). common /laptcois/ cthsq (64) c---- Index in local array. common /laptcois/ n c---- First index of subset of data. common /laptcois/ n1 c---- Last index of subset of data. common /laptcois/ n2 c---- Index in external array. common /laptcois/ nn c---- Size of current subset of data. common /laptcois/ ns c---- Mathematical constant, pi. common /laptcois/ pi c---- Square of sin (theta). common /laptcois/ sthsq (64) cbugc***DEBUG begins. cbug 9901 format (/ 'aptcois finding the quadric coefficients for a cone.', cbug & ' tol=',1pe13.5,' kang=',i2) cbug 9902 format (i3,' theta=',1pe22.14,' aw=',1pe22.14) cbug write ( 3, 9901) tol, kang cbug write ( 3, 9902) (n, theta(n), aw(n), n = 1, np) cbugc***DEBUG ends. c.... Initialize. c---- Mathematical constant, pi. pi = 3.14159265358979323 nerr = 0 c.... Test for input errors. if (np .le. 0) then nerr = 1 cbugc***DEBUG begins. cbug 9903 format (/ "aptcois fatal. Non-positive np =",i5) cbug write ( 3, 9903) np cbugc***DEBUG ends. go to 210 endif if ((kang .lt. 0) .or. (kang .gt. 1)) then nerr = 2 cbugc***DEBUG begins. cbug write ( 3, '(/ "aptcois fatal. Bad kang=",i3)') kang 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 vertex angle functions. c---- Angles are in degrees. if (kang .eq. 0) then c---- Loop over subset of data. do 120 n = 1, ns nn = n + n1 - 1 angle = theta(nn) * pi / 180.0 cthsq(n) = (cos (angle))**2 sthsq(n) = (sin (angle))**2 c---- End of loop over subset of data. 120 continue c---- Angles are in radians. else c---- Loop over subset of data. do 130 n = 1, ns nn = n + n1 - 1 angle = theta(nn) cthsq(n) = (cos (angle))**2 sthsq(n) = (sin (angle))**2 c---- End of loop over subset of data. 130 continue endif c.... See if truncation tests should be made. c---- Test for truncation to zero. if (tol .gt. 0.0) then c---- Loop over subset of data. do 140 n = 1, ns if (abs (sthsq(n)) .lt. tol) then sthsq(n) = 0.0 endif if (abs (cthsq(n)) .lt. tol) then cthsq(n) = 0.0 endif c---- End of loop over subset of data. 140 continue c---- Tested tol. endif c.... Find the coefficients of the equation of the cone. c---- Loop over subset of data. do 150 n = 1, ns nn = n + n1 - 1 qc(nn) = -aw(nn)**2 * sthsq(n) qu(nn) = 0.0 qv(nn) = 0.0 qw(nn) = 2.0 * aw(nn) * sthsq(n) quv(nn) = 0.0 qvw(nn) = 0.0 qwu(nn) = 0.0 quu(nn) = cthsq(n) qvv(nn) = quu(nn) qww(nn) = -sthsq(n) c---- End of loop over subset of data. 150 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 9905 format (/ 'aptcois 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, 9905) (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 aptcois. (+1 line.) end UCRL-WEB-209832