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