subroutine aptrkcl (pr, ur, ut, uz, sr, dr, dintmn, dintmx, np, & tol, nint, dint, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTRKCL c c call aptrkcl (pr, ur, ut, uz, sr, dr, dintmn, dintmx, np, tol, c & nint, dint, nerr) c c Version: aptrkcl Updated 1993 December 6 15:20. c aptrkcl 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, uz), with the cylindrical 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, ut, uz, 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, ut, uz, sr, dr, c dintmn, dintmx, nint. c c pr Input Cylindrical radial coordinate of initial point on c track. Size np. c c sr Input Cylindrical 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 cylindrical radial component of unit direction c vector along track. Size np. c c ut Input Initial theta component of unit direction vector c along track. Theta is angle in x-y plane c counterclockwise from x axis. Size np. c c uz Input Initial axial z component of unit direction vector c 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 cylindrical surface. dimension sr (1) c---- Component r of direction vector. dimension ur (1) c---- Component theta of direction vector. dimension ut (1) c---- Component z of direction vector. dimension uz (1) c.... Local variables. c---- Coefficient of dint**2 in quadratic. common /laptrkcl/ aa (64) c---- Coefficient of dint in quadratic. common /laptrkcl/ bb (64) c---- A very big number. common /laptrkcl/ big c---- Constant term in quadratic. common /laptrkcl/ cc (64) c---- Distance to intersection 1. common /laptrkcl/ dint1 (64) c---- Distance to intersection 2. common /laptrkcl/ dint2 (64) c---- A very small number. common /laptrkcl/ fuz c---- 1 = tangent, 2 = bad qq. common /laptrkcl/ itrun (64) c---- Index in local array. common /laptrkcl/ n c---- First index of subset of data. common /laptrkcl/ n1 c---- Last index of subset of data. common /laptrkcl/ n2 c---- Index in external array. common /laptrkcl/ nn c---- 1 if root1 an exit. common /laptrkcl/ nroot1 c---- 1 if root2 an exit. common /laptrkcl/ nroot2 c---- Number of real roots of quadratic. common /laptrkcl/ nroots (64) c---- Size of current subset of data. common /laptrkcl/ ns c---- Factor bb**2 - 4.0 * aa * cc. common /laptrkcl/ qq (64) c---- New ur at intersection 1. common /laptrkcl/ ur1 c---- New ur at intersection 2. common /laptrkcl/ ur2 cbugc***DEBUG begins. cbug 9901 format (/ 'aptrkcl finding intersections. np=',i3, cbug & ' tol=',1pe13.5) cbug 9902 format (i3,' pr= ',1pe22.14 / cbug & ' ur,ut,uz= ',1p3e22.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), ut(n), uz(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) = ur(nn)**2 + ut(nn)**2 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), due to terms cancelling. qq(n) = 4.0 * (sr(nn)**2 * aa(n) - pr(nn)**2 * ut(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 nroots(n) = 1 endif if (nroot1 .eq. 1) then nint(nn) = 1 dint(nn) = dint1(n) 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 (/ 'aptrkcl 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 aptrkcl. (+1 line.) end UCRL-WEB-209832