subroutine aptdirt (ax, ay, az, npa, bx, by, bz, npb, tol,
& ux, uy, uz, dab, ndir, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTDIRT
c
c call aptdirt (ax, ay, az, npa, bx, by, bz, npb, tol,
c & ux, uy, uz, dab, ndir, nerr)
c
c Version: aptdirt Updated 2003 June 10 12:00.
c aptdirt Originated 2003 May 21 17:30.
c
c Authors: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To find all unique direction vectors from each point
c a(na) = (ax(na), ay(na), az(na)), na = 1, npa, to each point
c b(nb) = (bx(nb), by(nb), bz(nb)), nb = 1, npb, within precision
c tol, and return them as unit direction vectors
c v(nab) = (vx(nab), vy(nab), vz(nab)), with minimum distances
c dab(nab), nab = 1, ndir, sorted in order of increasing distance.
c Flag nerr indicates any input error.
c
c Input: ax, ay, az, npa, bx, by, bz, npb, tol
c
c Output: ux, uy, uz, dab, ndir, nerr
c
c Calls: aptdist, aptvaxb
c
c Glossary:
c
c ax,ay,az Input The x, y and z coordinates of each point
c a(na) = (ax(na), ay(na), az(na)), na = 1, npa.
c Size npa or greater.
c If array b is sufficiently regular, all unique
c directions in b may be found by using only one
c or two points in b as array a.
c
c bx,by,bz Input The x, y and z coordinates of each point
c b(nb) = (bx(nb), by(nb), bz(nb)), nb = 1, npb.
c Size npb or greater,
c
c dab(nab) Output The minimum distance between pairs of points for each
c unique direction vector u(nab), nab = 1, ndir.
c Size npa * npb is sufficient, but see ndir.
c
c ndir Output The number of unique unit direction vectors u(nab),
c n = 1, ndir, connecting pairs of points.
c The maximum possible value of ndir is npa * npb.
c
c For a rectangular array of nx points in the
c x direction and ny <= nx points in the y direction,
c the maximum possible value of ndir is ndir(max) =
c 2 * nx * ny - 2 * nx - 4 * ny + 8.
c The minimum possible value is about F * ndir(max),
c where F is given by the product of all the terms
c (1 - 1/n**2), where n ranges over all of the prime
c numbers from 2 to the largest possible common prime
c factor of na and nb.
c For very large nx and ny, F = 0.607927565940.
c
c nerr Output Indicates an input error, if not zero:
c 1 if npa is less than 1.
c 2 if npb is less than 1.
c
c npa Input The number of points a(na) = (ax(na), ay(na), az(na)).
c
c npb Input The number of points b(nb) = (bx(nb), by(nb), bz(nb)).
c
c tol Input Numerical precision limit. Used to determine if
c points are coincident, or if components of vectors
c are zero.
c For 64-bit real numbers, recommend 1.e-5 to 1.e-11.
c
c ux,uy,uz Output The x, y and z components of each unique unit direction
c vector u(nab) = (ux(nab), uy(nab), uz(nab),
c nab = 1, ndir, connecting any point "a" to any point
c "b".
c The components ux, uy and ux are the cosines of the
c angles between the direction vectors and the positive
c x, y and z axes, respectively. The sign of each
c vector is arbitrary chosen to make the last
c non-zero component of ux, uy, uz positive.
c Size npa * npb is sufficient, but see ndir.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Declarations for arguments.
dimension ax(1) ! Coordinate x of point a.
dimension ay(1) ! Coordinate y of point a.
dimension az(1) ! Coordinate z of point a.
dimension bx(1) ! Coordinate x of point b.
dimension by(1) ! Coordinate y of point b.
dimension bz(1) ! Coordinate z of point b.
dimension ux(1) ! Component x of unit direction vector u.
dimension uy(1) ! Component y of unit direction vector u.
dimension uz(1) ! Component z of unit direction vector u.
dimension dab(1) ! Minimum distance between points.
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptdirt finding unique directions between points.',
cbug & ' npa =',i8,', npb =',i8,'.')
cbug 9902 format (/ ('n,ax,ay,az=',i6,1p3e20.12))
cbug 9912 format (/ ('n,bx,by,bz=',i6,1p3e20.12))
cbug write ( 3, 9901) npa, npb
cbug write ( 3, 9902) (n, ax(n), ay(n), az(n), n = 1, npa)
cbug write ( 3, 9912) (n, bx(n), by(n), bz(n), n = 1, npb)
cbugc***DEBUG ends.
c.... Initialize.
ndir = 0
c.... Test for input errors.
nerr = 0
if (npa .lt. 1) then
nerr = 1
go to 210
endif
if (npb .le. 1) then
nerr = 2
go to 210
endif
c.... Find the unique unit direction vectors between points.
do 140 na = 1, npa ! Loop over the first point.
do 130 nb = 1, npb ! Loop over the second point.
c.... Find the unit direction vector and distance between the two points.
call aptdist (0, 1, ax(na), ay(na), az(na),
& bx(nb), by(nb), bz(nb), 1, tol,
& abx, aby, abz, dabt, nerr)
if (dabt .eq. 0.0) go to 130 ! Points are coincident.
if (ndir .eq. 0) go to 120
do 110 nd = 1, ndir ! Loop over unique directions.
c.... See if directions are the same (cross product zero).
call aptvaxb (abx, aby, abz, ux(nd), uy(nd), uz(nd), 1, tol,
& cx, cy, cz, vlen, nerr)
if (vlen .eq. 0.0) then ! Directions are the same.
c.... See if the new distance is shorter.
if (dabt .lt. dab(nd)) then ! New distance is shorter.
dab(nd) = dabt
endif ! New distance is shorter.
go to 130
endif ! Directions are the same.
110 continue ! End of loop over unique directions.
c.... Found a new direction.
120 ndir = ndir + 1
ux(ndir) = abx
uy(ndir) = aby
uz(ndir) = abz
dab(ndir) = dabt
c.... Reverse direction if last non-zero component of ux, uy, uz is < 0.
if ((uz(ndir) .lt. 0.0) .or.
& ((uz(ndir) .eq. 0.0) .and.
& (uy(ndir) .lt. 0.0)) .or.
& ((uz(ndir) .eq. 0.0) .and.
& (uy(ndir) .eq. 0.0) .and.
& (ux(ndir) .lt. 0.0))) then
ux(ndir) = -ux(ndir)
uy(ndir) = -uy(ndir)
uz(ndir) = -uz(ndir)
endif
130 continue ! End of loop over the first point.
140 continue ! End of loop over the second point.
c.... Reorder in order of increasing minimum distance between points,
c.... or if that distance is the same, in order of increasing angle from the
c.... positive z direction, or if that angle is zero, in order of increasing
c.... angle from the positive y direction.
do 160 ns1 = 1, ndir - 1 ! Loop over directions.
do 150 ns2 = ns1 + 1, ndir ! Loop over rest of directions.
if ( (dab(ns2) .lt. dab(ns1)) .or.
& ((dab(ns2) .eq. dab(ns1)) .and.
& ((uz(ns2) .gt. uz(ns1) ) .or.
& ((uz(ns2) .eq. uz(ns1) ) .and.
& (uy(ns2) .gt. uy(ns1) ))))) then
uxsave = ux(ns1)
uysave = uy(ns1)
uzsave = uz(ns1)
dabsave = dab(ns1)
ux(ns1) = ux(ns2)
uy(ns1) = uy(ns2)
uz(ns1) = uz(ns2)
dab(ns1) = dab(ns2)
ux(ns2) = uxsave
uy(ns2) = uysave
uz(ns2) = uzsave
dab(ns2) = dabsave
endif
150 continue ! End of loop over directions.
160 continue ! End of loop over rest of directions.
210 continue
cbugc***DEBUG begins.
cbug 9903 format (/ 'aptdirt results: ndir=',i6,
cbug & ', nerr=',i2,'.' / )
cbug 9904 format ('n,ux,y,z,dab=',i6,1p4e15.7)
cbug write ( 3, 9903) ndir, nerr
cbug if (ndir .gt. 0) then
cbug write ( 3, 9904) (n, ux(n), uy(n), uz(n), dab(n), n = 1, ndir)
cbug endif
cbugc***DEBUG ends.
return
c.... End of subroutine aptdirt. (+1 line.)
end
UCRL-WEB-209832