subroutine aptprsp (nopt, ax, ay, az, bx, by, bz, rsph,
& px, py, pz, np, tol, qx, qy, qz, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTPRSP
c
c call aptprsp (nopt, ax, ay, az, bx, by, bz, rsph,
c & px, py, pz, np, tol, qx, qy, qz, nerr)
c
c Version: aptprsp Updated 1990 November 29 10:50.
c aptprsp Originated 1990 October 4 17:30.
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 spherical
c surface centered at the point a = (ax, ay, az), with the polar
c axis in the direction of the vector b = (bx, by, bz), and with
c radius rsph. 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, rsph, 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 the center of the sphere.
c
c bx,by,bz Input The x, y, z components of the vector parallel to the
c polar axis of the sphere.
c
c nerr Output Indicates an input error, if not 0.
c 1 if np is not positive.
c 2 if nopt is not 0, 1 or 2.
c 3 if rsph is not positive.
c 4 if the vector "b" is too short, based on tol,
c and nopt is not 0.
c
c nopt Input Each vector "ap", from point "a" to each point "p", has
c a radial component "apr" perpendicular to the polar
c axis, and an axial component "apz" parallel to the
c polar axis. The projection options are:
c 0 to project in the spherical radial direction.
c This makes the vector "aq" parallel to "ap", with a
c length rsph. Any point "p" at point "a" is
c projected to very large coordinates.
c 1 to project in the direction perpendicular to the
c polar axis. This does not change the vector "apz",
c except that magnitudes greater than rsph are
c reduced to rsph, putting the point "q" on a pole.
c Any point "p" on the polar axis is projected to
c very large coordinates.
c 2 to project in the axial direction. This does not
c change the vector "apr", except that magnitudes
c greater than rsph are reduced to rsph, putting the
c point "q" on the equator. Any point "p" on the
c equatorial plane is projected to very large
c coordinates.
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 sphere
c centered at point "a", with a polar axis parallel to
c vector "b", and a radius rsph. The coordinates of
c of point "q" will be very large for any point "p"
c which, based on tol, is:
c For nopt = 0, at point "a".
c For nopt = 1, on the polar axis.
c For nopt = 2, on the equatorial plane.
c Size np.
c
c rsph Input The radius of the spherical projection surface,
c measured from point "a". Must be positive.
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---- Argument of square root.
common /laptprsp/ arg
c---- A very big number.
common /laptprsp/ big
c---- Distance between points "a" and "p".
common /laptprsp/ dap
c---- Distance from the axis to point "p".
common /laptprsp/ daxis
c---- Distance from the equatorial plane.
common /laptprsp/ deq
c---- Projection factor.
common /laptprsp/ fact
c---- A very small number.
common /laptprsp/ fuz
c---- Index in arrays.
common /laptprsp/ n
c---- Rotation matrix.
common /laptprsp/ rotm (3,3)
c---- Component x of unit vector "ub".
common /laptprsp/ ubx
c---- Component y of unit vector "ub".
common /laptprsp/ uby
c---- Component z of unit vector "ub".
common /laptprsp/ ubz
c---- Length of axis vector "b".
common /laptprsp/ vlenb
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptprsp projecting points onto a sphere:' /
cbug & ' rsph=',1pe22.14,' nopt=',i2,' tol=',1pe13.5 /
cbug & ' The center is at the point:' /
cbug & ' ax,ay,az= ',1p3e22.14 /
cbug & ' the polar axis is 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) rsph, nopt, tol, 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++++ Disallowed value of nopt.
if ((nopt .lt. 0) .or. (nopt .gt. 2)) then
nerr = 2
cbugc***DEBUG begins.
cbug write ( 3, '(/ "nop not between 0 and 2")')
cbugc***DEBUG ends.
go to 210
c---- Tested nopt.
endif
if (rsph .le. 0.0) then
nerr = 3
cbugc***DEBUG begins.
cbug write ( 3, '(/ "rsph .le. 0.0")')
cbugc***DEBUG ends.
go to 210
endif
if (nopt .ne. 0) then
c.... Find the length of the polar axis vector "b".
call aptvunb (bx, by, bz, 1, 0., ubx, uby, ubz, vlenb, nerr)
c---- Vector "b" is too short.
if (vlenb .le. tol) then
nerr = 4
cbugc***DEBUG begins.
cbug write ( 3, '(/ "Vector b is too short.")')
cbugc***DEBUG ends.
go to 210
c---- Tested vlenb.
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---- Tested nopt.
endif
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)
if (nopt .ne. 0) then
c.... Rotate the polar axis vector "b" to the z axis.
call aptmopv (rotm, 0, 0., 0., 0., qx, qy, qz, np, tol, nerr)
c---- Tested nopt.
endif
c.... Find the projection option.
if (nopt .eq. 0) then
c.... Project in the spherical radial direction.
c---- Loop over data.
do 120 n = 1, np
dap = sqrt (qx(n)**2 + qy(n)**2 + qz(n)**2)
qx(n) = rsph * qx(n) / (dap + fuz)
qy(n) = rsph * qy(n) / (dap + fuz)
qz(n) = rsph * qz(n) / (dap + fuz)
if (dap .le. tol) then
qx(n) = big
qy(n) = big
qz(n) = big
endif
c---- End of loop over data.
120 continue
elseif (nopt .eq. 1) then
c.... Project in the direction perpendicular to the polar axis.
c.... Put points projected outside the sphere onto the poles.
c---- Loop over data.
do 130 n = 1, np
daxis = sqrt (qx(n)**2 + qy(n)**2)
arg = amax1 (0.0, rsph**2 - qz(n)**2)
fact = sqrt (arg) / (daxis + fuz)
qx(n) = fact * qx(n)
qy(n) = fact * qy(n)
qz(n) = amin1 (qz(n), rsph)
qz(n) = amax1 (qz(n), -rsph)
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 in the direction of the polar axis.
c.... Put points projected outside the sphere onto the equator.
c---- Loop over data.
do 140 n = 1, np
daxis = sqrt (qx(n)**2 + qy(n)**2)
arg = amax1 (0.0, rsph**2 - daxis**2)
fact = amin1 (rsph / (daxis + fuz), 1.0)
deq = abs (qz(n))
qx(n) = fact * qx(n)
qy(n) = fact * qy(n)
qz(n) = sign (sqrt (arg), qz(n))
if (deq .le. tol) then ! Do not know which way to go.
qx(n) = big
qy(n) = big
qz(n) = big
endif
c---- End of loop over data.
140 continue
c---- Tested nopt.
endif
if (nopt .ne. 0) then
c.... Restore the initial z axis.
call aptmopv (rotm, 1, 0., 0., 0., qx, qy, qz, np, tol, nerr)
c---- Tested nopt.
endif
c.... Restore the original origin.
call apttran (-ax, -ay, -az, qx, qy, qz, np, tol, nerr)
cbugc***DEBUG begins.
cbug 9902 format (/ 'aptprsp 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 aptprsp. (+1 line.)
end
UCRL-WEB-209832