subroutine aptrefl (pu, pv, pw, au, av, bu, bv, np, tol,
& qu, qv, qw, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTREFL
c
c call aptrefl (pu, pv, pw, au, av, bu, bv, np, tol,
c & qu, qv, qw, nerr)
c
c Version: aptrefl Updated 1990 November 29 15:30.
c aptrefl 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 vector
c q = (qu, qv, qw) resulting from the reflection of vector
c p = (pu, pv, pw) from the surface perpendicular to the uv plane
c and through the points a = (au, av, 0), and b = (bu, bv, 0),
c when incident on the surface on the line "ab".
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, v and w 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, pw, au, av, bu, bv, np, tol.
c
c Output: qu, qv, qw, 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,pw Input The u, v, w components of incident vector "p".
c Size np.
c
c qu,qv,qw Output The u, v, w components of 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 Note that qw = pw. 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 w of incident vector "p".
dimension pw (1)
c---- Component u of reflected vector "q".
dimension qu (1)
c---- Component v of reflected vector "q".
dimension qv (1)
c---- Component w of reflected vector "q".
dimension qw (1)
c.... Local variables.
c---- Ratio of magnitudes of "p", "q".
common /laptrefl/ fact
c---- Factor 2.0 * pdots / vlen2.
common /laptrefl/ fan
c---- A very small number.
common /laptrefl/ fuz
c---- First index of subset of data.
common /laptrefl/ n1
c---- Last index of subset of data.
common /laptrefl/ n2
c---- Index in external array.
common /laptrefl/ nn
c---- Size of current subset of data.
common /laptrefl/ ns
c---- Dot product of "p" and "sn".
common /laptrefl/ pdots
c---- Estimated error in pdots.
common /laptrefl/ pdserr
c---- Estimated error in qu.
common /laptrefl/ querr
c---- Estimated error in qv.
common /laptrefl/ qverr
c---- Component u of normal vector.
common /laptrefl/ snu
c---- Estimated error in snu.
common /laptrefl/ snuerr
c---- Component v of normal vector.
common /laptrefl/ snv
c---- Estimated error in snv.
common /laptrefl/ snverr
c---- Square of magnitude of vector "sn".
common /laptrefl/ vlen2
c---- Square of magnitude of vector "p".
common /laptrefl/ vlenp2
c---- Square of magnitude of vector "q".
common /laptrefl/ vlenq2
cbugc***DEBUG begins.
cbugc---- Index in arrays.
cbug common /laptrefl/ n
cbug 9901 format (/ 'aptrefl finding reflection of the vector',
cbug & ' (tol=',1pe13.5,'):' /
cbug & (i3,' pu,pv,pw=',1p3e22.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), pw(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 vector "p" in the plane with normal
c.... vector "sn", which is 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)
qw(nn) = pw(nn)
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 (av(nn)) + abs (bv(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 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
qw(nn) = pw(nn)
c.... Make the magnitude of vector "q" the same as vector "p".
vlenp2 = pu(nn)**2 + pv(nn)**2 + pw(nn)**2
vlenq2 = qu(nn)**2 + qv(nn)**2 + qw(nn)**2
fact = sqrt ((vlenp2 + fuz) / (vlenq2 + fuz))
qu(nn) = fact * qu(nn)
qv(nn) = fact * qv(nn)
qw(nn) = fact * qw(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 (/ 'aptrefl results:' /
cbug & (i3,' qu,qv,qw=',1p3e22.14))
cbug write ( 3, 9902) (n, qu(n), qv(n), qw(n), n = 1, np)
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptrefl. (+1 line.)
end
UCRL-WEB-209832