subroutine aptmvcy (kth, psz, psr, cths, sths, usz, usr, ust,
& dpmove, np, tol,
& pz, pr, cth, sth, uz, ur, ut, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTMVCY
c
c call aptmvcy (kth, psz, psr, cths, sths, usz, usr, ust,
c & dpmove, np, tol,
c & pz, pr, cth, sth, uz, ur, ut, nerr)
c
c Version: aptmvcy Updated 1990 November 28 14:50.
c aptmvcy Originated 1989 December 4 11:00.
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 new
c point p = (pz, pr, cth, sth) and unit direction vector
c u = (uz, ur, ut), resulting from moving from the point
c ps = (psz, psr, cths, sths) in the direction of the unit
c vector us = (usz, usr, ust) for a distance dpmove, in
c cylindrical coordinates.
c If kth = 0, all cths = 1.0, all sths = 0.0, and cth and sth
c will not be calculated, and none need be dimensioned.
c Any component of point "p" or vector "u" less than the
c estimated error in its calculation, based on tol, will be
c truncated to zero. If tol = 0, no truncation tests are done.
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 1990 March 30. Fixed array index error which affected problems
c with np .gt. 64.
c
c Input: kth, psz, psr, cths, sths, usz, usr, ust, dpmove, np, tol.
c
c Output: pz, pr, cth, sth, uz, ur, ut, nerr.
c
c Calls: aptvuna
c
c Glossary:
c Size np.
c
c cth Output The cosine of the final value of theta (angle in the
c xy plane counterclockwise from x axis).
c May be truncated to zero, if less than the estimated
c numerical error in their calculation based on tol.
c Size np, if kth = 1. Otherwise, not calculated.
c
c cths Input The cosine of the initial value of theta. See cth.
c Size np, if kth = 1. Otherwise, scalar 1.0.
c
c dpmove Input The distance from point "ps" to point "p". Size np.
c (Assuming vector "us" is a unit vector.)
c
c kth Input Indicates size of arrays cth, cths, sth, sths:
c 0 if array size is 1, with cths = 1.0, sths = 0.0,
c and cth and sth are not to be calculated.
c 1 if array size is np, input values of cths and sths
c will be used, and cth and sth will be calculated.
c
c nerr Output Indicates an input error, if not 0.
c 1 if np is not positive.
c 2 if kth is not 0 or 1.
c
c np Input Size of arrays.
c
c pr, pz Output The r and z components of final point "p". Size np.
c May be truncated to zero, if less than the estimated
c numerical error in their calculation based on tol.
c
c psr, psz Input The r and z coordinates of initial point "ps".
c Size np.
c
c sth Output The sine of the final value of theta. See cth.
c May be truncated to zero, if less than the estimated
c numerical error in their calculation based on tol.
c Size np, if kth = 1. Otherwise, not calculated.
c
c sths Input The sine of the initial value of theta. See cth.
c Size np, if kth = 1. Otherwise, scalar 0.0.
c
c tol Input Numerical tolerance limit. If zero, no tests done.
c On Cray computers, recommend 1.e-5 to 1.e-11.
c
c ur Output The r component of final unit direction vector "u".
c May be truncated to zero, if less than the estimated
c numerical error in their calculation based on tol.
c Size np.
c
c usr Input The r component of initial unit direction vector "us".
c Size np.
c
c ust Input The theta component of initial unit direction vector
c "us". See cth. Size np.
c
c usz Input The z component of initial unit direction vector "us".
c Size np.
c
c ut Output The theta component of final unit direction vector "u".
c May be truncated to zero, if less than the estimated
c numerical error in their calculation based on tol.
c See cth. Size np.
c
c uz Output The z component of final unit direction vector "u".
c Size np.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c---- Cosine of final theta.
dimension cth (1)
c---- Cosine of initial theta.
dimension cths (1)
c---- The distance from "ps" to "p".
dimension dpmove (1)
c---- Coordinate r of point "p".
dimension pr (1)
c---- Coordinate r of point "ps".
dimension psr (1)
c---- Coordinate z of point "ps".
dimension psz (1)
c---- Coordinate z of point "p".
dimension pz (1)
c---- Sine of final theta.
dimension sth (1)
c---- Sine of initial theta.
dimension sths (1)
c---- Component r of unit vector "u".
dimension ur (1)
c---- Component r of unit vector "us".
dimension usr (1)
c---- Theta component of unit vector "us".
dimension ust (1)
c---- Component z of unit vector "us".
dimension usz (1)
c---- Theta component of unit vector "u".
dimension ut (1)
c---- Component z of unit vector "u".
dimension uz (1)
c.... Local variables.
c---- Cosine of angle theta.
common /laptmvcy/ costh
c---- A very small number.
common /laptmvcy/ fuz
c---- Index in arrays.
common /laptmvcy/ n
c---- First index of subset of data.
common /laptmvcy/ n1
c---- Last index of subset of data.
common /laptmvcy/ n2
c---- Size of current subset of data.
common /laptmvcy/ ns
c---- Coordinate x at point "ps".
common /laptmvcy/ px
c---- Estimated error in px.
common /laptmvcy/ pxerr
c---- Coordinate y at point "ps".
common /laptmvcy/ py
c---- Estimated error in py.
common /laptmvcy/ pyerr
c---- Estimated error in pz.
common /laptmvcy/ pzerr
c---- Sine of angle theta.
common /laptmvcy/ sinth
c---- Estimated error in ur.
common /laptmvcy/ urerr
c---- Estimated error in ut.
common /laptmvcy/ uterr
c---- Component x of unit vector "us".
common /laptmvcy/ ux
c---- Estimated error in ux.
common /laptmvcy/ uxerr
c---- Component y of unit vector "us".
common /laptmvcy/ uy
c---- Estimated error in uy.
common /laptmvcy/ uyerr
c---- Magnitude of new direction vector.
common /laptmvcy/ vlen (64)
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptmvcy initial point, direction and distance',
cbug & ' (kth=',i2,'):' /
cbug & (i3,' pz,pr= ',1p2e22.14 /
cbug & ' uz,ur,ut=',1p3e22.14 /
cbug & ' dpmove= ',1pe22.14))
cbug 9902 format (i3,' cth,sth= ',1p2e22.14)
cbug write ( 3, 9901) kth, (n, psz(n), psr(n),
cbug & usz(n), usr(n), ust(n), dpmove(n), n = 1, np)
cbug if (kth .eq. 1) then
cbug write ( 3, 9902) (n, cths(n), sths(n), n = 1, np)
cbug endif
cbugc***DEBUG ends.
c.... Initialize.
c---- A very small number.
fuz = 1.e-99
nerr = 0
c.... Test for input errors.
if (np .le. 0) then
nerr = 1
go to 210
endif
if ((kth .lt. 0) .or. (kth .gt. 1)) then
nerr = 2
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.... Find new point "p" and new unit direction vector "u".
c---- No truncation.
if (tol .le. 0.0) then
c---- No calculation of cth, sth.
if (kth .eq. 0) then
c---- Loop over subset of data.
do 120 n = n1, n2
pz(n) = psz(n) + dpmove(n) * usz(n)
px = psr(n) + dpmove(n) * usr(n)
py = dpmove(n) * ust(n)
pr(n) = sqrt (px**2 + py**2)
uz(n) = usz(n)
ur(n) = ( usr(n) * (px + fuz) + ust(n) * py) /
& (pr(n) + fuz)
ut(n) = (-usr(n) * py + ust(n) * (px + fuz)) /
& (pr(n) + fuz)
c---- End of loop over subset of data.
120 continue
c---- Calculate cth, sth.
else
c---- Loop over subset of data.
do 130 n = n1, n2
pz(n) = psz(n) + dpmove(n) * usz(n)
ux = usr(n) * cths(n) - ust(n) * sths(n)
uy = usr(n) * sths(n) + ust(n) * cths(n)
px = psr(n) * cths(n) + dpmove(n) * ux
py = psr(n) * sths(n) + dpmove(n) * uy
pr(n) = sqrt (px**2 + py**2)
cth(n) = (px + fuz) / (pr(n) + fuz)
sth(n) = py / (pr(n) + fuz)
uz(n) = usz(n)
ur(n) = ux * cth(n) + uy * sth(n)
ut(n) = -ux * sth(n) + uy * cth(n)
c---- End of loop over subset of data.
130 continue
c---- Tested kth.
endif
c---- Truncate small values to zero.
else
c---- No calculation of cth, sth.
if (kth .eq. 0) then
c---- Loop over subset of data.
do 140 n = n1, n2
pz(n) = psz(n) + dpmove(n) * usz(n)
pzerr = tol * (abs (psz(n)) + abs (dpmove(n) * usz(n)))
if (abs (pz(n)) .lt. pzerr) then
pz(n) = 0.0
endif
px = psr(n) + dpmove(n) * usr(n)
pxerr = tol * (abs (psr(n)) + abs (dpmove(n) * usr(n)))
if (abs (px) .lt. pxerr) then
px = 0.0
endif
py = dpmove(n) * ust(n)
pr(n) = sqrt (px**2 + py**2)
costh = (px + fuz) / (pr(n) + fuz)
if (abs (costh) .lt. tol) then
costh = 0.0
endif
sinth = py / (pr(n) + fuz)
if (abs (sinth) .lt. tol) then
sinth = 0.0
endif
uz(n) = usz(n)
ur(n) = usr(n) * costh + ust(n) * sinth
urerr = tol * (abs (usr(n) * costh) +
& abs (ust(n) * sinth))
if (abs (ur(n)) .lt. urerr) then
ur(n) = 0.0
endif
ut(n) = -usr(n) * sinth + ust(n) * costh
uterr = tol * (abs (usr(n) * sinth) +
& abs (ust(n) * costh))
if (abs (ut(n)) .lt. uterr) then
ut(n) = 0.0
endif
c---- End of loop over subset of data.
140 continue
c---- Calculate cth, sth.
else
c---- Loop over subset of data.
do 150 n = n1, n2
pz(n) = psz(n) + dpmove(n) * usz(n)
pzerr = tol * (abs (psz(n)) + abs (dpmove(n) * usz(n)))
if (abs (pz(n)) .lt. pzerr) then
pz(n) = 0.0
endif
costh = cths(n)
if (abs (costh) .lt. tol) then
costh = 0.0
endif
sinth = sths(n)
if (abs (sinth) .lt. tol) then
sinth = 0.0
endif
ux = usr(n) * costh - ust(n) * sinth
uxerr = tol * (abs (usr(n) * costh) +
& abs (ust(n) * sinth))
if (abs (ux) .lt. uxerr) then
ux = 0.0
endif
uy = usr(n) * sinth + ust(n) * costh
uyerr = tol * (abs (usr(n) * sinth) +
& abs (ust(n) * costh))
if (abs (uy) .lt. uyerr) then
uy = 0.0
endif
px = psr(n) * costh + dpmove(n) * ux
pxerr = tol * (abs (psr(n) * costh) +
& abs (dpmove(n) * ux))
if (abs (px) .lt. pxerr) then
px = 0.0
endif
py = psr(n) * sinth + dpmove(n) * uy
pyerr = tol * (abs (psr(n) * sinth) +
& abs (dpmove(n) * uy))
if (abs (py) .lt. pyerr) then
py = 0.0
endif
pr(n) = sqrt (px**2 + py**2)
uz(n) = usz(n)
costh = (px + fuz) / (pr(n) + fuz)
if (abs (costh) .lt. tol) then
costh = 0.0
endif
sinth = py / (pr(n) + fuz)
if (abs (sinth) .lt. tol) then
sinth = 0.0
endif
cth(n) = costh
sth(n) = sinth
ur(n) = ux * costh + uy * sinth
urerr = tol * (abs (ux * costh) + abs (uy * sinth))
if (abs (ur(n)) .lt. urerr) then
ur(n) = 0.0
endif
ut(n) = -ux * sinth + uy * costh
uterr = tol * (abs (ux * sinth) + abs (uy * costh))
if (abs (ut(n)) .lt. uterr) then
ut(n) = 0.0
endif
c---- End of loop over subset of data.
150 continue
c---- Tested kth.
endif
c---- Tested tol.
endif
c.... Renormalize unit direction vector.
call aptvuna (ur(n1), ut(n1), uz(n1), ns, 0., vlen, nerr)
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 9903 format (/ 'aptmvcy results: new point and direction:' /
cbug & (i3,' pz,pr= ',1p2e22.14 /
cbug & ' uz,ur,ut=',1p3e22.14))
cbug write ( 3, 9903) (n, pz(n), pr(n),
cbug & uz(n), ur(n), ut(n), n = 1, np)
cbug if (kth .eq. 1) then
cbug write ( 3, 9902) (n, cths(n), sths(n), n = 1, np)
cbug endif
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptmvcy. (+1 line.)
end
UCRL-WEB-209832