subroutine aptprcy (nopt, ax, ay, az, bx, by, bz, rcyl,
& px, py, pz, np, tol, qx, qy, qz, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTPRCY
c
c call aptprcy (nopt, ax, ay, az, bx, by, bz, rcyl,
c & px, py, pz, np, tol, qx, qy, qz, nerr)
c
c Version: aptprcy Updated 1990 November 29 10:50.
c aptprcy Originated 1990 October 3 11:40.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To project the np points p = (px, py, pz) onto the cylindrical
c surface which has an axis passing through point a = (ax, ay, az)
c in the direction of the vector b = (bx, by, bz), and a radius
c of rcyl. The projection option is specified by nopt.
c The projection of each point "p" is point q = (qx, qy, qz).
c
c Flag nerr indicates any input error.
c
c Input: nopt, ax, ay, az, bx, by, bz, rcyl, px, py, pz, np, tol.
c
c Output: qx, qy, qz, nerr.
c
c Calls: aptmopv, aptrotv, apttran, aptvunb
c
c
c Glossary:
c
c ax,ay,az Input The x, y, z coordinates of a point on the axis of
c the cylindrical projection surface.
c
c bx,by,bz Input The x, y, z components of the vector parallel to the
c axis of the cylindrical projection surface.
c
c nerr Output Indicates an input error, if not 0.
c 1 if np is not positive.
c 2 if the vector "b" is too short, based on tol.
c 3 if nopt is not 0, 1 or 2.
c
c nopt Input Each vector "ap", from point "a" to each point "p", has
c a component in the radial direction from the axis of
c the cylinder and a component in the axial direction,
c relative to point "a". The projection options are:
c 0 to project in the cylindrical radial direction.
c This gives the radial component a length rcyl, and
c leaves the axial component unchanged.
c 1 to project in the spherical radial direction, based
c on point "a". This gives the radial component a
c length rcyl, and multiplies the axial component by
c the ratio of rcyl to the initial radial component.
c 2 to project the radial component in the cylindrical
c radial direction, giving it a length rcyl, and to
c replace the axial component with the cosine of the
c angle formed by the axis, point "a", and point "p".
c This is a projection which preserves relative areas
c when projecting points from a sphere centered at
c point "a" onto the cylinder.
c
c np Input Size of arrays px, py, pz, qx, qy, qz.
c The number of points projected.
c
c px,py,pz Input The x, y, z coordinates of point "p". Size np.
c
c qx,qy,qz Output The x, y, z coordinates of the projection of the point
c "p" on the projection surface, which is a cylinder of
c radius rcyl with an axis through point "a", parallel
c to the vector "b". For any point "p" which is on
c the axis of the cylinder, based on tol, the
c coordinates of point "q" will be very large.
c Size np.
c
c rcyl Input The radius of the cylindrical projection surface,
c measured from the axis, which passes through point
c "a" in the direction of vector "b".
c
c tol Input Numerical tolerance limit.
c On Cray computers, recommend 1.e-5 to 1.e-11.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
c---- Coordinate x of point "p".
dimension px (1)
c---- Coordinate y of point "p".
dimension py (1)
c---- Coordinate z of point "p".
dimension pz (1)
c---- Coordinate x of point "q".
dimension qx (1)
c---- Coordinate y of point "q".
dimension qy (1)
c---- Coordinate z of point "q".
dimension qz (1)
c.... Local variables.
c---- A very big number.
common /laptprcy/ big
c---- Distance between points "a" and "p".
common /laptprcy/ dap
c---- Distance from the axis to point "p".
common /laptprcy/ daxis
c---- A very small number.
common /laptprcy/ fuz
c---- Index in arrays.
common /laptprcy/ n
c---- Rotation matrix.
common /laptprcy/ rotm (3,3)
c---- Component x of unit vector "ub".
common /laptprcy/ ubx
c---- Component y of unit vector "ub".
common /laptprcy/ uby
c---- Component z of unit vector "ub".
common /laptprcy/ ubz
c---- Length of axis vector "b".
common /laptprcy/ vlenb
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptprcy projecting points onto a cylinder:' /
cbug & ' rcyl=',1pe22.14,' nopt=',i2 /
cbug & ' The axis goes through the point:' /
cbug & ' ax,ay,az= ',1p3e22.14 /
cbug & ' parallel to the vector:' /
cbug & ' bx,by,bz= ',1p3e22.14 /
cbug & ' The points projected are:' /
cbug & (i3," px,py,pz=",1p3e22.14))
cbug write ( 3, 9901) rcyl, nopt, ax, ay, az, bx, by, bz,
cbug & (n, px(n), py(n), pz(n), n = 1, np)
cbugc***DEBUG ends.
c.... Initialize.
c---- A very big number.
big = 1.e+99
c---- A very small number.
fuz = 1.e-99
nerr = 0
c.... Test for input errors.
if (np .le. 0) then
nerr = 1
cbugc***DEBUG begins.
cbug write ( 3, '(/ "np = 0")')
cbugc***DEBUG ends.
go to 210
endif
c.... Find the length of the axis vector "b".
call aptvunb (bx, by, bz, 1, 0., ubx, uby, ubz, vlenb, nerr)
if (vlenb .le. tol) then
nerr = 2
cbugc***DEBUG begins.
cbug write ( 3, '(/ "Vector b is too short.")')
cbugc***DEBUG ends.
go to 210
endif
if ((nopt .lt. 0) .or. (nopt .gt. 2)) then
nerr = 3
cbugc***DEBUG begins.
cbug write ( 3, '(/ "nop not between 0 and 2")')
cbugc***DEBUG ends.
go to 210
endif
c.... Find the rotation matrix to rotate vector "b" onto the z axis.
call aptrotv (bx, by, bz, 0., 0., 1., tol, rotm, nerr)
c.... Initialize the projected points "q".
c---- Loop over data.
do 110 n = 1, np
qx(n) = px(n)
qy(n) = py(n)
qz(n) = pz(n)
c---- End of loop over data.
110 continue
c.... Move the origin to point "a".
call apttran (ax, ay, az, qx, qy, qz, np, tol, nerr)
c.... Rotate the vector "b" to the z axis.
call aptmopv (rotm, 0, 0., 0., 0., qx, qy, qz, np, tol, nerr)
c.... Find the projection option.
if (nopt .eq. 0) then
c.... Project in the cylindrical radial direction.
c.... Does not change the axial component, relative to point "a".
c---- Loop over data.
do 120 n = 1, np
daxis = sqrt (qx(n)**2 + qy(n)**2)
qx(n) = rcyl * qx(n) / (daxis + fuz)
qy(n) = rcyl * qy(n) / (daxis + fuz)
if (daxis .le. tol) then
qx(n) = big
qy(n) = big
endif
c---- End of loop over data.
120 continue
elseif (nopt .eq. 1) then
c.... Project in the spherical radial direction, based on point "a".
c---- Loop over data.
do 130 n = 1, np
daxis = sqrt (qx(n)**2 + qy(n)**2)
qx(n) = rcyl * qx(n) / (daxis + fuz)
qy(n) = rcyl * qy(n) / (daxis + fuz)
qz(n) = rcyl * qz(n) / (daxis + fuz)
if (daxis .le. tol) then
qx(n) = big
qy(n) = big
qz(n) = big
endif
c---- End of loop over data.
130 continue
elseif (nopt .eq. 2) then
c.... Project the radial components in the cylindrical radial direction.
c.... Replace the axial component with the cosine of the angle formed
c.... by the axis, point "a", and point "p".
c---- Loop over data.
do 140 n = 1, np
daxis = sqrt (qx(n)**2 + qy(n)**2)
dap = sqrt (qx(n)**2 + qy(n)**2 + qz(n)**2)
qx(n) = rcyl * qx(n) / (daxis + fuz)
qy(n) = rcyl * qy(n) / (daxis + fuz)
qz(n) = qz(n) / (dap + fuz)
if (daxis .le. tol) then
qx(n) = big
qy(n) = big
endif
if (dap .le. tol) then
qz(n) = big
endif
c---- End of loop over data.
140 continue
c---- Tested nopt.
endif
c.... Rotate the z axis back to vector "b".
call aptmopv (rotm, 1, 0., 0., 0., qx, qy, qz, np, tol, nerr)
c.... Restore the original origin.
call apttran (-ax, -ay, -az, qx, qy, qz, np, tol, nerr)
cbugc***DEBUG begins.
cbug 9902 format (/ 'aptprcy results:' /
cbug & (i3,' qx,qy,qz=',1p3e22.14))
cbug write ( 3, 9902) (n, qx(n), qy(n), qz(n), n = 1, np)
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptprcy. (+1 line.)
end
UCRL-WEB-209832