subroutine aptrkcy (pz, pr, uz, ur, ut, az, ar, bz, br, & dintmn, dintmx, np, noptd, tolf, tols, & nint, pinz, pinr, dint, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTRKCY c c call aptrkcy (pz, pr, uz, ur, ut, az, ar, bz, br, c & dintmn, dintmx, np, noptd, tolf, tols, c & nint, pinz, pinr, dint, nerr) c c Version: aptrkcy Updated 1993 December 6 15:20. c aptrkcy Originated 1989 December 7 16:40. 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 acceptable c exit intersection point pin = (pinz, pinr) of the linear track c through point p = (pz, pr) with the unit direction vector c u = (uz, ur, ut), with the conical surface symmetric around the c z axis, and through the points a = (az, ar) and b = (bz, br), c for which: (1) the distance from point "p" to point "pin" is c between the limits dintmn and dintmx, (2) the intersection is c between points "a" and "b", and (3) the crossing at the c intersection is from left to right in the z-r plane (right to c left in the r-z plane). This is a zone exit if the points "a" c and "b" are the vertices of a zone edge, counterclockwise around c the zone in the the z-r plane (clockwise in the r-z plane). c Flag nerr indicates any input error, if not zero. c c The coordinate system is 3-D cylindrical, with the axial c coordinate z, the radial coordinate r = sqrt (x**2 + y**2), c and the angular coordinate theta (counterclockwise in the x-y c plane from the x axis. c c History: 1990 January 23 10:10. Implemented use of aptfdav to adjust c values of fint1 and fint2 near 0.0 or 1.0. Vectorized the c intersection test loop. c c Input: pz, pr, uz, ur, ut, az, ar, bz, br, dintmn, dintmx, np, c noptd, tolf, tols. c c Output: nint, pinz, pinr, dint, nerr. c c Calls: aptqrtv, aptvlim, aptfdav c c c Glossary: c c az, ar Input Coordinates z and r of the beginning of the surface. c Size np. c c bz, br Input Coordinates z and r of the end of the surface. c Size np. c c dint Output Distance from point "p" to the intersection at point c "pin" along the track (not distance in z-r plane). c Size np. c c dintmn Input Minimum allowed value of the distance to intersection c dint. Size np. c c dintmx Input Maximum allowed value of the distance to intersection c dint. 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 acceptable exit intersection was found. c 1 if an acceptable exit intersection was found on the c surface segment between points "a" and "b", with a c distance to the intersection dint between dintmn and c dintmx. Size np. c c noptd Input Option to limit the minimum magnitude of the components c of the direction vector u = (uz, ur, ut). The method c used here will fail if uz or ut is zero. c 0 for no magnitude test. Must be done elsewhere. c 1 to limit the minimum magnitudes of ut and uz to no c less than tols. c c np Input Size of arrays pz, pr, uz, ur, ut, az, ar, bz, br, c dintmn, dintmx, nint, pinz, pinr. c c pinr Output Coordinate r of the intersection point "pin". Size np. c c pinz Output Coordinate z of the intersection point "pin". Size np. c c pz, pr Input Coordinates z and r of the initial point on the linear c track. Size np. c c uz,ur,ut Input Initial z, r and theta components of the unit direction c vector along the linear track. Size np. c c tolf Input Truncation error limit to be imposed on the fractional c distance of the intersection point "pin" along c the line segment from point "a" to point "b". Values c less than -tolf or greater than 1.0 + tolf will not c be accepted. Values from -tolf to tolf will be c changed to tolf. Values from 1.0 - tolf to c 1.0 + tolf will be changed to 1.0 - tolf. c Also used to test for the intersection being nearly c tangent, and for the accuracy of the intersection. c Must not be zero. c On Cray computers, recommend 1.e-11. c c tols Input Truncation error limit to be imposed on uz, ut and c pinr. Magnitudes of ut and uz less than tols c will be increased to tols. Values of pinr less c than tols * pr will be increased to tols * pr. c On Cray computers, recommend 1.e-5. c A value of zero may produce unpredictable results. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccend. c.... Dimensioned arguments. c---- Coordinate r at start of edge "ab". dimension ar (1) c---- Coordinate z at start of edge "ab". dimension az (1) c---- Coordinate r at end of edge "ab". dimension br (1) c---- Coordinate z at end of edge "ab". dimension bz (1) c---- Distance to intersection "pin". dimension dint (1) c---- Minimum distance to intersection. dimension dintmn (1) c---- Maximum distance to intersection. dimension dintmx (1) c---- Number of intersections (0 or 1). dimension nint (1) c---- Coordinate r of intersection "pin". dimension pinr (1) c---- Coordinate z of intersection "pin". dimension pinz (1) c---- Coordinate r of beginning of track. dimension pr (1) c---- Coordinate z of beginning of track. dimension pz (1) c---- Component r of direction vector "u". dimension ur (1) c---- Component theta of "u". dimension ut (1) c---- Component z of direction vector "u". dimension uz (1) c.... Local variables. c---- Coefficient of fint**2 in quadratic. common /laptrkcy/ aa (64) c---- Coefficient of fint in quadratic. common /laptrkcy/ bb (64) c---- A very big number. common /laptrkcy/ big c---- Factor 2.0 * pr * ur / uz. common /laptrkcy/ c1 (64) c---- Factor (ur**2 + ut**2) / uz**2. common /laptrkcy/ c2 (64) c---- Constant term in quadratic. common /laptrkcy/ cc (64) c---- Distance to intersection 1. common /laptrkcy/ dint1 c---- Distance to intersection 2. common /laptrkcy/ dint2 c---- Difference br - ar. common /laptrkcy/ dr (64) c---- Difference bz - az. common /laptrkcy/ dz (64) c---- First root of quadratic. common /laptrkcy/ fint1 (64) c---- Second root of quadratic. common /laptrkcy/ fint2 (64) c---- A very small number. common /laptrkcy/ fuz c---- 1 = tangent, 2 = bad qq. common /laptrkcy/ itrun (64) c---- Index in local array. common /laptrkcy/ n c---- First index of subset of data. common /laptrkcy/ n1 c---- Last index of subset of data. common /laptrkcy/ n2 c---- 1 if intersection 1 good exit. common /laptrkcy/ nexit1 c---- 1 if intersection 2 good exit. common /laptrkcy/ nexit2 c---- 1 if fint1 or fint 2 adjusted. common /laptrkcy/ nlim (64) c---- Index in external array. common /laptrkcy/ nn c---- Number of real roots of quadratic. common /laptrkcy/ nroots (64) c---- Size of current subset of data. common /laptrkcy/ ns c---- Coordinate r at intersection 1. common /laptrkcy/ pinr1 c---- Coordinate r at intersection 2. common /laptrkcy/ pinr2 c---- Coordinate z at intersection 1. common /laptrkcy/ pinz1 c---- Coordinate z at intersection 2. common /laptrkcy/ pinz2 c---- Factor bb**2 - 4.0 * aa * cc. common /laptrkcy/ qq (64) c---- New ur at intersection 1. common /laptrkcy/ ur1 c---- New ur at intersection 2. common /laptrkcy/ ur2 c---- Magnitude of direction vector. common /laptrkcy/ vlen (64) cbugc***DEBUG begins. cbug 9901 format (/ 'aptrkcy finding intersections. np=',i3, cbug & ' noptd=',i2 / ' tolf=',1pe22.14,' tols=',1pe22.14) cbug 9902 format (i3,' pz,pr= ',1p2e22.14 / cbug & ' uz,ur,ut= ',1p3e22.14 / cbug & ' az,ar= ',1p2e22.14 / cbug & ' bz,br= ',1p2e22.14 / cbug & ' dintmn= ',1pe22.14,' dintmx=',1pe22.14) cbug write ( 3, 9901) np, noptd, tolf, tols cbug write ( 3, 9902) (nn, pz(nn), pr(nn), uz(nn), ur(nn), ut(nn), cbug & az(nn), ar(nn), bz(nn), br(nn), dintmn(nn), dintmx(nn), cbug & nn = 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 c.... Test for input errors. nerr = 0 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.... Do not allow any tracks that are linear in the z-r plane. c---- Limit minimum size of uz, ut. if (noptd .ne. 0) then call aptvlim (uz(n1), ur(n1), ut(n1), ns, tols, 0., tols, & vlen, nerr) c---- Tested noptd. endif c.... Find the coefficients of aa * fint**2 + bb * fint + cc = 0, where c.... fint = ((pinz - az) * dz + (pinr - ar) * dr) / (dr**2 + dz**2), c.... by substituting pinz = az + fint * dz, pinr = ar + fint * dr, in c.... pinr = sqrt (pr**2 + c1 * (pinz - pz) + c2 * (pinz - pz)**2). c---- Loop over local arrays. do 120 n = 1, ns nn = n + n1 - 1 c1(n) = 2.0 * pr(nn) * ur(nn) / uz(nn) c2(n) = (ur(nn)**2 + ut(nn)**2) / uz(nn)**2 dr(n) = br(nn) - ar(nn) dz(n) = bz(nn) - az(nn) aa(n) = c2(n) * dz(n)**2 - dr(n)**2 bb(n) = (c1(n) + 2.0 * c2(n) * (az(nn) - pz(nn))) * dz(n) - & 2.0 * ar(nn) * dr(n) cc(n) = pr(nn)**2 - ar(nn)**2 + (az(nn) - pz(nn)) * & (c1(n) + c2(n) * (az(nn) - pz(nn))) 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 * dr(n) * (dr(n) * & (pr(nn)**2 + (az(nn) - pz(nn)) * & (c1(n) + c2(n) * (az(nn) - pz(nn)))) - & dz(n) * ar(nn) * (c1(n) + & 2.0 * c2(n) * (az(nn) - pz(nn)))) + dz(n)**2 * & (c1(n)**2 - 4.0 * c2(n) * (pr(nn)**2 - ar(nn)**2)) c---- End of loop over local arrays. 120 continue c.... Find any intersections. call aptqrtv (1, aa, bb, cc, qq, ns, tolf, & nroots, fint1, fint2, itrun, nerr) c.... Adjust any values of fint1 and fint2 near 0.0 or 1.0. call aptfdav (fint1, ns, 1, tolf, nlim, nerr) call aptfdav (fint2, ns, 1, tolf, nlim, nerr) c.... See if either intersection is on zone edge, on track, c.... and an exit from zone. c---- Loop over local arrays. do 130 n = 1, ns nn = n + n1 - 1 pinz1 = az(nn) + fint1(n) * dz(n) pinz2 = az(nn) + fint2(n) * dz(n) pinr1 = amax1 (pr(nn) * tols, ar(nn) + fint1(n) * dr(n)) pinr2 = amax1 (pr(nn) * tols, ar(nn) + fint2(n) * dr(n)) dint1 = (pinz1 - pz(nn)) / uz(nn) dint2 = (pinz2 - pz(nn)) / uz(nn) ur1 = (pr(nn) * ur(nn) + (ur(nn)**2 + ut(nn)**2) * dint1) / & (pinr1 + fuz) ur2 = (pr(nn) * ur(nn) + (ur(nn)**2 + ut(nn)**2) * dint2) / & (pinr2 + fuz) if ((nroots(n) .ge. 1) .and. & (fint1(n) .ge. 0.0) .and. & (fint1(n) .le. 1.0) .and. & (dint1 .ge. dintmn(nn)) .and. & (dint1 .le. dintmx(nn)) .and. & ((uz(nn) * dr(n)) .gt. (ur1 * dz(n)))) then nexit1 = 1 else nexit1 = 0 endif if (nexit1 .eq. 1) then nroots(n) = 1 endif if (nexit1 .eq. 1) then nint(nn) = 1 pinr(nn) = pinr1 pinz(nn) = pinz1 dint(nn) = dint1 else nint(nn) = 0 pinr(nn) = -big pinz(nn) = -big dint(nn) = -big endif if ((nroots(n) .eq. 2) .and. & (fint2(n) .ge. 0.0) .and. & (fint2(n) .le. 1.0) .and. & (dint2 .ge. dintmn(nn)) .and. & (dint2 .le. dintmx(nn)) .and. & ((uz(nn) * dr(n)) .gt. (ur2 * dz(n)))) then nexit2 = 1 else nexit2 = 0 endif if (nexit2 .eq. 1) then nint(nn) = 1 pinr(nn) = pinr2 pinz(nn) = pinz2 dint(nn) = dint2 endif c---- End of loop over local arrays. 130 continue cbugc***DEBUG begins. cbug 9903 format (/ 'aptrkcy results:') cbug 9904 format (i3,' nint=',i2 / cbug & ' z,r,dint= ',1p3e22.14 / cbug & ' nn=',i3 / cbug & ' fint(1,2)=',1p2e22.14 / cbug & ' c1,c2,qq= ',1p3e22.14 / cbug & ' aa,bb,cc= ',1p3e22.14) cbug write ( 3, 9903) cbug do 140 n = 1, ns cbug nn = n + n1 - 1 cbug write ( 3, 9904) nn, nint(nn), pinz(nn), pinr(nn), dint(nn), cbug & n, fint1(n), fint2(n), cbug & c1(n), c2(n), qq(n), aa(n), bb(n), cc(n) cbug 140 continue 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 aptrkcy. (+1 line.) end UCRL-WEB-209832