subroutine aptcris (ax, ay, az, qc, qx, qy, qz, qxy, qyz, qzx,
& qxx, qyy, qzz, tol, bx, by, bz,
& nrad, rcurve, cx, cy, cz, ux, uy, uz, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTCRIS
c
c call aptcris (ax, ay, az, qc, qx, qy, qz, qxy, qyz, qzx,
c & qxx, qyy, qzz, tol, bx, by, bz,
c & nrad, rcurve, cx, cy, cz, ux, uy, uz, nerr)
c
c Version: aptcris Updated 1997 December 3 13:50.
c aptcris Originated 1997 October 31 14:00.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c Source: In directory ~edwards/work/apt/src on the DEC Compass Cluster.
c
c Purpose: To find, at point a = (ax, ay, az) on the quadric surface
c specified by the coefficients qx to qzz (with qc adjusted to
c make the quadric surface pass through point "a"), the normal
c vector b = (bx, by, bz), the number nrad of tangent directions
c in the surface, with extreme values of the radius of curvature
c rcurve, the corresponding centers of curvature c = (cx, cy, cz),
c and the tangent unit vectors u = (ux, uy, uz).
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 Flag nerr indicates any error.
c
c Method:
c
c The normal vector "b" at point "a" has the components:
c
c bx = qx + 2.0 * qxx * ax + qxy * ay + qzx * az
c by = qy + 2.0 * qyy * ay + qyz * az + qxy * ax
c bz = qz + 2.0 * qzz * az + qzx * ax + qyz * ay
c
c If the system is rotated to point vector "b" in the z direction,
c then all tangent unit vectors through point "a" will satisfy
c u = (cos(t), sin(t), 0),
c where t is an arbitrary angle, and
c rcurve = -|b| / (2*W),
c where |b| is the magnitude of vector "b", and
c W = qxx*cos(t)**2 + qxy*cos(t)*sin(t) + qyy*sin(t)**2.
c The extreme values of rcurve occur at W = 0, if there are any
c real values of t that give that result, and at the extreme
c values of W, found by differentiating the equation for W with
c respect to t, setting the result to zero, and solving for t:
c (qyy - qxx) * cos(2*t) + qxy * sin(2*t) = 0
c For each solution t, the values of W, rcurve and U may be found.
c The center of curvature is at c = a - b / (2*W).
c
c Input: ax, ay, az, qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz, tol.
c
c Output: rcurve, cx, cy, cz, ux, uy, uz, nerr.
c
c Calls: aptmopv, aptqnor, aptqrts, aptrotv
c
c
c Glossary:
c
c ax,ay,az Input The x, y and z coordinates of point "a", at which to
c find the normal vector, extreme radii of curvature
c and the corresponding tangent vectors, for the
c quadric surface though point "a".
c
c bx,by,bz Output The x, y and z components at point "a" of the normal
c vector of the quadric surface through point "a".
c
c cx,cy,cz Output The x, y and z coordinates of point "c", the center of
c curvature corresponding to each value of rcurve.
c Size 4. c = a + rcurve * b / |b|.
c
c nerr Output Indicates an indeterminate point "a", if not zero:
c 1 if the normal vector "b" is zero at point "a".
c
c nrad Output Number of values of extreme radii of curvature rcurve
c and unit vectors "u" found.
c If zero, the normal vector "b" is zero at point "a",
c and the values are indeterminate.
c
c rcurve Output The nrad extreme values of the radius of curvature,
c in the direction of vector "b", of all possible
c curves in the quadric surface through point "a".
c Size 4.
c Actually, values for which the reciprocal of rcurve
c has extreme values, or is zero.
c
c qc-qzz Input Coefficients of the implicit equation of one quadric
c surface of the family of quadric surfaces represented
c by all possible values of qc. Any value of qc may be
c specified, and it will not be changed.
c
c tol Input Numerical tolerance limit.
c On 64-bit computers, recommend 1.e-11.
c On 32-bit computers, recommend 1.e-6.
c
c ux,uy,uz Output The x, y and z components of nrad unit vectors "u"
c tangent to the curves lying in the quadric surface
c through point "a", having the nrad extreme values
c rcurve of the radius of curvature. Size 4.
c
c History:
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
dimension cx(4) ! Coordinate x of a center of curvature.
dimension cy(4) ! Coordinate y of a center of curvature.
dimension cz(4) ! Coordinate z of a center of curvature.
dimension rcurve(4) ! An extreme radius of curvature.
dimension ux(4) ! Component x of a tangent unit vector.
dimension uy(4) ! Component y of a tangent unit vector.
dimension uz(4) ! Component z of a tangent unit vector.
c.... Local variables.
common /laptcris/ rotm(3,3) ! Rotation matrix, vector "d" to z axis.
common /laptcris/ theta(4) ! Angle of unit vector from x axis.
common /laptcris/ pi ! Mathematic constant, pi.
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptcris 001 finding extreme radii of curvature.' /
cbug & ' tol=',1pe13.5 /
cbug & ' ax,ay,az= ',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, ax, ay, az
cbug write ( 3, 9902) qc, qx, qy, qz, qxy, qyz, qzx, qxx, qyy, qzz
cbugc***DEBUG ends.
c.... Initialize.
pi = 3.14159265358979323 ! Mathematical constant, pi.
nerr = 0
nrad = 0
do 110 n = 1, 4
rcurve(n) = 0.0
cx(n) = -1.e99
cy(n) = -1.e99
cz(n) = -1.e99
ux(n) = -1.e99
uy(n) = -1.e99
uz(n) = -1.e99
110 continue
c.... Find the normal vector "b" at point "a" of the quadric surface.
call aptqnor (ax, ay, az, qc, qx, qy, qz, qxy, qyz, qzx,
& qxx, qyy, qzz, tol, bc, bx, by, bz, blen, nerra)
cbugcbugc***DEBUG begins.
cbugcbug 9961 format ('aptcris DEBUG. blen=',1pe22.14)
cbugcbug write ( 3, 9961) blen
cbugcbugc***DEBUG ends.
if (nerra .ne. 0) then
nerr = 1
go to 999
endif
c.... Find the rotation operator to rotate vector "b" to the z axis.
call aptrotv (bx, by, bz, 0., 0., 1., tol, rotm, nerra)
cbugcbugc***DEBUG begins.
cbugcbug 9990 format (' op11,op12,op13 ',1p3e20.12 /
cbugcbug & ' op21,op22,op23 ',1p3e20.12 /
cbugcbug & ' op31,op32,op33 ',1p3e20.12)
cbugcbug write ( 3, 9990) ((rotm(i,j), j = 1, 3), i = 1, 3)
cbugcbugc***DEBUG ends.
if (nerra .ne. 0) then
nerr = 1
go to 999
endif
c.... Rotate the quadric surface around point "a" to make vector "b"
c.... point in the z direction.
sc = qc
sx = qx
sy = qy
sz = qz
sxy = qxy
syz = qyz
szx = qzx
sxx = qxx
syy = qyy
szz = qzz
cbugcbugc***DEBUG begins.
cbugcbug 9960 format (' sc ',1pe20.12 )
cbugcbug 9980 format (' sx, sy, sz ',1p3e20.12 )
cbugcbug 9985 format (' sxy, syz, szx ',1p3e20.12 )
cbugcbug 9991 format (' sxx, syy, szz ',1p3e20.12 )
cbugcbug write ( 3, 9960) sc
cbugcbug write ( 3, 9980) sx, sy, sz
cbugcbug write ( 3, 9985) sxy, syz, szx
cbugcbug write ( 3, 9991) sxx, syy, szz
cbugcbugc***DEBUG ends.
call aptrois (rotm, 0, ax, ay, az, sc, sx, sy, sz,
& sxy, syz, szx, sxx, syy, szz, tol, nerra)
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9960) sc
cbugcbug write ( 3, 9980) sx, sy, sz
cbugcbug write ( 3, 9985) sxy, syz, szx
cbugcbug write ( 3, 9991) sxx, syy, szz
cbugcbugc***DEBUG ends.
c.... Find any infinite radii of curvature.
c.... sxx + sxy * tan(theta) + syy * tan(theta)**2 = 0.
c.... Solve equation for tan(theta).
call aptqrts (0, syy, sxy, sxx, qq, tol,
& nroots, root1, root2, itrun)
if (nroots .eq. -1) then
c.... Solve equation for cot(theta).
call aptqrts (0, sxx, sxy, syy, qq, tol,
& nroots, root1, root2, itrun)
if (nroots .eq. -1) then ! All angles same. No curvature.
nrad = 4
theta(1) = 0.0
theta(2) = 0.25 * pi
theta(3) = 0.50 * pi
theta(4) = 0.75 * pi
rcurve(1) = 1.e99
rcurve(2) = 1.e99
rcurve(3) = 1.e99
rcurve(4) = 1.e99
go to 120
elseif (nroots .eq. 1) then
nrad = 1
theta(1) = atan2 (1.0, root1)
rcurve(1) = 1.e99
elseif (nroots .eq. 2) then
nrad = 2
theta(1) = atan2 (1.0, root1)
theta(2) = atan2 (1.0, root2)
rcurve(1) = 1.e99
rcurve(2) = 1.e99
go to 120
endif
elseif (nroots .eq. 1) then
nrad = 1
theta(1) = atan (root1)
rcurve(1) = 1.e99
elseif (nroots .eq. 2) then
nrad = 2
theta(1) = atan (root1)
theta(2) = atan (root2)
rcurve(1) = 1.e99
rcurve(2) = 1.e99
endif
cbugcbugc***DEBUG begins.
cbugcbug 9981 format ('aptcris DEBUG. # of infinite radii =',i3)
cbugcbug if (nrad .ne. 0) then
cbugcbug write ( 3, 9981) nrad
cbugcbug endif
cbugcbugc***DEBUG ends.
c.... Find any extreme values of curvature.
c.... (syy - sxx) * sin(2*theta) + sxy * cos(2*theta) = 0
if ((syy .eq. sxx) .and. (sxy .eq. 0.0)) then ! All angles same.
nrad = 4
theta(1) = 0.0
theta(2) = 0.25 * pi
theta(3) = 0.50 * pi
theta(4) = 0.75 * pi
else ! Two angles 90 degrees apart.
sinth = sxy
costh = sxx - syy
angle = 0.5 * atan2 (sinth, costh)
nrad = nrad + 1
theta(nrad) = angle
angle = angle + 0.5 * pi
do 115 n = 1, nrad
dang = abs (angle - theta(n)) - pi
if (abs (dang) .le. tol) go to 120
115 continue
nrad = nrad + 1
theta(nrad) = angle
endif ! Tested coefficients.
c.... Find the tangent unit vectors and the radii and centers of curvature.
120 if (nrad .gt. 0) then ! Extreme radii of curvature found.
cbugcbugc***DEBUG begins.
cbugcbug 9971 format ('theta=',1p4e15.8)
cbugcbug write ( 3, 9971) (theta(n), n = 1, nrad)
cbugcbugc***DEBUG ends.
do 130 n = 1, nrad ! Loop over extreme radii.
costh = cos (theta(n))
sinth = sin (theta(n))
cbugcbugc***DEBUG begins.
cbugcbug 9983 format ('costh, sinth=',1p2e22.14)
cbugcbug write ( 3, 9983) costh, sinth
cbugcbugc***DEBUG ends.
if (abs (costh) .le. tol) costh = 0.0
if (abs (sinth) .le. tol) sinth = 0.0
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9983) costh, sinth
cbugcbugc***DEBUG ends.
ux(n) = costh
uy(n) = sinth
uz(n) = 0.0
den = sxx * costh**2 + sxy * costh*sinth + syy * sinth**2
denerr = 3.0 * tol * (abs (sxx) * costh**2 +
& abs (sxy * costh * sinth) +
& abs (syy) * sinth**2)
if (abs (den) .le. denerr) den = 0.0
if (den .eq. 0.0) then ! No curvature.
rcurve(n) = 1.e99
else ! Some curvature.
if (rcurve(n) .ne. 1.e99) then
rcurve(n) = -0.5 * blen / den
endif
endif
cbugcbugc***DEBUG begins.
cbugcbug 9972 format ('n, den, denerr, rcurve=',i3,2x,1p3e15.7)
cbugcbug write ( 3, 9972) n, den, denerr, rcurve(n)
cbugcbugc***DEBUG ends.
call aptvsum (0, 1.0, ax, ay, az, rcurve(n), 0., 0., 1.,
& 1, tol, cx(n), cy(n), cz(n), clen, nerr)
130 continue ! End of loop over extreme radii.
endif ! Tested nrad.
c.... Rotate any centers of curvature and tangent unit vectors back to the
c.... original direction.
if (nrad .gt. 0) then
call aptmopv (rotm, 1, 0., 0., 0., ux, uy, uz, nrad, tol, nerra)
call aptmopv (rotm, 1, ax, ay, az, cx, cy, cz, nrad, tol, nerra)
endif
c.... Remove any duplicate results.
if (nrad .gt. 1) then ! Remove any duplicate results.
nrads = nrad ! Number of original results.
nrad = 1 ! Initial number of unique results.
do 150 n = 2, nrads ! Loop over original results.
do 140 nn = 1, nrad ! Loop over unique results.
if ((ux(n) .eq. ux(nn)) .and.
& (uy(n) .eq. uy(nn)) .and.
& (uz(n) .eq. uz(nn)) .and.
& (cx(n) .eq. cx(nn)) .and.
& (cy(n) .eq. cy(nn)) .and.
& (cz(n) .eq. cz(nn)) .and.
& (rcurve(n) .eq. rcurve(nn))) go to 150 ! Not unique.
140 continue ! End of loop over unique results.
c.... This result is unique. Save.
nrad = nrad + 1
ux(nrad) = ux(n)
uy(nrad) = uy(n)
uz(nrad) = uz(n)
cx(nrad) = cx(n)
cy(nrad) = cy(n)
cz(nrad) = cz(n)
rcurve(nrad) = rcurve(n)
150 continue ! End of loop over original results.
endif ! Remove any duplicate results.
999 continue
cbugc***DEBUG begins.
cbug 9910 format (/ 'aptcris results: nrad=',i2,' nerr=',i2)
cbug 9911 format ( 'Normal (xyz)=',1p3e22.14 )
cbug 9912 format (/ 'Tangent (xyz)=',1p3e22.14 /
cbug & 'Center (xyz)=',1p3e22.14 /
cbug & 'Radius =',1pe22.14 )
cbug write ( 3, 9910) nrad, nerr
cbug write ( 3, 9911) bx, by, bz
cbug if (nrad .gt. 0) then
cbug do 991 n = 1, nrad
cbug write ( 3, 9912) ux(n), uy(n), uz(n),
cbug & cx(n), cy(n), cz(n), rcurve(n)
cbug 991 continue
cbug endif
cbugc***DEBUG ends.
return
c.... End of subroutine aptcris. (+1 line)
end
UCRL-WEB-209832