subroutine aptsclr (scale, ax, ay, az, bx, by, bz, px, py, pz, np,
& tol, refm, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTSCLR
c
c call aptsclr (scale, ax, ay, az, bx, by, bz, px, py, pz, np,
c & tol, refm, nerr)
c
c Version: aptsclr Updated 1994 February 4 11:20.
c aptsclr Originated 1994 February 4 11:20.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To find the matrix operator refm for cylindrical scaling with
c the factor scale, in the direction perpendicular to the axis in
c the direction of vector a = (ax, ay, az), through the point
c b = (bx, by, bz), and to do the scaling on the np points or
c vectors p = (px, py, pz). np may be 0. If p = (px, py, pz) are
c unbound vectors, make sure invariant point "b" is at the origin.
c Any components of refm within tol of -1.0, 0.0, or 1.0, will be
c truncated to those values.
c Flag nerr indicates any input error.
c
c History:
c
c Input: scale, ax, ay, az, bx, by, bz, px, py, pz, np, tol.
c
c Output: px, py, pz, refm, nerr.
c
c Calls: aptvunb, aptmopv
c
c
c Glossary:
c
c ax,ay,az Input The x, y, z components of vector "a" in the direction
c of the axis of cylindrical scaling, which passes
c through invariant point "b".
c
c bx,by,bz Input The x, y, z coordinates of invariant point "b", on the
c axis of cylindrical scaling.
c
c nerr Output Indicates an input error, if not 0.
c 1 if the magnitude of vector "a" is too
c small, relative to tol.
c
c np Input Number of points or vectors "p". May be 0.
c
c px,py,pz In/Out The x, y, z coordinates or components of point or
c vector "p", before and after scaling. Size np.
c
c refm Output Cylindrical scaling operator (a unitary 3 x 3 matrix).
c Must be sized refm(3,3).
c
c scale Input Cylindrical scaling factor. A negative value is
c equivalent to a positive cylindrical scaling,
c followed by a
c reflection in the plane with the normal vector "a".
c
c tol Input Numerical tolerance limit. Used to test and adjust
c unit vector, matrix element, and point
c 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---- Cylindrical scaling matrix operator.
dimension refm (3,3)
c.... Local variables.
c---- Row index of cylindrical scaling matrix.
common /laptsclr/ i
c---- Column index of cylindrical scaling matrix.
common /laptsclr/ j
c---- Component x of unit normal vector.
common /laptsclr/ ux
c---- Component y of unit normal vector.
common /laptsclr/ uy
c---- Component z of unit normal vector.
common /laptsclr/ uz
c---- Magnitude of a vector.
common /laptsclr/ vlen
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptsclr. Cylindrical scaling by factor scale=',
cbug & 1pe22.14 /
cbug & ' perpendicular to the axis ax,ay,az=' / 5x,1p3e22.14 /
cbug & ' through point bx,by,bz=' / 5x,1p3e22.14)
cbug 9902 format (i3,2x,'px,py,pz=',1p3e22.14)
cbug write ( 3, 9901) scale, ax, ay, az, bx, by, bz
cbug if (np .gt. 0) then
cbug write ( 3, 9902) (n, px(n), py(n), pz(n), n = 1, np)
cbug endif
cbugc***DEBUG ends.
c.... Initialize.
nerr = 0
c.... Find the unit vector in the direction of the vector "a".
call aptvunb (ax, ay, az, 1, 0., ux, uy, uz, vlen, nerr)
c---- Vector magnitude too small.
if (vlen .le. tol) then
nerr = 1
go to 210
endif
c.... Form the components of the cylindrical scaling matrix.
refm(1,1) = scale + (1.0 - scale) * ux**2
refm(1,2) = + (1.0 - scale) * ux * uy
refm(1,3) = + (1.0 - scale) * ux * uz
refm(2,1) = + (1.0 - scale) * ux * uy
refm(2,2) = scale + (1.0 - scale) * uy**2
refm(2,3) = + (1.0 - scale) * uy * uz
refm(3,1) = + (1.0 - scale) * ux * uz
refm(3,2) = + (1.0 - scale) * uy * uz
refm(3,3) = scale + (1.0 - scale) * 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)) .le. tol) then
refm(i,j) = 0.0
elseif ((abs (abs (refm(i,j)) - 1.0)) .le. tol) then
refm(i,j) = sign (1.0, refm(i,j))
endif
110 continue
120 continue
c---- Tested tol.
endif
cbugc***DEBUG begins.
cbug 9903 format (' ux,uy,uz=',1p3e22.14)
cbug 9904 format (/ ' refm=',3(/ 1p3e22.14))
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 cylindrical scaling operation on the np points or vectors "p".
c---- Scale the points or vectors.
if (np .gt. 0) then
call aptmopv (refm, 0, bx, by, bz, px, py, pz, np, tol, nerr)
cbugc***DEBUG begins.
cbug write ( 3, 9902) (n, px(n), py(n), pz(n), n = 1, np)
cbugc***DEBUG ends.
c---- Tested np.
endif
210 return
c.... End of subroutine aptsclr. (+1 line.)
end
UCRL-WEB-209832