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