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