subroutine aptmovs (rhos, cths, sths, cphs, sphs,
& usrh, usth, usph, dpmove, np, tol,
& rho, cth, sth, cph, sph, urh, uth, uph, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTMOVS
c
c call aptmovs (rhos, cths, sths, cphs, sphs,
c & usrh, usth, usph, dpmove, np, tol,
c & rho, cth, sth, cph, sph, urh, uth, uph, nerr)
c
c Version: aptmovs Updated 1990 November 28 14:50.
c aptmovs Originated 1989 December 4 17: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 = (rho, cth, sth, cph, sph) and unit direction vector
c u = (urh, uth, uph), resulting from moving from the point
c ps = (rhos, cths, sths, cphs, sphs) in the direction of the
c unit vector us = (usrh, usth, usph) for a distance dpmove,
c in spherical coordinates.
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: rhos, cths, sths, cphs, sphs, usrh, usth, usph, dpmove, np, tol.
c
c Output: rho, cth, sth, cph, sph, urh, uth, uph, nerr.
c
c Calls: aptvuna
c
c Glossary:
c
c cph Output The cosine of the final value of phi (angle from the
c z axis).
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 cphs Input The cosine of the initial value of phi. See uph.
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.
c
c cths Input The cosine of the initial value of theta. See cth.
c Size np.
c
c dpmove Input The distance from point "ps" to point "p". Size np.
c (Assuming vector "us" is a unit vector.)
c
c nerr Output Indicates an input error, if not 0.
c 1 if np is not positive.
c
c np Input Size of arrays.
c
c rho Output The spherical radial component of final point "p".
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 rhos Input The spherical radial component of initial point "ps".
c Size np.
c
c sphs Input The sine of the initial value of phi. See cph.
c Size np.
c
c sph Output The sine of the final value of phi. See cph.
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 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.
c
c sths Input The sine of the initial value of theta. See cth.
c Size np.
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 uph Output The phi 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 cph. Size np.
c
c urh Output The rho 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 usph Input The phi component of initial unit direction vector
c "us". See cph. Size np.
c
c usrh Input The rho component of initial unit direction vector
c "us". Size np.
c
c usth Input The theta component of initial unit direction vector
c "us". See cth. Size np.
c
c uth 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
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c---- Cosine of final phi.
dimension cph (1)
c---- Cosine of initial phi.
dimension cphs (1)
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---- Distance from origin to point "p".
dimension rho (1)
c---- Distance from origin to point "ps".
dimension rhos (1)
c---- Sine of final phi.
dimension sph (1)
c---- Sine of initial phi.
dimension sphs (1)
c---- Sine of final theta.
dimension sth (1)
c---- Sine of initial theta.
dimension sths (1)
c---- Phi component of unit vector "u".
dimension uph (1)
c---- Radial component of unit vector "u".
dimension urh (1)
c---- Phi component of unit vector "us".
dimension usph (1)
c---- Radial component of unit vector "us".
dimension usrh (1)
c---- Theta component of unit vector "us".
dimension usth (1)
c---- Theta component of unit vector "u".
dimension uth (1)
c.... Local variables.
c---- Cosine of angle phi.
common /laptmovs/ cosph
c---- Cosine of angle theta.
common /laptmovs/ costh
c---- A very small number.
common /laptmovs/ fuz
c---- Index in arrays.
common /laptmovs/ n
c---- First index of subset of data.
common /laptmovs/ n1
c---- Last index of subset of data.
common /laptmovs/ n2
c---- Size of current subset of data.
common /laptmovs/ ns
c---- Coordinate x at point "ps".
common /laptmovs/ px
c---- Estimated error in px.
common /laptmovs/ pxerr
c---- Coordinate y at point "ps".
common /laptmovs/ py
c---- Estimated error in py.
common /laptmovs/ pyerr
c---- Coordinate z at point "ps".
common /laptmovs/ pz
c---- Estimated error in pz.
common /laptmovs/ pzerr
c---- Distance from z axis.
common /laptmovs/ rcyl
c---- Sine of angle phi.
common /laptmovs/ sinph
c---- Sine of angle theta.
common /laptmovs/ sinth
c---- Estimated error in uph.
common /laptmovs/ upherr
c---- Estimated error in urh.
common /laptmovs/ urherr
c---- Estimated error in uth.
common /laptmovs/ utherr
c---- Component x of unit vector "us".
common /laptmovs/ ux
c---- Estimated error in ux.
common /laptmovs/ uxerr
c---- Component y of unit vector "us".
common /laptmovs/ uy
c---- Estimated error in uy.
common /laptmovs/ uyerr
c---- Component z of unit vector "us".
common /laptmovs/ uz
c---- Estimated error in uz.
common /laptmovs/ uzerr
c---- Magnitude of final direction vector.
common /laptmovs/ vlen (64)
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptmovs initial point, distance and direction:' /
cbug & (i3,' rho= ',1pe22.14,' dpmove=',1pe22.14 /
cbug & ' cth,sth= ',1p2e22.14 /
cbug & ' cph,sph= ',1p2e22.14 /
cbug & ' urh,uth,uph=',1p3e22.14))
cbug write ( 3, 9901) (n, dpmove(n), rhos(n), cths(n), sths(n),
cbug & cphs(n), sphs(n), usrh(n), usth(n), usph(n), n = 1, np)
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
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---- Loop over subset of data.
do 120 n = n1, n2
ux = usrh(n) * cths(n) * sphs(n) - usth(n) * sths(n) +
& usph(n) * cths(n) * cphs(n)
uy = usrh(n) * sths(n) * sphs(n) + usth(n) * cths(n) +
& usph(n) * sths(n) * cphs(n)
uz = usrh(n) * cphs(n) - usph(n) * sphs(n)
cbugc***DEBUG begins.
cbug px = rhos(n) * cths(n) * sphs(n)
cbug py = rhos(n) * sths(n) * sphs(n)
cbug pz = rhos(n) * cphs(n)
cbug 9801 format ("initial x,y,z=",1p3e22.14 /
cbug & " ux,uy,uz=",1p3e22.14)
cbug write ( 3, 9801) px, py, pz, ux, uy, uz
cbugc***DEBUG ends.
px = rhos(n) * cths(n) * sphs(n) + dpmove(n) * ux
py = rhos(n) * sths(n) * sphs(n) + dpmove(n) * uy
pz = rhos(n) * cphs(n) + dpmove(n) * uz
cbugc***DEBUG begins.
cbug 9802 format ("final x,y,z=",1p3e22.14)
cbug write ( 3, 9802) px, py, pz
cbugc***DEBUG ends.
rcyl = sqrt (px**2 + py**2)
cth(n) = (px + fuz) / (rcyl + fuz)
sth(n) = py / (rcyl + fuz)
rho(n) = sqrt (px**2 + py**2 + pz**2)
cph(n) = pz / (rho(n) + fuz)
sph(n) = (rcyl + fuz) / (rho(n) + fuz)
urh(n) = ux * cth(n) * sph(n) + uy * sth(n) * sph(n) +
& uz * cph(n)
uth(n) = -ux * sth(n) + uy * cth(n)
uph(n) = ux * cth(n) * cph(n) + uy * sth(n) * cph(n) -
& uz * sph(n)
c---- End of loop over subset of data.
120 continue
c---- Truncate small values to zero.
else
c---- Loop over subset of data.
do 130 n = n1, n2
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
cosph = cphs(n)
if (abs (cosph) .lt. tol) then
cosph = 0.0
endif
sinph = sphs(n)
if (abs (sinph) .lt. tol) then
sinph = 0.0
endif
ux = usrh(n) * costh * sinph - usth(n) * sinth +
& usph(n) * costh * cosph
uxerr = tol * (abs (usrh(n) * costh * sinph) +
& abs (usth(n) * sinth) +
& abs (usph(n) * costh * cosph))
if (abs (ux) .lt. uxerr) then
ux = 0.0
endif
uy = usrh(n) * sinth * sinph + usth(n) * costh +
& usph(n) * sinth * cosph
uyerr = tol * (abs (usrh(n) * sinth * sinph) +
& abs (usth(n) * costh) +
& abs (usph(n) * sinth * cosph))
if (abs (uy) .lt. uyerr) then
uy = 0.0
endif
uz = usrh(n) * cosph - usph(n) * sinph
uzerr = tol * (abs (usrh(n) * cosph) +
& abs (usph(n) * sinph))
if (abs (uz) .lt. uzerr) then
uz = 0.0
endif
cbugc***DEBUG begins.
cbug px = rhos(n) * cths(n) * sphs(n)
cbug py = rhos(n) * sths(n) * sphs(n)
cbug pz = rhos(n) * cphs(n)
cbug write ( 3, 9801) px, py, pz, ux, uy, uz
cbugc***DEBUG ends.
px = rhos(n) * costh * sinph + dpmove(n) * ux
pxerr = tol * (abs (rhos(n) * costh * sinph) +
& abs (dpmove(n) * ux))
if (abs (px) .lt. pxerr) then
px = 0.0
endif
py = rhos(n) * sinth * sinph + dpmove(n) * uy
pyerr = tol * (abs (rhos(n) * sinth * sinph) +
& abs (dpmove(n) * uy))
if (abs (py) .lt. pyerr) then
py = 0.0
endif
pz = rhos(n) * cosph + dpmove(n) * uz
pzerr = tol * (abs (rhos(n) * cosph) + abs (dpmove(n) * uz))
if (abs (pz) .lt. pzerr) then
pz = 0.0
endif
cbugc***DEBUG begins.
cbug write ( 3, 9802) px, py, pz
cbugc***DEBUG ends.
rcyl = sqrt (px**2 + py**2)
costh = (px + fuz) / (rcyl + fuz)
if (abs (costh) .lt. tol) then
costh = 0.0
endif
sinth = py / (rcyl + fuz)
if (abs (sinth) .lt. tol) then
sinth = 0.0
endif
cth(n) = costh
sth(n) = sinth
rho(n) = sqrt (px**2 + py**2 + pz**2)
cosph = pz / (rho(n) + fuz)
if (abs (cosph) .lt. tol) then
cosph = 0.0
endif
sinph = (rcyl + fuz) / (rho(n) + fuz)
if (abs (sinph) .lt. tol) then
sinph = 0.0
endif
cph(n) = cosph
sph(n) = sinph
urh(n) = ux * costh * sinph + uy * sinth * sinph +
& uz * cosph
urherr = tol * (abs (ux * costh * sinph) +
& abs (uy * sinth * sinph) +
& abs (uz * cosph))
if (abs (urh(n)) .lt. urherr) then
urh(n) = 0.0
endif
uth(n) = -ux * sinth + uy * costh
utherr = tol * (abs (ux * sinth) + abs (uy * costh))
if (abs (uth(n)) .lt. utherr) then
uth(n) = 0.0
endif
uph(n) = ux * costh * cosph + uy * sinth * cosph -
& uz * sinph
upherr = tol * (abs (ux * costh * cosph) +
& abs (uy * sinth * cosph) +
& abs (uz * sinph))
if (abs (uph(n)) .lt. upherr) then
uph(n) = 0.0
endif
c---- End of loop over subset of data.
130 continue
c---- Tested tol.
endif
c.... Renormalize the unit direction vector.
call aptvuna (urh(n1), uth(n1), uph(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 9902 format (/ 'aptmovs results: final point and direction:' /
cbug & (i3,' rho= ', 1pe22.14 /
cbug & ' cth,sth= ',1p2e22.14 /
cbug & ' cph,sph= ',1p2e22.14 /
cbug & ' urh,uth,uph=',1p3e22.14))
cbug write ( 3, 9902) (n, rho(n), cth(n), sth(n),
cbug & cph(n), sph(n), urh(n), uth(n), uph(n), n = 1, np)
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptmovs. (+1 line.)
end
UCRL-WEB-209832