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