subroutine aptrefc (pu, pv, au, av, bu, bv, np, tol, qu, qv, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTREFC
c
c call aptrefc (pu, pv, au, av, bu, bv, np, tol, qu, qv, nerr)
c
c Version: aptrefc Updated 1990 November 29 15:30.
c aptrefc Originated 1990 January 10 13:50.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To find, for each of np sets of input data, the 2-D vector
c q = (qu, qv) resulting from the reflection of the 2-D vector
c p = (pu, pv) from the line through the points a = (au, av)
c and b = (bu, bv), all in the uv plane.
c When point "a" coincides with point "b", based on tol,
c no reflection takes place, and vector "q" equals vector "p".
c Coordinates u and v may be any orthogonal coordinates.
c Flag nerr indicates any input error.
c
c Method: The vector "sn" normal to the reflecting surface is the cross
c product of the vector normal to the uv plane, (0, 0, 1), and
c the vector parallel to the line "ab". From simple geometry,
c "q" = "p" - 2.0 * ("p" dot "u") * "u", where "u" is the unit
c vector parallel to "sn".
c
c Input: pu, pv, au, av, bu, bv, np, tol.
c
c Output: qu, qv, nerr.
c
c Glossary:
c
c au, av Input The u and v coordinates of point "a" in the uv plane.
c Must differ from "b", based on tol. Size np.
c
c bu, bv Input The u and v coordinates of point "b" in the uv plane.
c Must differ from "a", based on tol. Size np.
c
c nerr Output Indicates an input error, if not 0.
c 1 if np is not positive.
c
c np Input Size of arrays.
c
c pu, pv Input The u and v components of the incident vector "p".
c Size np.
c
c qu, qv Output The u and v components of the reflected vector "q".
c Components of "q" less than the estimated error in
c their calculation, based on tol, will be truncated
c to zero.
c Size np.
c
c tol Input Numerical tolerance limit.
c On Cray computers, recommend 1.e-5 to 1.e-11.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c---- Coordinate u of point "a".
dimension au (1)
c---- Coordinate v of point "a".
dimension av (1)
c---- Coordinate u of point "b".
dimension bu (1)
c---- Coordinate v of point "b".
dimension bv (1)
c---- Component u of incident vector "p".
dimension pu (1)
c---- Component v of incident vector "p".
dimension pv (1)
c---- Component u of reflected vector "q".
dimension qu (1)
c---- Component v of reflected vector "q".
dimension qv (1)
c.... Local variables.
c---- Ratio of magnitudes of "p", "q".
common /laptrefc/ fact
c---- Factor 2.0 * pdots / vlen2.
common /laptrefc/ fan
c---- A very small number.
common /laptrefc/ fuz
c---- First index of subset of data.
common /laptrefc/ n1
c---- Last index of subset of data.
common /laptrefc/ n2
c---- Index in external array.
common /laptrefc/ nn
c---- Size of current subset of data.
common /laptrefc/ ns
c---- Dot product of "p" and "sn".
common /laptrefc/ pdots
c---- Estimated error in pdots.
common /laptrefc/ pdserr
c---- Estimated error in qu.
common /laptrefc/ querr
c---- Estimated error in qv.
common /laptrefc/ qverr
c---- Component u of normal vector.
common /laptrefc/ snu
c---- Estimated error in snu.
common /laptrefc/ snuerr
c---- Component v of normal vector.
common /laptrefc/ snv
c---- Estimated error in snv.
common /laptrefc/ snverr
c---- Magnitude of normal vector "sn".
common /laptrefc/ vlen2
c---- Square of magnitude of vector "p".
common /laptrefc/ vlenp2
c---- Square of magnitude of vector "q".
common /laptrefc/ vlenq2
cbugc***DEBUG begins.
cbugc---- Index in arrays.
cbug common /laptrefc/ n
cbug 9901 format (/ 'aptrefc finding reflection of the vector',
cbug & ' (tol=',1pe13.5,'):' /
cbug & (i3,' pu,pv=',1p2e22.14 /
cbug & ' from the line (in the uv plane) through:' /
cbug & ' au,av=',1p2e22.14 /
cbug & ' bu,bv=',1p2e22.14))
cbug write ( 3, 9901) tol, (n, pu(n), pv(n),
cbug & au(n), av(n), bu(n), bv(n), n = 1, np)
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
c.... Set up the indices of the first subset of data.
n1 = 1
n2 = min (np, 64)
110 ns = n2 - n1 + 1
c.... Find the reflection of incident 2-D vector "p", in the plane with
c.... normal vector "sn" perpendicular to line "ab" and the uv plane.
c---- No truncation tests.
if (tol .le. 0.0) then
c---- Loop over subset of data.
do 120 nn = n1, n2
c.... Find the vector "sn" normal to the reflecting surface.
snu = av(nn) - bv(nn)
snv = bu(nn) - au(nn)
vlen2 = snu**2 + snv**2
pdots = pu(nn) * snu + pv(nn) * snv
c.... Find the reflected vector "q".
qu(nn) = pu(nn) - 2.0 * pdots * snu / (vlen2 + fuz)
qv(nn) = pv(nn) - 2.0 * pdots * snv / (vlen2 + fuz)
c---- End of loop over subset of data.
120 continue
c---- Test for numerical truncation.
else
c---- Loop over subset of data.
do 130 nn = n1, n2
c.... Find the vector "sn" normal to the reflecting surface.
snu = av(nn) - bv(nn)
snuerr = tol * (abs (av(nn)) + abs (bv(nn)))
if (abs (snu) .lt. snuerr) then
snu = 0.0
endif
snv = bu(nn) - au(nn)
snverr = tol * (abs (au(nn)) + abs (bu(nn)))
if (abs (snv) .lt. snverr) then
snv = 0.0
endif
vlen2 = snu**2 + snv**2
pdots = pu(nn) * snu + pv(nn) * snv
pdserr = tol * (abs (pu(nn) * snu) + abs (pv(nn) * snv))
if (abs (pdots) .lt. pdserr) then
pdots = 0.0
endif
fan = 2.0 * pdots / (vlen2 + fuz)
c.... Find the reflected 2-D vector "q".
qu(nn) = pu(nn) - fan * snu
querr = tol * (abs (pu(nn)) + abs (fan * snu))
if (abs (qu(nn)) .lt. querr) then
qu(nn) = 0.0
endif
qv(nn) = pv(nn) - fan * snv
qverr = tol * (abs (pv(nn)) + abs (fan * snv))
if (abs (qv(nn)) .lt. qverr) then
qv(nn) = 0.0
endif
c.... Make the magnitude of 2-D vector "q" the same as 2-D vector "p".
vlenp2 = pu(nn)**2 + pv(nn)**2
vlenq2 = qu(nn)**2 + qv(nn)**2
fact = sqrt ((vlenp2 + fuz) / (vlenq2 + fuz))
qu(nn) = fact * qu(nn)
qv(nn) = fact * qv(nn)
c---- End of loop over subset of data.
130 continue
c---- Tested tol.
endif
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 9902 format (/ 'aptrefc results:' /
cbug & (i3,' qu,qv=',1p2e22.14))
cbug write ( 3, 9902) (n, qu(n), qv(n), n = 1, np)
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptrefc. (+1 line.)
end
UCRL-WEB-209832