subroutine aptintq (px, py, pz,
& ac, ax, ay, az, axy, ayz, azx, axx, ayy, azz,
& bc, bx, by, bz, bxy, byz, bzx, bxx, byy, bzz,
& tol, aquint, tax, tay, taz, tbx, tby, tbz,
& dx, dy, dz, dist, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTINTQ
c
c call aptintq (px, py, pz,
c & ac, ax, ay, az, axy, ayz, azx, axx, ayy, azz,
c & bc, bx, by, bz, bxy, byz, bzx, bxx, byy, bzz,
c & tol, aquint, tax, tay, taz, tbx, tby, tbz,
c & dx, dy, dz, dist, nerr)
c
c Version: aptintq Updated 1998 May 14 17:20.
c aptintq Originated 1997 December 18 11:50.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To search for, starting with point p = (px, py, pz), the
c proximal points ta = (tax, tay, taz) and tb = (tbx, tby, tbz),
c and the direction d = (dx, dy, dz) and distance dist between
c ta and tb, or an intersection point (ta = tb), of the two
c quadric surfaces with the implicit equations:
c
c A = ac + ax * x + ay * y + az * z
c + axy * x * y + ayz * y * z + azx * z * x
c + axx * x**2 + ayy * y**2 + azz * z**2 = 0
c B = bc + bx * x + by * y + bz * z
c + bxy * x * y + byz * y * z + bzx * z * x
c + bxx * x**2 + byy * y**2 + bzz * z**2 = 0
c
c The two quadric surfaces have the normal vectors:
c
c an = (anx, any, anz) = grad (A)
c anx = ax + 2 * axx * x + axy * y + azx * z
c any = ay + axy * x + 2 * ayy * y + ayz * z
c anz = az + azx * x + ayz * y + 2 * azz * z
c
c bn = (bnx, bny, bnz) = grad (B)
c bnx = bx + 2 * bxx * x + bxy * y + bzx * z
c bny = by + bxy * x + 2 * byy * y + byz * z
c bnz = bz + bzx * x + byz * y + 2 * bzz * z
c
c The method is iterative, starting with the initial guess
c (x, y, z) = (px, py, pz). Each step of the iteration, the
c estimated intersection or proximal point is successively
c replaced by its proximal point on the other surface (which
c works well for surfaces that are highly curved or meet nearly
c orthogonally).
c
c Convergence is determined by the precision limit tol, which
c should be in the range from 1.e-5 to 1.e-15 for 64-bit
c floating point arithmetic.
c
c Input: px, py, pz,
c ac, ax, ay, az, axy, ayz, azx, axx, ayy, azz,
c bc, bx, by, bz, bxy, byz, bzx, bxx, byy, bzz, tol.
c
c Output: aquint, tax, tay, taz, tbx, tby, tbz, nerr.
c
c Calls: aptplpl, aptvdis, aptwhis, aptvaxb
c
c
c Glossary:
c
c ac, ... Input Coefficients of the implicit equation for the first
c quadric surface:
c ac + ax * x + ay * y + az * z +
c axy * x * y + ayz * y * z + azx * z * x +
c axx * x**2 + ayy * y**2 + azz * z**2 = 0
c
c aquint Output An 8-character ASCII word indicating the nature of
c points "ta" and "tb" on the two quadric surfaces:
c
c "intcurve" Two coincident points on an intersection
c curve of the two surfaces.
c "proximal" Two nearest points on non-intersecting
c surfaces.
c "tangent" Two coincident points at a tangency of
c the two surfaces.
c "unknown" If nerr is not zero, method failed.
c If nerr = 0, ambiguous results (probably
c a tangency)
c
c bc, ... Input Coefficients of the implicit equation for the second
c quadric surface:
c bc + bx * x + by * y + bz * z +
c bxy * x * y + byz * y * z + bzx * z * x +
c bxx * x**2 + byy * y**2 + bzz * z**2 = 0
c
c dist Output The distance from point "ta" to point "tb".
c Zero if the surfaces intersect at point "ta".
c
c dx,dy,dz Output The x, y and z components of the vector from point
c "ta" to point "tb". All zero if the surfaces
c intersect at point "ta".
c
c nerr Output Indicates an input error, if not zero:
c 1 if an equation of a quadric surface is bad.
c 2 if the method fails to converge in 1000 iterations,
c which may indicate too small a value of tol.
c 3 if the the method diverges.
c
c px,py,pz Input The x, y, z coordinates of point "p", the initial guess
c for the intersection point.
c
c tax,y,z Output Coordinates x, y and z at an intersection point of
c surfaces A and B (dist = 0) or proximal point
c on surface A (dist not zero), if nerr = 0.
c Unless the two surfaces are tangent at a single
c point, an intersection point is not unique, but
c only one point on a continuous intersection curve,
c usually one near point "p".
c
c tbx,y,z Output Coordinates x, y and z at an intersection point of
c surfaces A and B (dist = 0) or proximal point
c on surface B (dist not zero), if nerr = 0.
c Unless the two surfaces are tangent at a single
c point, an intersection point is not unique, but
c only one point on a continuous intersection curve,
c usually one near point "p".
c
c tol Input Numerical tolerance limit. With 64-bit floating point
c numbers, 1.e-5 to 1.e-11 is recommended.
c Used both as an absolute and relative convergence
c criterion. Changes in point coordinates less then
c 10 * tol are treated as insignificant.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Arguments.
character*8 aquint ! "intcurve", "proximal" or "tangent".
c.... Local variables.
common /ltintq/ dprox(4) ! Distance to a proximal point.
common /ltintq/ xx(4) ! Coordinate x of a proximal point.
common /ltintq/ yy(4) ! Coordinate x of a proximal point.
common /ltintq/ zz(4) ! Coordinate x of a proximal point.
c=======================================================================********
c.... Initialize.
nerr = 0
tolx = tol
if (tolx .eq. 0.0) tolx = 1.e-11
tax = -1.e99
tay = -1.e99
taz = -1.e99
tbx = -1.e99
tby = -1.e99
tbz = -1.e99
dx = -1.e99
dy = -1.e99
dz = -1.e99
dist = -1.e99
aquint = "unknown"
c=======================================================================********
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptintq finding intersection point of two',
cbug & ' quadric surfaces.',
cbug & ' tol=',1pe13.5)
cbug 9902 format ('Initial guess:' /
cbug & ' px, py, pz =',1p3e22.14 )
cbug 9903 format ('Quadric ac =',1pe22.14 /
cbug & ' ax,ay,az =',1p3e22.14 /
cbug & ' axy,ayz,azx=',1p3e22.14 /
cbug & ' axx,ayy,azz=',1p3e22.14 )
cbug 9904 format ('Quadric bc =',1pe22.14 /
cbug & ' bx,by,bz =',1p3e22.14 /
cbug & ' bxy,byz,bzx=',1p3e22.14 /
cbug & ' bxx,byy,bzz=',1p3e22.14 )
cbug write ( 3, 9901) tolx
cbug write ( 3, 9902) px, py, pz
cbug write ( 3, 9903) ac, ax, ay, az, axy, ayz, azx, axx, ayy, azz
cbug write ( 3, 9904) bc, bx, by, bz, bxy, byz, bzx, bxx, byy, bzz
cbugc***DEBUG ends.
c.... Initialize for the iteration.
tax = px
tay = py
taz = pz
tbx = px
tby = py
tbz = pz
sidea = -1.e99
sideb = -1.e99
kprox = 0 ! 1 if two proximal points found.
kintr = 0 ! 1 if an intersection found.
nitex = 100 ! Iteration to begin extrapolation method.
nitstop = 1000 ! Iteration to stop.
cbugcbugc***DEBUG begins.
cbugcbug nit = 0
cbugcbug write ( 3, 9911) nit, tax, tay, taz
cbugcbug write ( 3, 9912) nit, tbx, tby, tbz
cbugcbugc***DEBUG ends.
c=======================================================================********
c.... Iterate to find the solution.
do 380 nit = 1, 1000 ! Iterate to find solution.
c.... Replace point "ta" with the point on surface A proximal to point "tb".
taxs = tax
tays = tay
tazs = taz
sidebs = sideb
call aptwhis (tbx, tby, tbz,
& ac, ax, ay, az, axy, ayz, azx, axx, ayy, azz,
& tolx, sideb, anx, any, anz,
& nprox, xx, yy, zz, dprox, nerra)
if ((nerra .ne. 0) .or.
& (nprox .le. 0)) then
cbugcbugc***DEBUG begins.
cbugcbug 9913 format (/ 'DEBUG: nerra=',i2,' nprox=',i2)
cbugcbug write ( 3, 9913) nerra, nprox
cbugcbugc***DEBUG ends.
nerr = 1
go to 999
endif
tax = xx(1)
tay = yy(1)
taz = zz(1)
c.... Find the vector "an" normal to surface A at new point "ta".
call aptwhis (tax, tay, taz,
& ac, ax, ay, az, axy, ayz, azx, axx, ayy, azz,
& tolx, sideaa, anx, any, anz,
& nprox, xx, yy, zz, dprox, nerra)
cbugcbugc***DEBUG begins.
cbugcbug 9911 format (/ i4,1x,'tax,y,z =',1p3e22.14 )
cbugcbug write ( 3, 9911) nit, tax, tay, taz
cbugcbugc***DEBUG ends.
c.... Replace point "tb" with the point on surface B proximal to point "ta".
tbxs = tbx
tbys = tby
tbzs = tbz
sideas = sidea
call aptwhis (tax, tay, taz,
& bc, bx, by, bz, bxy, byz, bzx, bxx, byy, bzz,
& tolx, sidea, bnx, bny, bnz,
& nprox, xx, yy, zz, dprox, nerra)
if ((nerra .ne. 0) .or.
& (nprox .le. 0)) then
cbugcbugc***DEBUG begins.
cbugcbug write ( 3, 9913) nerra, nprox
cbugcbugc***DEBUG ends.
nerr = 1
go to 999
endif
tbx = xx(1)
tby = yy(1)
tbz = zz(1)
c.... Find the vector "bn" normal to surface B at new point "tb".
call aptwhis (tbx, tby, tbz,
& bc, bx, by, bz, bxy, byz, bzx, bxx, byy, bzz,
& tolx, sidebb, bnx, bny, bnz,
& nprox, xx, yy, zz, dprox, nerra)
cbugcbugc***DEBUG begins.
cbugcbug 9912 format (i4,1x,'tbx,y,z =',1p3e22.14 )
cbugcbug write ( 3, 9912) nit, tbx, tby, tbz
cbugcbugc***DEBUG ends.
cbugcbugc***DEBUG begins.
cbugcbug 9801 format (i4,1x,'dtax,y,z=',1p3e22.14 )
cbugcbug 9802 format (i4,1x,'dtbx,y,z=',1p3e22.14 )
cbugcbug 9803 format (i4,1x,'erax,y,z=',1p3e22.14 )
cbugcbug 9804 format (i4,1x,'erbx,y,z=',1p3e22.14 )
cbugcbug dtax = tax - taxs
cbugcbug dtay = tay - tays
cbugcbug dtaz = taz - tazs
cbugcbug dtbx = tbx - tbxs
cbugcbug dtby = tby - tbys
cbugcbug dtbz = tbz - tbzs
cbugcbug write ( 3, 9801) nit, dtax, dtay, dtaz
cbugcbug write ( 3, 9802) nit, dtbx, dtby, dtbz
cbugcbug erax = tolx * (10.0 + abs (taxs))
cbugcbug eray = tolx * (10.0 + abs (tays))
cbugcbug eraz = tolx * (10.0 + abs (tazs))
cbugcbug erbx = tolx * (10.0 + abs (tbxs))
cbugcbug erby = tolx * (10.0 + abs (tbys))
cbugcbug erbz = tolx * (10.0 + abs (tbzs))
cbugcbug write ( 3, 9803) nit, erax, eray, eraz
cbugcbug write ( 3, 9804) nit, erbx, erby, erbz
cbugcbugc***DEBUG ends.
c.... See if the test points have converged.
if (kprox .eq. 2) go to 410
if ((abs (tax - taxs) .le. tolx * (10.0 + abs (taxs))) .and.
& (abs (tay - tays) .le. tolx * (10.0 + abs (tays))) .and.
& (abs (taz - tazs) .le. tolx * (10.0 + abs (tazs))) .and.
& (abs (tbx - tbxs) .le. tolx * (10.0 + abs (tbxs))) .and.
& (abs (tby - tbys) .le. tolx * (10.0 + abs (tbys))) .and.
& (abs (tbz - tbzs) .le. tolx * (10.0 + abs (tbzs)))) then
cbugcbugc***DEBUG begins.
cbugcbug 9915 format ('Converged in',i5,' iterations (ta and tb).')
cbugcbug write ( 3, 9915) nit
cbugcbugc***DEBUG ends.
kprox = kprox + 1
endif
c.... Try method of extrapolating to an intersection point.
if ((nit .eq. nitex) .and. (kprox .eq. 0)) then ! Try extrapolation.
cold call aptlnln (taxs, tays, tazs, tax, tay, taz,
cold & tbxs, tbys, tbzs, tbx, tby, tbz, 1, tolx,
cold & dpmin, fracab, fraccd,
cold & taxe, taye, taze, tbxe, tbye, tbze,
cold & itrun, ipar, nerr)
call aptplpl (tax, tay, taz, anx, any, anz,
& tbx, tby, tbz, bnx, bny, bnz,
& 1, tolx,
& taxe, taye, taze, tbxe, tbye, tbze,
& ux, uy, uz, ipar, dpmin, itrun, nerra)
if ((ipar .eq. 0) .and.
& (itrun .eq. 0) .and.
& (nerra .eq. 0)) then
cbugcbugc***DEBUG begins.
cbugcbug 9711 format (i4,1x,'tbx,y,ze=',1p3e22.14 )
cbugcbug 9712 format (i4,1x,'tbx,y,ze=',1p3e22.14 )
cbugcbug write ( 3, 9711) nit, taxe, taye, taze
cbugcbug write ( 3, 9712) nit, tbxe, tbye, tbze
cbugcbugc***DEBUG ends.
tax = taxe
tay = taye
taz = taze
tbx = tbxe
tby = tbye
tbz = tbze
sideas = -1.e99
sidebs = -1.e99
nitex = nitex + 2
endif
endif ! Tried extrapolation method.
c.... See if the test points are at an intersection point.
cbugcbugc***DEBUG begins.
cbugcbug 9921 format (i4,1x,'sidea,b =',1p2e22.14 )
cbugcbug write ( 3, 9921) nit, sidea, sideb
cbugcbugc***DEBUG ends.
if ((sidea .eq. 0.0) .and.
& (sideb .eq. 0.0)) then
cbugcbugc***DEBUG begins.
cbugcbug 9914 format ('Converged in',i5,' iterations (sidea and sideb).')
cbugcbug write ( 3, 9914) nit
cbugcbugc***DEBUG ends.
kintr = 1
if (kprox .gt. 0) go to 410
endif
c.... See if the method is diverging.
if ((abs (sidea) .gt. abs (sideas)) .and.
& (abs (sideb) .gt. abs (sidebs))) then
cbugcbugc***DEBUG begins.
cbugcbug 9931 format ('aptintq diverging. sidea and sideb increasing.')
cbugcbug write ( 3, 9931)
cbugcbugc***DEBUG ends.
nerr = 3
go to 999
endif
380 continue ! End of each iteration.
nit = nit - 1
cbugcbugc***DEBUG begins.
cbugcbug 9961 format (/ 'aptintq failed to converge in',i5,' iterations.')
cbugcbug write ( 3, 9961) nit
cbugcbugc***DEBUG ends.
nerr = 2 ! Failed to converge.
go to 999
410 continue
c.... Find the distance between points "ta" and "tb".
call aptvdis (tax, tay, taz, tbx, tby, tbz, 1, tolx,
& dx, dy, dz, dist, nerra)
c.... Find the cross product of the two normal vectors at "ta" and "tb".
call aptvaxb (anx, any, anz, bnx, bny, bnz, 1, tolx,
& dnx, dny, dnz, vlength, nerra)
c.... Find the nature of points "ta" and "tb".
if ((sidea .eq. 0.0) .and. (sideb .eq. 0.0)) then
if (dist .eq. 0.0) then
if (vlength .eq. 0.0) then
aquint = "tangent"
else
aquint = "intcurve"
endif
endif
elseif ((sidea .ne. 0.0) .and. (sideb .ne. 0.0)) then
if ((dist .ne. 0.0) .and. (vlength .eq. 0.0)) then
aquint = "proximal"
endif
endif
c=======================================================================********
999 continue
cbugc***DEBUG begins.
cbug 9991 format (/ 'aptintq results. nerr=',i2,' nit=',i4,
cbug & ' kprox=',i2,' kintr=',i2,' aquint=',a8)
cbug 9992 format ('tax,tay,taz =',1p3e22.14 )
cbug 9993 format ('tbx,tby,tbz =',1p3e22.14 )
cbug 9994 format ('sidea, sideb =',1p2e22.14)
cbug 9995 format ('vlength =',1pe22.14 )
cbug 9996 format ('dx, dy, dz =',1p3e22.14 /
cbug & 'dist =',1pe22.14 )
cbug 9997 format ('anx, any, anz=',1p3e22.14 )
cbug 9998 format ('bnx, bny, bnz=',1p3e22.14 )
cbug write ( 3, 9991) nerr, nit, kprox, kintr, aquint
cbug write ( 3, 9992) tax, tay, taz
cbug write ( 3, 9993) tbx, tby, tbz
cbug write ( 3, 9994) sidea, sideb
cbug write ( 3, 9995) vlength
cbug write ( 3, 9996) dx, dy, dz, dist
cbug write ( 3, 9997) anx, any, anz
cbug write ( 3, 9998) bnx, bny, bnz
cbugc***DEBUG ends.
return
c.... End of subroutine aptintq. (+1 line)
end
UCRL-WEB-209832