subroutine aptrefs (itype, ax, ay, az, bx, by, bz, cx, cy, cz,
& px, py, pz, np, tol, refm, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTREFS
c
c call aptrefs (itype, ax, ay, az, bx, by, bz, cx, cy, cz,
c & px, py, pz, np, tol, refm, nerr)
c
c Version: aptrefs Updated 1990 March 14 16:00.
c aptrefs 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 the matrix operator refm for reflection from a plane:
c itype = 0: the plane contains point a = (ax, ay, az), and
c has normal vector b = (bx, by, bz).
c itype = 1: the reflection exchanges the points a = (ax, ay, az)
c and (bx, by, bz).
c itype = 2: the plane contains the 3 points a = (ax, ay, az),
c b = (bx, by, bz), and c = (cx, cy, cz).
c
c To do the reflection operation on the np points
c or vectors p = (px, py, pz). If p = (px, py, pz) are unbound
c vectors, either use module aptmopv, or make sure the reflection
c plane contains the origin. Size np may be 0.
c The new values of (px, py, pz) will be truncated to zero
c if less than the estimated error in their calculation,
c based on tol.
c Flag nerr indicates any input error, if not zero.
c
c History: 1990 March 14. Changed tol to 0.0 in call to unit vector
c subroutine. Allows small magnitudes.
c
c Input: itype, ax, ay, az, bx, by, bz, cx, cy, cz, px, py, pz, np, tol.
c
c Output: px, py, pz, refm, nerr.
c
c Calls: aptvxun, aptvunb
c
c
c Glossary:
c
c ax,ay,az Input The x, y, z coordinates or components of a point or
c vector.
c
c bx,by,bz Input The x, y, z coordinates or components of a point or
c vector.
c
c cx,cy,cz Input The x, y, z coordinates of a point.
c
c itype Input Defines option for describing reflection plane:
c 0 if plane contains point a = (ax, ay, az), and
c and has normal vector b = (bx, by, bz).
c 1 if reflection exchanges (ax, ay, az), (bx, by, bz).
c 2 if plane contains points (ax, ay, az), (bx, by, bz)
c and (cx, cy, cz).
c
c nerr Output Indicates an input error, if not 0.
c 1 if the magnitude of any input vector is too small,
c 2 if any 2 of the points (a, b, c) are congruent,
c or if the points (a, b, c) are colinear, and
c itype = 3.
c 3 if itype is not 0, 1 or 2.
c
c np Input Number of points or vectors (px, py, pz). May be 0.
c
c px,py,pz In/Out The x, y, z coordinates of a point, or
c the x, y, z components of a vector. Size np.
c Will be truncated to zero, if less than the estimated
c error in their calculation, based on tol.
c
c refm Output Reflection operator (a unitary 3 x 3 matrix).
c Must be sized refm(3,3).
c
c tol Input Numerical tolerance limit. Used to test and adjust
c unit vector and point components.
c On Cray computers, recommend 1.e-5 to 1.e-11.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c---- Coordinate or component x. Size np.
dimension px (1)
c---- Coordinate or component y. Size np.
dimension py (1)
c---- Coordinate or component z. Size np.
dimension pz (1)
c---- Reflection matrix operator.
dimension refm (3,3)
c.... Local variables.
c---- Row index of reflection matrix.
common /laptrefs/ i
c---- Column index of reflection matrix.
common /laptrefs/ j
c---- Coordinate x of point in plane.
common /laptrefs/ qx
c---- Coordinate y of point in plane.
common /laptrefs/ qy
c---- Coordinate z of point in plane.
common /laptrefs/ qz
c---- Component x of unit normal vector.
common /laptrefs/ ux
c---- Component y of unit normal vector.
common /laptrefs/ uy
c---- Component z of unit normal vector.
common /laptrefs/ uz
c---- Magnitude of a vector.
common /laptrefs/ vlen
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptrefs. Reflecting in plane defined by itype=',i2 /
cbug & ' 0 for a=pt in plane, b=normal, c=null.' /
cbug & ' 1 for a and b exchanged, c=null.' /
cbug & ' 2 for a, b and c in plane.' /
cbug & ' tol=',1pe13.5 /
cbug & ' ax,ay,az=',1p3e22.14 /
cbug & ' bx,by,bz=',1p3e22.14 /
cbug & " cx,cy,cz=",1p3e22.14)
cbug write ( 3, 9901) itype, tol, ax, ay, az, bx, by, bz, cx, cy, cz
cbugc***DEBUG ends.
c.... Initialize.
nerr = 0
c.... Find the unit vector normal to the plane.
c---- Plane through "a", normal is "b".
if (itype .eq. 0) then
c.... Find the unit vector for vector "b".
call aptvunb (bx, by, bz, 1, 0., ux, uy, uz, vlen, nerr)
c---- Vector magnitude too small.
if (vlen .le. tol) then
nerr = 1
go to 210
endif
qx = ax
qy = ay
qz = az
c---- Reflection exchanges "a", "b".
elseif (itype .eq. 1) then
c.... Find the unit vector for vector "b - a".
call aptvunb (bx - ax, by - ay, bz - az, 1, 0.,
& ux, uy, uz, vlen, nerr)
c---- Points "a" and "b" are congruent.
if (vlen .le. tol) then
nerr = 1
go to 210
endif
qx = 0.5 * (ax + bx)
qy = 0.5 * (ay + by)
qz = 0.5 * (az + bz)
c---- Plane is through "a", "b", "c".
elseif (itype .eq. 2) then
c.... Find the unit vector of the cross product of b - a and c - a.
call aptvxun (bx - ax, by - ay, bz - az,
& cx - ax, cy - ay, cz - az, 1, tol,
& ux, uy, uz, vlen, nerr)
c---- 2 or 3 points are congruent.
if (vlen .le. tol) then
nerr = 1
go to 210
endif
qx = ax
qy = ay
qz = az
c---- Itype is not 0, 1, or 2.
else
nerr = 2
go to 210
c---- Tested itype.
endif
c.... Form the components of the reflection matrix.
c.... Row and column vectors are non-orthogonal unit vectors.
refm(1,1) = 1.0 - 2.0 * ux**2
refm(1,2) = - 2.0 * ux * uy
refm(1,3) = - 2.0 * ux * uz
refm(2,1) = - 2.0 * ux * uy
refm(2,2) = 1.0 - 2.0 * uy**2
refm(2,3) = - 2.0 * uy * uz
refm(3,1) = - 2.0 * ux * uz
refm(3,2) = - 2.0 * uy * uz
refm(3,3) = 1.0 - 2.0 * uz**2
c---- Adjust components near 0 or 1.
if (tol .gt. 0.0) then
do 120 i = 1, 3
do 110 j = 1, 3
if (abs (refm(i,j)) .lt. tol) then
refm(i,j) = 0.0
elseif ((abs (abs (refm(i,j)) - 1.0)) .lt. tol) then
refm(i,j) = sign (1.0, refm(i,j))
endif
110 continue
120 continue
c---- Tested tol.
endif
c.... See if the reflection plane contains the origin.
vlen = sqrt (qx**2 + qy**2 + qz**2)
if (abs (qx * ux + qy * uy + qz * uz) .le. tol * vlen) then
qx = 0.0
qy = 0.0
qz = 0.0
endif
cbugc***DEBUG begins.
cbug 9902 format (/ ' invariant point:' /
cbug & ' qx,qy,qz=',1p3e22.14)
cbug 9903 format (' ux,uy,uz=',1p3e22.14)
cbug 9904 format (/ ' refm=',3(/ 11x,1p3e22.14))
cbug write ( 3, 9902) qx, qy, qz
cbug write ( 3, 9903) ux, uy, uz
cbug write ( 3, 9904) ((refm(i,j), j = 1, 3), i = 1, 3)
cbugc***DEBUG ends.
c.... Do the reflection operation on the np points or vectors (px, py, pz).
c.... Truncate very small coordinates or components to zero.
c---- Reflect the points or vectors.
if (np .gt. 0) then
call aptmopv (refm, 0, qx, qy, qz, px, py, pz, np, tol, nerr)
c---- Tested np.
endif
210 return
c.... End of subroutine aptrefs. (+1 line.)
end
UCRL-WEB-209832