subroutine aptwist (ax, ay, az, vx, vy, vz, pitch,
& px, py, pz, np, tol, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTWIST
c
c call aptwist (ax, ay, az, vx, vy, vz, pitch,
c & px, py, pz, np, tol, nerr)
c
c Version: aptwist Updated 1995 May 8 13:30.
c aptwist Originated 1995 May 8 13:30.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To twist the np points p = (px, py, pz) around the axis through
c point a = (ax, ay, az) parallel to vector v = (vx, vy, vz).
c For a point at axial distance d from point "a", the twist is
c clockwise 2 * pi * d / pitch (radians).
c Note: a clockwise twist looks counterclockwise when looking in
c the negative direction along the axis.
c Flag nerr indicates any input error.
c
c For example, points on a straight line parallel to the axis
c would be twisted to lie on a curve shaped like a coil spring.
c
c Calls: aptmovp, aptrotv, apttran, aptvusz
c
c Glossary:
c
c ax,ay,az Input Coordinates x, y and z of point "a" on the twist axis.
c
c nerr Output Indicates an input error, if not zero.
c 1 if vector "v" has zero length.
c 2 if pitch is zero.
c 3 if np is not positive.
c
c np Input Size of arrays px, py, pz. Number of points "p".
c
c pitch Input The pitch of the twist, i.e., the distance along the
c axis vector for one complete revolution around the
c axis.
c
c px,py,pz In/Out Coordinates x, y and z of point "p", to be twisted
c around twist axis "v" through point "a".
c
c tol Input Numerical tolerance limit.
c On Cray computers, recommend 1.e-5 to 1.e-11.
c
c vx,vy,vz Input Components x, y and z of twist axis vector "v".
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c____ Mathematical constant pi. Dimensionless.
parameter (pi = 3.14159265358979323)
c.... Dimensioned arguments.
dimension px(1) ! Coordinate x of point "p".
dimension py(1) ! Coordinate y of point "p".
dimension pz(1) ! Coordinate z of point "p".
c.... Local variables.
common /laptwist/ amx ! Negative of coordinate ax.
common /laptwist/ amy ! Negative of coordinate ay.
common /laptwist/ amz ! Negative of coordinate az.
common /laptwist/ angle ! Angle around twist axis.
common /laptwist/ cosang ! Cosine of angle around twist axis.
common /laptwist/ n ! Index of a point.
common /laptwist/ pxerr ! Error in calculating px.
common /laptwist/ pxold ! Coordinate px before twist.
common /laptwist/ pyerr ! Error in calculating py.
common /laptwist/ pyold ! Coordinate py before twist.
common /laptwist/ rotm (3,3) ! Rotation matrix
common /laptwist/ sinang ! Sine of angle around twist axis.
common /laptwist/ ux ! Component x of unit vector.
common /laptwist/ uy ! Component y of unit vector.
common /laptwist/ uz ! Component z of unit vector.
common /laptwist/ vlength ! Length of vector "v".
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptwist twisting points around an axis. np=',i3,
cbug & ' tol= ',1pe22.14 /
cbug & ' pitch= ',1pe22.14 /
cbug & ' ax,ay,az=',1p3e22.14 /
cbug & ' vx,vy,vz=',1p3e22.14 )
cbug 9902 format (i3,' px,py,pz=',1p3e22.14 )
cbug
cbug write ( 3, 9901) np, tol, pitch, ax, ay, az, vx, vy, vz
cbug do 110 n = 1, np
cbug write ( 3, 9902) n, px(n), py(n), pz(n)
cbug 110 continue
cbugc***DEBUG ends.
c.... initialize.
nerr = 0
c.... Test for input errors.
call aptvusz (vx, vy, vz, tol, ux, uy, uz, vlength)
if (vlength .eq. 0.0) then
nerr = 1
go to 200
endif
if (pitch .eq. 0.0) then
nerr = 2
go to 200
endif
if (np .le. 0) then
nerr = 3
go to 200
endif
c.... Do a translation on points "p" that translates point "a" to the origin.
call apttran (ax, ay, az, px, py, pz, np, tol, nerr)
c.... Do a rotation on points "p" that rotates vector "v" to the z axis.
call aptrotv (vx, vy, vz, 0.0, 0.0, 1.0, tol, rotm, nerr)
call aptmopv (rotm, 0, 0.0, 0.0, 0.0, px, py, pz, np, tol, nerr)
c.... Twist all of the points.
do 100 n = 1, np ! Loop over points.
c.... Find the angle of twist.
angle = 2.0 * pi * pz(n) / pitch
c.... Twist this point.
cosang = cos (angle)
if (abs (cosang - 1.0) .le. tol) cosang = 1.0
if (abs (cosang + 1.0) .le. tol) cosang = -1.0
if (abs (cosang) .le. tol) cosang = 0.0
sinang = sin (angle)
if (abs (sinang - 1.0) .le. tol) sinang = 1.0
if (abs (sinang + 1.0) .le. tol) sinang = -1.0
if (abs (sinang) .le. tol) sinang = 0.0
pxold = px(n)
pyold = py(n)
px(n) = pxold * cosang - pyold * sinang
pxerr = tol * (abs (pxold * cosang) + abs (pyold * sinang))
if (abs (px(n)) .le. pxerr) px(n) = 0.0
py(n) = pxold * sinang + pyold * cosang
pyerr = tol * (abs (pxold * sinang) + abs (pyold * cosang))
if (abs (py(n)) .le. pyerr) py(n) = 0.0
100 continue ! End of loop over points.
c.... Do a rotation on point "p" that rotates the z axis to vector "v".
call aptmopv (rotm, 1, 0.0, 0.0, 0.0, px, py, pz, np, tol, nerr)
if (nerr .ne. 0) nerr = 4
c.... Do a translation on point "p" that translates the origin to point "a".
amx = -ax
amy = -ay
amz = -az
call apttran (amx, amy, amz, px, py, pz, np, tol, nerr)
200 continue
cbugc***DEBUG begins.
cbug 9903 format (/ 'aptwist results: nerr=',i3)
cbug write ( 3, 9903) nerr
cbug do 205 n = 1, np
cbug write ( 3, 9902) n, px(n), py(n), pz(n)
cbug 205 continue
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptwist. (+1 line.)
end
UCRL-WEB-209832