subroutine aptrkrl (pr, ur, sr, dr, dintmn, dintmx, np, tol, & nint, dint, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTRKRL c c call aptrkrl (pr, ur, sr, dr, dintmn, dintmx, np, tol, c & nint, dint, nerr) c c Version: aptrkrl Updated 1993 December 6 15:20. c aptrkrl Originated 1990 January 19 16:20. 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, any exit c intersection of the linear track through point p = (pr) with c unit direction vector u = (ur, ut, up), with the spherical c surface with fixed radius sr, for which (1) the distance dint c from point "p" to the intersection is between the limits dintmn c and dintmx, and (2) the radial component ur of the direction c vector "u" at the intersection has the same sign as dr. c Flag nerr indicates any input error. c c Input: pr, ur, sr, dr, dintmn, dintmx, np, tol. c c Output: nint, dint, nerr. c c Calls: aptqrtv c c c Glossary: c c dint Output Distance from point "p" to intersection at radius sr c along track (not radial distance). Size np. c c dintmn Input Minimum allowed value of distance to intersection. c Size np. c c dintmx Input Maximum allowed value of distance to intersection. c Size np. c c dr Input Sign of exit direction through surface at sr. Size np. c c nerr Output If not 0, indicates an input error. c 1 if np is not positive. c c nint Output O if no exit intersection was found. c 1 if an exit intersection was found at the surface c at radius sr, with a distance to intersection dint c between dintmn and dintmx. Size np. c c np Input Size of arrays pr, ur, sr, dr, dintmn, dintmx, nint. c c pr Input Spherical radial coordinate of initial point on c track. Size np. c c sr Input Spherical radial coordinate of surface. Size np. c c tol Input Truncation error limit. c Used to test for intersection being nearly c tangent, and for accuracy of intersection. c Must not be zero. c On Cray computers, recommend 1.e-11. c c ur Input Initial spherical radial component of unit direction c vector along track. Size np. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccend. c.... Dimensioned arguments. c---- Distance to intersection. dimension dint (1) c---- Minimum distance to intersection. dimension dintmn (1) c---- Maximum distance to intersection. dimension dintmx (1) c---- Sign of exit direction at sr. dimension dr (1) c---- Number of intersections (0 or 1). dimension nint (1) c---- Coordinate r of beginning of track. dimension pr (1) c---- Coordinate r of spherical surface. dimension sr (1) c---- Component r of direction vector. dimension ur (1) c.... Local variables. c---- Coefficient of dint**2 in quadratic. common /laptrkrl/ aa (64) c---- Coefficient of dint in quadratic. common /laptrkrl/ bb (64) c---- A very big number. common /laptrkrl/ big c---- Constant term in quadratic. common /laptrkrl/ cc (64) c---- Distance to intersection 1. common /laptrkrl/ dint1 (64) c---- Distance to intersection 2. common /laptrkrl/ dint2 (64) c---- A very small number. common /laptrkrl/ fuz c---- 1 = tangent, 2 = bad qq. common /laptrkrl/ itrun (64) c---- Index in local array. common /laptrkrl/ n c---- First index of subset of data. common /laptrkrl/ n1 c---- Last index of subset of data. common /laptrkrl/ n2 c---- Index in external array. common /laptrkrl/ nn c---- 1 if root1 an exit. common /laptrkrl/ nroot1 c---- 1 if root2 an exit. common /laptrkrl/ nroot2 c---- Number of real roots of quadratic. common /laptrkrl/ nroots (64) c---- Size of current subset of data. common /laptrkrl/ ns c---- Factor bb**2 - 4.0 * aa * cc. common /laptrkrl/ qq (64) c---- New ur at intersection 1. common /laptrkrl/ ur1 c---- New ur at intersection 2. common /laptrkrl/ ur2 cbugc***DEBUG begins. cbug 9901 format (/ 'aptrkrl finding intersections. np=',i3, cbug & ' tol=',1pe13.5) cbug 9902 format (i3,' pr= ',1pe22.14 / cbug & ' ur= ',1pe22.14 / cbug & ' sr,dr= ',1p2e22.14 / cbug & ' dintmn=',1pe22.14,' dintmx=',1pe22.14) cbug write ( 3, 9901) np, tol cbug write ( 3, 9902) (n, pr(n), ur(n), cbug & sr(n), dr(n), dintmn(n), dintmx(n), n = 1, np) cbugc***DEBUG ends. c.... Initialize. c---- A very big number. big = 1.e+99 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.... Do data sets in groups of 64. n1 = 1 n2 = min (np, 64) 110 ns = n2 - n1 + 1 c.... Find coefficients of aa * dint**2 + bb * dint + cc = 0. c---- Loop over local arrays. do 120 n = 1, ns nn = n + n1 - 1 aa(n) = 1.0 bb(n) = 2.0 * pr(nn) * ur(nn) cc(n) = pr(nn)**2 - sr(nn)**2 c.... The following calculation of qq is more accurate than c.... qq(n) = bb(n)**2 - 4.0 * aa(n) * cc(n). qq(n) = 4.0 * (sr(nn)**2 - pr(nn)**2 * (1.0 - ur(nn)**2)) c---- End of loop over local arrays. 120 continue c.... Find any intersections. call aptqrtv (1, aa, bb, cc, qq, ns, tol, & nroots, dint1, dint2, itrun, nerr) c.... See if either solution has dint between dintmn and dintmx, and c.... has an intersection ur with same sign as dr. c---- Loop over local arrays. do 130 n = 1, ns nn = n + n1 - 1 ur1 = (pr(nn) * ur(nn) + aa(n) * dint1(n)) / (sr(nn) + fuz) ur2 = (pr(nn) * ur(nn) + aa(n) * dint2(n)) / (sr(nn) + fuz) if ((nroots(n) .ge. 1) .and. & (dint1(n) .ge. dintmn(nn)) .and. & (dint1(n) .le. dintmx(nn)) .and. & ((ur1 * dr(nn)) .gt. 0.0) ) then nroot1 = 1 else nroot1 = 0 endif if (nroot1 .eq. 1) then nint(nn) = 1 dint(nn) = dint1(n) nroots(n) = 1 else nint(nn) = 0 dint(nn) = -big endif if ((nroots(n) .eq. 2) .and. & (dint2(n) .ge. dintmn(nn)) .and. & (dint2(n) .le. dintmx(nn)) .and. & ((ur2 * dr(nn)) .gt. 0.0) ) then nroot2 = 1 else nroot2 = 0 endif if (nroot2 .eq. 1) then nint(nn) = 1 dint(nn) = dint2(n) endif c---- End of loop over local arrays. 130 continue cbugc***DEBUG begins. cbug 9903 format (/ 'aptrkrl results:' / cbug & (i3,' nint=',i2,' dint=',1pe22.14)) cbug write ( 3, 9903) (n, nint(n), dint(n), n = 1, np) cbugc***DEBUG ends. c.... See if all sets are done. c---- Do another data set. if (n2 .lt. np) then n1 = n2 + 1 n2 = min (np, n1 + 63) go to 110 endif 210 return c.... End of subroutine aptrkrl. (+1 line.) end UCRL-WEB-209832