subroutine aptwirl (ax, ay, az, vx, vy, vz, pitch, rinv,
& px, py, pz, np, tol, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTWIRL
c
c call aptwirl (ax, ay, az, vx, vy, vz, pitch, rinv,
c & px, py, pz, np, tol, nerr)
c
c Version: aptwirl Updated 1995 May 8 13:30.
c aptwirl 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 radius r from the axis, the twist is clockwise
c 2 * pi * (r - rinv) / 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 perpendicular to the axis
c would be twisted to lie on an Archimedean spiral.
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 rinv is negative.
c 4 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 increment in radial
c distance from the axis for one complete revolution
c around the 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 rinv Input Radius from the twist axis at which the twist is zero,
c or an integer number of complete revolutions.
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 /laptwirl/ amx ! Negative of coordinate ax.
common /laptwirl/ amy ! Negative of coordinate ay.
common /laptwirl/ amz ! Negative of coordinate az.
common /laptwirl/ angle ! Angle around twist axis.
common /laptwirl/ cosang ! Cosine of angle around twist axis.
common /laptwirl/ n ! Index of a point.
common /laptwirl/ pxerr ! Error in calculating px.
common /laptwirl/ pxold ! Coordinate px before twist.
common /laptwirl/ pyerr ! Error in calculating py.
common /laptwirl/ pyold ! Coordinate py before twist.
common /laptwirl/ rotm (3,3) ! Rotation matrix
common /laptwirl/ sinang ! Sine of angle around twist axis.
common /laptwirl/ ux ! Component x of unit vector.
common /laptwirl/ uy ! Component y of unit vector.
common /laptwirl/ uz ! Component z of unit vector.
common /laptwirl/ vlength ! Length of vector "v".
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptwirl twisting points around an axis. np=',i3,
cbug & ' tol= ',1pe22.14 /
cbug & ' pitch= ',1pe22.14 /
cbug & ' rinv= ',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, rinv, 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 (rinv .lt. 0.0) then
nerr = 3
go to 200
endif
if (np .le. 0) then
nerr = 4
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 radial distance from the twist axis.
raxis = sqrt (px(n)**2 + py(n)**2)
c.... Find the angle of twist.
angle = 2.0 * pi * (raxis - rinv) / 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 (/ 'aptwirl 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 aptwirl. (+1 line.)
end
UCRL-WEB-209832