subroutine aptdist (nsys, iunit, au, av, aw, bu, bv, bw,
& np, tol, cu, cv, cw, dab, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTDIST
c
c call aptdist (nsys, iunit, au, av, aw, bu, bv, bw,
c & np, tol, cu, cv, cw, dab, nerr)
c
c Version: aptdist Updated 1990 March 14 16:00.
c aptdist Originated 1989 November 27 16:50.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To find, for each of the np sets of input data, the distance
c dab, and initial unit direction vector c = (cu, cv, cw),
c from the point a = (au, av, aw) to the point b = (bu, bv, bw).
c Option nsys specifies the coordinate system:
c 0 for Cartesian, 1 for cylindrical, 2 for spherical.
c Option iunit indicates the units for angles: 0 for degrees,
c 1 for radians.
c Any component of vector "c" less than the estimated error in
c its calculation, based on tol, will be truncated to zero.
c Vector "c" will be zero, if point "a" is coincident with point
c "b", based on tol.
c Flag nerr indicates any input error.
c
c History: 1990 March 14. Changed tol to 0.0 in call to unit vector
c subroutine. Allows small magnitudes.
c
c Input: nsys, iunit, au, av, aw, bu, bv, bw, np, tol.
c
c Output: cu, cv, cw, dab, nerr.
c
c Calls: aptcsys, aptcsyv, aptvdis, aptvuna
c
c
c Glossary:
c
c au,av,aw Input The u, v, w coordinates of point "a". Size np.
c
c bu,bv,bw Input The u, v, w coordinates of point "b". Size np.
c
c cu,cv,cw Output The u, v, w components of vector "c". Size np.
c May be truncated to zero, if less than the estimated
c numerical error in their calculation based on tol.
c
c dab Output The distance from point "a" to point "b". Size np.
c
c iunit Input Indicates unit to be used for angles:
c 0 if angles are in degrees.
c 1 if angles are in radians.
c
c nerr Output Indicates an input error, if not 0.
c 1 if np is not positive.
c 2 if nsys is not 1, 2, or 3.
c 3 if iunit is not 0 or 1.
c
c np Input Size of arrays.
c
c nsys Input Indicates coordinate system type:
c 0 for cartesian coordinates. u = x, v = y, w = z.
c 1 for cylindrical coordinates. u = radius from z
c axis, v = angle in xy plane, counterclockwise from
c x axis, w = z.
c 2 for spherical coordinates. u = radius from origin,
c v = angle in xy plane, counterclockwise from
c x axis, w = angle from z axis.
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.
c---- Coordinate u of input point "a".
dimension au (1)
c---- Coordinate v of input point "a".
dimension av (1)
c---- Coordinate w of input point "a".
dimension aw (1)
c---- Coordinate u of input point "b".
dimension bu (1)
c---- Coordinate v of input point "b".
dimension bv (1)
c---- Coordinate w of input point "b".
dimension bw (1)
c---- Component u of unit vector "c".
dimension cu (1)
c---- Component v of unit vector "c".
dimension cv (1)
c---- Component w of unit vector "c".
dimension cw (1)
c---- The distance from "a" to "b".
dimension dab (1)
c.... Local variables.
c---- Temporary coordinate u of point "a".
common /laptdist/ aus (64)
c---- Temporary coordinate v of point "a".
common /laptdist/ avs (64)
c---- Temporary coordinate w of point "a".
common /laptdist/ aws (64)
c---- Temporary coordinate u of point "b".
common /laptdist/ bus (64)
c---- Temporary coordinate v of point "b".
common /laptdist/ bvs (64)
c---- Temporary coordinate w of point "b".
common /laptdist/ bws (64)
c---- Index in arrays.
common /laptdist/ n
c---- First index of subset of data.
common /laptdist/ n1
c---- Last index of subset of data.
common /laptdist/ n2
c---- Size of current subset of data.
common /laptdist/ ns
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptdist finding distance and direction:' /
cbug & ' nsys=',i2,' iunit=',i2 /
cbug & (i3,' au,av,aw=',1p3e22.14 /
cbug & ' bu,bv,bw=',1p3e22.14))
cbug write ( 3, 9901) nsys, iunit, (n, au(n), av(n), aw(n),
cbug & bu(n), bv(n), bw(n), n = 1, np)
cbugc***DEBUG ends.
c.... Initialize.
nerr = 0
c.... Test for input errors.
if (np .le. 0) then
nerr = 1
go to 210
endif
if ((nsys .lt. 0) .or. (nsys .gt. 2)) then
nerr = 2
go to 210
endif
if ((iunit .lt. 0) .or. (iunit .gt. 1)) then
nerr = 3
go to 210
endif
c.... Set up the indices of the first subset of data.
n1 = 1
n2 = min (np, 64)
c.... Loop over subsets of data.
110 ns = n2 - n1 + 1
c.... Store temporary values of input.
c---- Loop over subset of data.
do 120 n = 1, ns
aus(n) = au(n+n1-1)
avs(n) = av(n+n1-1)
aws(n) = aw(n+n1-1)
bus(n) = bu(n+n1-1)
bvs(n) = bv(n+n1-1)
bws(n) = bw(n+n1-1)
c---- End of loop over subset of data.
120 continue
c.... See if conversion to Cartesian coordinates is needed.
c---- Convert points "a" and "b".
if (nsys .ne. 0) then
call aptcsys (nsys, 0, iunit, aus, avs, aws, ns, tol, nerr)
call aptcsys (nsys, 0, iunit, bus, bvs, bws, ns, tol, nerr)
endif
c.... Find the distance dab and vector "c" from point "a" to point "b".
call aptvdis (aus, avs, aws, bus, bvs, bws, ns, tol,
& cu(n1), cv(n1), cw(n1), dab, nerr)
c.... Make vector "c" into a unit vector.
call aptvuna (cu(n1), cv(n1), cw(n1), ns, 0., dab, nerr)
c.... See if conversion back to initial coordinates is needed.
c---- Convert vector "c" at point "a".
if (nsys .ne. 0) then
call aptcsyv (0, nsys, iunit, aus, avs, aws,
& cu(n1), cv(n1), cw(n1), ns, tol, nerr)
endif
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 9902 format (/ 'aptdist distance and direction:' /
cbug & (i3,' dab= ',1pe22.14 /
cbug & ' cu,cv,cw=',1p3e22.14))
cbug write ( 3, 9902) (n, dab(n),
cbug & cu(n), cv(n), cw(n), n = 1, np)
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptdist. (+1 line.)
end
UCRL-WEB-209832