subroutine aptqexc (qc, qu, qv, quv, quu, qvv, tol,
& np, eu, ev, adir, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTQEXC
c
c call aptqexc (qc, qu, qv, quv, quu, qvv, tol,
c np, eu, ev, adir, nerr)
c
c Version: aptqexc Updated 1994 April 11 11:10.
c aptqexc Originated 1994 April 7 15:00
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To find the coordinates (eu, ev), type adir, and number np, of
c any points on a quadric curve in the uv plane with extreme
c values of u or v, or at an intersection of two straight lines.
c Up to four extreme points may exist.
c
c The implicit equation of the quadric curve is:
c
c f(u, v) = qc + qu * u + qv * v + quv * u * v +
c quu * u**2 + qvv * v**2 = 0 .
c
c Results may be unpredictable if the coefficients differ by
c many orders of magnitude.
c
c Method:
c
c The components of the normal vector to the curve, N, are the
c partial derivatives of f with respect to u and v:
c
c N = (df/du, df/dv) (partial derivatives)
c fu = df/du = qu + quv * v + 2 * quu * u (partial derivative)
c fv = df/dv = qv + quv * u + 2 * qvv * v (partial derivative)
c
c The complete first derivatives, u' and v', are:
c
c u' = du/dv = -fv / fu
c v' = dv/du = -fu / fv
c
c From zero to two extreme values of u occur where fv = 0, fu is
c not zero, and u' = 0. The extreme value is a minimum if u'' is
c positive, and a maximum if u'' is negative at the extreme point:
c
c u'' = d(u')/dv = -2.0 * qvv / fu, when u' = 0.
c
c From zero to two extreme values of v occur where fu = 0, fv is
c not zero, and v' = 0. The extreme value is a minimum if v'' is
c positive, and a maximum if v'' is negative at the extreme point:
c
c v'' = d(v')/du = -2.0 * quu / fv, when v' = 0.
c
c An intersection of two straight lines occurs where fu = fv = 0.
c
c General cases:
c
c Null curves (all coefficients zero): no extreme points.
c Skew straight lines (simple, coincident or parallel):
c no extreme points.
c Horizontal straight lines (simple or coincident):
c one constant value of v.
c Horizontal parallel lines: minimum and maximum values of v.
c Vertical straight lines (simple or coincident):
c one constant value of u.
c Vertical parallel lines: minimum and maximum values of u.
c Intersecting straight lines: one intersection point.
c Hyperbolas: zero, one or two extreme points. With two extreme
c points, the minimum of one curve will be larger than the
c maximum of the other curve.
c Parabolas: one extreme point if aligned with an axis, and two
c extreme points (one u, one v), if not.
c Circles or ellipses: four extreme points (two u, two v).
c
c Input: qc, qu, qv, quv, quu, qvv, tol.
c
c Output: np, eu, ev, adir, nerr.
c
c Calls: aptqrts
c
c Glossary:
c
c adir Output The type of extreme point:
c Minima: "umin" or "vmin".
c Maxima: "umax" or "vmax".
c Horizontal or vertical lines: "ucon" or "vcon".
c Intersection of two straight lines:
c "uint" or "vint".
c If none, return "unknown".
c Size = 4.
c
c eu, ev Output The u and v coordinates of an extreme point or an
c intersection of two straight lines.
c If none, return (-1.e99,-1.e99).
c Size = 4.
c
c nerr Output Indicates an input error, if not zero:
c 1 if all coefficients, except possible qc, are zero.
c
c np Output Number of extreme or intersection points found.
c
c qc In Constant term in the equation of the quadric curve.
c
c qu,qv In Coefficients of u and v, respectively, in the equation
c of the quadric curve.
c
c qxx In Coefficient of x**2 in the equation of the quadric
c curve.
c
c qxy In Coefficient of x * y in the equation of the quadric
c curve.
c
c qyy In Coefficient of y**2 in the equation of the quadric
c curve.
c
c tol Input Numerical tolerance limit.
c On Cray computers, recommend 1.e-11.
c On 32-bit computers, recommend 1.e-6.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
dimension eu(4) ! Coordinate u of extreme point.
dimension ev(4) ! Coordinate v of extreme point.
dimension adir(4) ! Normal direction of extreme point.
character*8 adir ! Normal direction of extreme point.
c.... Local variables.
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptqexc finding the extrema of quadric curve.' /
cbug & ' tol=',1pe13.5 )
cbug 9902 format (
cbug & ' qc= ',1pe22.14 /
cbug & ' qu, qv= ',1p2e22.14 /
cbug & ' quv= ',1pe22.14 /
cbug & ' quu, qvv= ',1p2e22.14)
cbug write ( 3, 9901) tol
cbug write ( 3, 9902) qc, qu, qv, quv, quu, qvv
cbugc***DEBUG ends.
c.... Initialize.
nerr = 0
np = 0
do 110 n = 1, 4
eu(n) = -1.e99
ev(n) = -1.e99
adir(n) = 'unknown'
110 continue
c.... Test for input errors.
if ((qu .eq. 0.0) .and.
& (qv .eq. 0.0) .and.
& (quv .eq. 0.0) .and.
& (quu .eq. 0.0) .and.
& (qvv .eq. 0.0)) then
nerr = 1
go to 210
endif
c.... Test for simple straight lines.
if ((quv .eq. 0.0) .and.
& (quu .eq. 0.0) .and.
& (qvv .eq. 0.0)) then
cbugc***DEBUG begins.
cbug 9903 format ('aptqexc DEBUG. Simple straight line.')
cbug write ( 3, 9903)
cbugc***DEBUG ends.
if (qv .eq. 0.0) then ! Straight line u = -qc / qu.
cbugc***DEBUG begins.
cbug 9904 format ('aptqexc DEBUG. Vertical straight line.')
cbug write ( 3, 9904)
cbugc***DEBUG ends.
np = np + 1
eu(np) = -qc / qu
ev(np) = 0.0
vnormu = qu
adir(np) = 'ucon'
go to 210
endif
if (qu .eq. 0.0) then ! Straight line v = -qc / qv.
cbugc***DEBUG begins.
cbug 9905 format ('aptqexc DEBUG. Horizontal straight line.')
cbug write ( 3, 9905)
cbugc***DEBUG ends.
np = np + 1
eu(np) = 0.0
ev(np) = -qc / qv
unormv = qv
adir(np) = 'vcon'
go to 210
endif
go to 210
endif
c.... Test for horizontal and vertical parallel lines.
if ((qv .eq. 0.0) .and.
& (quv .eq. 0.0) .and.
& (qvv .eq. 0.0)) then
cbugc***DEBUG begins.
cbug 9906 format ('aptqexc DEBUG. Vertical parallel lines.')
cbug write ( 3, 9906)
cbugc***DEBUG ends.
call aptqrts (0, quu, qu, qc, qq, tol, nrootu, uroot1, uroot2,
& itrun)
if (nrootu .ge. 1) then
np = np + 1
eu(np) = uroot1
ev(np) = 0.0
unormu = qu + 2.0 * quu * eu(np)
adir(np) = 'ucon'
endif
if (nrootu .eq. 2) then
np = np + 1
eu(np) = uroot2
ev(np) = 0.0
unormu = qu + 2.0 * quu * eu(np)
if (eu(np) .ge. eu(np-1)) then
adir(np) = 'umax'
adir(np-1) = 'umin'
else
adir(np) = 'umin'
adir(np-1) = 'umax'
endif
endif
go to 210
endif
if ((qu .eq. 0.0) .and.
& (quv .eq. 0.0) .and.
& (quu .eq. 0.0)) then
cbugc***DEBUG begins.
cbug 9907 format ('aptqexc DEBUG. Horizontal parallel lines.')
cbug write ( 3, 9907)
cbugc***DEBUG ends.
call aptqrts (0, qvv, qv, qc, qq, tol, nrootv, vroot1, vroot2,
& itrun)
if (nrootv .ge. 1) then
np = np + 1
eu(np) = 0.0
ev(np) = vroot1
vnormv = qv + 2.0 * qvv * ev(np)
adir(np) = 'vcon'
endif
if (nrootv .eq. 2) then
np = np + 1
eu(np) = 0.0
ev(np) = vroot2
vnormv = qv + 2.0 * qvv * ev(np)
if (ev(np) .ge. ev(np-1)) then
adir(np) = 'vmax'
adir(np-1) = 'vmin'
else
adir(np) = 'vmin'
adir(np-1) = 'vmax'
endif
endif
go to 210
endif
c.... Test for a quadric curve aligned with the axes.
if (quv .eq. 0.0) then ! Aligned quadric curve.
cbugc***DEBUG begins.
cbug 9908 format ('aptqexc DEBUG. Quadric aligned with axes.')
cbug write ( 3, 9908)
cbugc***DEBUG ends.
c.... Find points with extreme values of u, for aligned quadric.
if (qvv .ne. 0.0) then
v = -0.5 * qv / qvv
aa = quu
bb = qu
cc = qc + qv * v + qvv * v**2
call aptqrts (0, aa, bb, cc, qq, tol, nrootu, uroot1, uroot2,
& itrun)
if (nrootu .ge. 1) then
np = np + 1
eu(np) = uroot1
ev(np) = v
vnormu = qu + 2.0 * quu * eu(np)
if (qvv * vnormu .gt. 0.0) then
adir(np) = 'umax'
elseif (qvv * vnormu .lt. 0.0) then
adir(np) = 'umin'
else
adir(np) = 'uint'
endif
endif
if (nrootu .eq. 2) then
np = np + 1
eu(np) = uroot2
ev(np) = v
vnormu = qu + 2.0 * quu * eu(np)
if (qvv * vnormu .gt. 0.0) then
adir(np) = 'umax'
elseif (qvv * vnormu .lt. 0.0) then ! Exchange two extrema.
adir(np) = adir(np-1)
adir(np-1) = 'umin'
eu(np) = uroot1
eu(np-1) = uroot2
else
adir(np) = 'uint'
endif
endif
endif
c.... Find points with extreme values of v, for aligned quadric.
if (quu .ne. 0.0) then
u = -0.5 * qu / quu
aa = qvv
bb = qv
cc = qc + qu * u + quu * u**2
call aptqrts (0, aa, bb, cc, qq, tol, nrootv, vroot1, vroot2,
& itrun)
if (nrootv .ge. 1) then
np = np + 1
eu(np) = u
ev(np) = vroot1
vnormv = qv + 2.0 * qvv * ev(np)
if (quu * vnormv .gt. 0.0) then
adir(np) = 'vmax'
elseif (quu * vnormv .lt. 0.0) then
adir(np) = 'vmin'
else
adir(np) = 'vint'
endif
endif
if (nrootv .eq. 2) then
np = np + 1
eu(np) = u
ev(np) = vroot2
vnormv = qv + 2.0 * qvv * ev(np)
if (quu * vnormv .gt. 0.0) then
adir(np) = 'vmax'
elseif (quu * vnormv .lt. 0.0) then ! Exchange two extrema.
adir(np) = adir(np-1)
adir(np-1) = 'vmin'
ev(np) = vroot1
ev(np-1) = vroot2
else
adir(np) = 'vint'
endif
endif
endif
go to 210
endif ! Tested for aligned quadric curve.
c.... Test for an unaligned quadric curve.
if (quv .ne. 0.0) then ! Unaligned quadric curve.
cbugc***DEBUG begins.
cbug 9909 format ('aptqexc DEBUG. Unaligned quadric.')
cbug write ( 3, 9909)
cbugc***DEBUG ends.
c.... Find points with extreme values of u, for unaligned quadric.
if (qvv .ne. 0.0) then
aa = qvv * (4.0 * qvv * quu - quv**2)
bb = 2.0 * qvv * (2.0 * qv * quu - qu * quv)
cc = qv**2 * quu - qv * qu * quv + qc * quv**2
call aptqrts (0, aa, bb, cc, qq, tol, nrootv, vroot1, vroot2,
& itrun)
if (nrootv .ge. 1) then
np = np + 1
eu(np) = -(qv + 2.0 * qvv * vroot1) / quv
ev(np) = vroot1
vnormu = qu + quv * ev(np) + 2.0 * quu * eu(np)
if (qvv * vnormu .gt. 0.0) then
adir(np) = 'umax'
elseif (qvv * vnormu .lt. 0.0) then
adir(np) = 'umin'
else
adir(np) = 'uint'
endif
endif
if (nrootv .eq. 2) then
np = np + 1
eu(np) = -(qv + 2.0 * qvv * vroot2) / quv
ev(np) = vroot2
vnormu = qu + quv * ev(np) + 2.0 * quu * eu(np)
if (qvv * vnormu .gt. 0.0) then
adir(np) = 'umax'
elseif (qvv * vnormu .lt. 0.0) then ! Exchange two extrema.
adir(np) = adir(np-1)
adir(np-1) = 'umin'
eux = eu(np)
eu(np) = eu(np-1)
eu(np-1) = eux
ev(np) = vroot1
ev(np-1) = vroot2
else
adir(np) = 'uint'
endif
endif
endif
c.... Find points with extreme values of v, for unaligned quadric.
if (quu .ne. 0.0) then
aa = quu * (4.0 * quu * qvv - quv**2)
bb = 2.0 * quu * (2.0 * qu * qvv - qv * quv)
cc = qu**2 * qvv - qu * qv * quv + qc * quv**2
call aptqrts (0, aa, bb, cc, qq, tol, nrootu, uroot1, uroot2,
& itrun)
if (nrootu .ge. 1) then
np = np + 1
eu(np) = uroot1
ev(np) = -(qu + 2.0 * quu * uroot1) / quv
vnormv = qv + quv * eu(np) + 2.0 * qvv * ev(np)
if (quu * vnormv .gt. 0.0) then
adir(np) = 'vmax'
elseif (quu * vnormv .lt. 0.0) then
adir(np) = 'vmin'
else
adir(np) = 'vint'
endif
endif
if (nrootu .eq. 2) then
np = np + 1
eu(np) = uroot2
ev(np) = -(qu + 2.0 * quu * uroot2) / quv
vnormv = qv + quv * eu(np) + 2.0 * qvv * ev(np)
if (quu * vnormv .gt. 0.0) then
adir(np) = 'vmax'
elseif (quu * vnormv .lt. 0.0) then ! Exchange two extrema.
adir(np) = adir(np-1)
adir(np-1) = 'vmin'
eu(np) = uroot1
eu(np-1) = uroot2
evx = ev(np)
ev(np) = ev(np-1)
ev(np-1) = evx
else
adir(np) = 'vint'
endif
endif
endif
go to 210
endif ! Tested for unaligned quadric curve.
210 continue
cbugc***DEBUG begins.
cbug 9910 format (/ 'aptqexc results: np=',i2,' nerr=',i2)
cbug 9911 format (i3,' eu, ev= ',1p2e22.14,2x,a8 )
cbug write ( 3, 9910) np, nerr
cbug do 991 n = 1, np
cbug write ( 3, 9911) n, eu(n), ev(n), adir(n)
cbug 991 continue
cbugc***DEBUG ends.
return
c.... End of subroutine aptqexc. (+1 line.)
end
UCRL-WEB-209832