subroutine aptpers (ax, ay, az, bx, by, bz, cx, cy, cz,
& dx, dy, dz, px, py, pz, np, tol,
& qu, qv, dap, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTPERS
c
c call aptpers (ax, ay, az, bx, by, bz, cx, cy, cz,
c & dx, dy, dz, px, py, pz, np, tol,
c & qu, qv, dap, nerr)
c
c Version: aptpers Updated 1990 November 28 14:50.
c aptpers Originated 1990 October 1 10: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) toward the view point
c a = (ax, ay, az), onto the plane through the focal point
c b = (bx, by, bz) and perpendicular to the line "ab", and to
c find the distance dap from view point "a" to each point "p".
c The projected points will then be rotated and translated into
c the uv plane, so that the projection of the horizon line through
c the points c = (cx, cy, cz) and d = (dx, dy, dz) will be
c parallel to the x axis, and the origin is at point "b".
c The projection of each point "p" is point q = (qu, qv).
c No points "p" behind the view point "a" will be projected.
c
c Flag nerr indicates any input error.
c
c Input: ax, ay, az, bx, by, bz, cx, cy, cz, dx, dy, dz,
c px, py, pz, np, tol.
c
c Output: qu, qv, dap, nerr.
c
c Calls: aptmopv, aptrotc, aptrotv, apttran, aptvdic, aptvdis
c
c
c Glossary:
c
c ax,ay,az Input The x, y, z coordinates of the view point "a".
c
c bx,by,bz Input The x, y, z coordinates of focal point "b" in the
c projection plane. The projection plane is
c perpendicular to the line "ab", and will have
c its origin at point "b".
c Must be distinct from point "a".
c
c cx,cy,cz Input A point on the horizon of the view.
c Must be distinct from point "d".
c The u (horizontal) axis of the projection plane will
c be parallel to the projection of line "cd".
c
c dx,dy,dz Input A point on the horizon of the view.
c Must be distinct from point "c".
c The u (horizontal) axis of the projection plane will
c be parallel to the projection of line "cd".
c
c dap Output The distance from point "a" to point "p". Size np.
c
c nerr Output Indicates an input error, if not 0.
c 1 if np is not positive.
c 2 if the line "ab" is too short, based on tol.
c 3 if the projection of line "cd" is too short,
c based on tol.
c
c np Input Size of arrays px, py, pz, qu, qv, dap.
c
c px,py,pz Input The x, y, z coordinates of point "p". Size np.
c No points "p" behind the view point "a", relative to
c the direction "ab", will be projected.
c
c qu, qv Output The u and v coordinates of the projection of point
c "p" in the plane that passes through point "b" and
c is perpendicular to the line "ab". The u axis is
c parallel to the projection of line "cd" in the plane,
c and the origin is at point "b". Size np.
c For points "p" behind point "a", relative to the
c direction "ab", point "q" will have very large
c coordinates.
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---- Distance from point "a" to point "p".
dimension dap (1)
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 u of point "q".
dimension qu (1)
c---- Coordinate v of point "q".
dimension qv (1)
c.... Local variables.
c---- Component x of vector "ab".
common /laptpers/ abx
c---- Component y of vector "ab".
common /laptpers/ aby
c---- Component z of vector "ab".
common /laptpers/ abz
c---- A very big number.
common /laptpers/ big
c---- Modified value of cx.
common /laptpers/ ccx
c---- Modified value of cy.
common /laptpers/ ccy
c---- Modified value of cz.
common /laptpers/ ccz
c---- Component x of vector "cd".
common /laptpers/ cdx
c---- Component y of vector "cd".
common /laptpers/ cdy
c---- Distance between points "c" and "d".
common /laptpers/ dcd
c---- Modified value of dx.
common /laptpers/ ddx
c---- Modified value of dy.
common /laptpers/ ddy
c---- Modified value of dz.
common /laptpers/ ddz
c---- Distance between points "a" and "b".
common /laptpers/ dab
c---- A very small number.
common /laptpers/ fuz
c---- Index in arrays.
common /laptpers/ n
c---- First index of subset of data.
common /laptpers/ n1
c---- Last index of subset of data.
common /laptpers/ n2
c---- Index in external array.
common /laptpers/ nn
c---- Size of current subset of data.
common /laptpers/ ns
c---- Modified value of px.
common /laptpers/ ppx (64)
c---- Modified value of py.
common /laptpers/ ppy (64)
c---- Modified value of pz.
common /laptpers/ ppz (64)
c---- Rotation matrix.
common /laptpers/ rotm (3,3)
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptpers finding the view from point a:' /
cbug & ' ax,ay,az= ',1p3e22.14 /
cbug & ' looking in the direction of point b:' /
cbug & ' bx,by,bz= ',1p3e22.14 /
cbug & ' projected on the plane thru b perpendicular to ab.' /
cbug & ' The horizon is defined by points c and d:' /
cbug & ' cx,cy,cz=',1p3e22.14 /
cbug & " dx,dy,dz=",1p3e22.14 /
cbug & " The points viewed are:" /
cbug & (i3," px,py,pz=",1p3e22.14))
cbug write ( 3, 9901) ax, ay, az, bx, by, bz,
cbug & cx, cy, cz, dx, dy, dz,
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 vector "ab".
call aptvdis (ax, ay, az, bx, by, bz, 1, tol,
& abx, aby, abz, dab, nerr)
if (dab .le. tol) then
nerr = 2
cbugc***DEBUG begins.
cbug write ( 3, '("Points a and b coincide.")')
cbugc***DEBUG ends.
go to 210
endif
c.... Initialize temporary points "c" and "d".
ccx = cx
ccy = cy
ccz = cz
ddx = dx
ddy = dy
ddz = dz
c.... Move the origin to point "a".
call apttran (ax, ay, az, ccx, ccy, ccz, 1, tol, nerr)
call apttran (ax, ay, az, ddx, ddy, ddz, 1, tol, nerr)
c.... Rotate the vector "ab" to the z axis, rotate points "c" and "d".
call aptrotv (abx, aby, abz, 0., 0., 1., tol, rotm, nerr)
call aptmopv (rotm, 0, 0., 0., 0., ccx, ccy, ccz, 1, tol, nerr)
call aptmopv (rotm, 0, 0., 0., 0., ddx, ddy, ddz, 1, tol, nerr)
c.... Project points "c" and "d" onto the projection plane.
ccx = dab * ccx / (ccz + fuz)
ccy = dab * ccy / (ccz + fuz)
if (ccz .lt. 0.0) then
ccy = big
endif
ddx = dab * ddx / (ddz + fuz)
ddy = dab * ddy / (ddz + fuz)
if (ddz .lt. 0.0) then
ddy = big
endif
c.... Find the projected vector "cd".
call aptvdic (ccx, ccy, ddx, ddy, 1, tol, cdx, cdy, dcd, nerr)
if (dcd .le. tol) then
nerr = 3
cbugc***DEBUG begins.
cbug write ( 3, '("Points c and d coincide.")')
cbugc***DEBUG ends.
go to 210
endif
c.... Set up the indices of the first subset of data.
n1 = 1
n2 = min (np, 64)
110 ns = n2 - n1 + 1
c.... Initialize temporary points "p".
c---- Loop over subset of data.
do 120 n = 1, ns
nn = n + n1 - 1
ppx(n) = px(nn)
ppy(n) = py(nn)
ppz(n) = pz(nn)
c---- End of loop over subset of data.
120 continue
cbugc***DEBUG begins.
cbug 9801 format (/ 'aptpers temporaries (initial):' /
cbug & (i3,' ppx,y,z=',1p3e22.14))
cbug write ( 3, 9801) (n, ppx(n), ppy(n), ppz(n), n = 1, np)
cbugc***DEBUG ends.
c.... Move the origin to point "a".
call apttran (ax, ay, az, ppx, ppy, ppz, ns, tol, nerr)
c.... Rotate the vector "ab" to the z axis.
call aptmopv (rotm, 0, 0., 0., 0., ppx, ppy, ppz, ns, tol, nerr)
cbugc***DEBUG begins.
cbug 9802 format (/ 'aptpers temporaries (-a, ab onto z):' /
cbug & (i3,' ppx,y,z=',1p3e22.14))
cbug write ( 3, 9802) (n, ppx(n), ppy(n), ppz(n), n = 1, np)
cbugc***DEBUG ends.
c.... Find distance dap, and project points "p" onto the projection plane.
c---- Loop over subset of data.
do 130 n = 1, ns
nn = n + n1 - 1
dap(nn) = sqrt (ppx(n)**2 + ppy(n)**2 + ppz(n)**2)
ppx(n) = dab * ppx(n) / (ppz(n) + fuz)
if (ppz(n) .lt. 0.0) then
ppx(n) = big
endif
ppy(n) = dab * ppy(n) / (ppz(n) + fuz)
if (ppz(n) .lt. 0.0) then
ppy(n) = big
endif
c---- End of loop over subset of data.
130 continue
cbugc***DEBUG begins.
cbug 9803 format (/ 'aptpers temporaries (projected):' /
cbug & (i3,' ppx,y =',1p2e22.14))
cbug write ( 3, 9803) (n, ppx(n), ppy(n), n = 1, np)
cbugc***DEBUG ends.
c.... Rotate the x axis to be parallel to the projection of line "cd".
call aptrotc (cdx, cdy, 1., 0., ppx, ppy, ns, tol, nerr)
cbugc***DEBUG begins.
cbug 9806 format (/ 'aptpers temporaries (cd parallel to x):' /
cbug & (i3,' ppx,y =',1p2e22.14))
cbug write ( 3, 9806) (n, ppx(n), ppy(n), n = 1, np)
cbugc***DEBUG ends.
c.... Store the final coordinates of points "q".
c---- Loop over subset of data.
do 140 n = 1, ns
nn = n + n1 - 1
qu(nn) = ppx(n)
qv(nn) = -ppy(n)
c---- End of loop over subset of data.
140 continue
c.... See if all data subsets are done.
c---- Do another subset of data.
if (n2 .lt. np) then
n1 = n2 + 1
n2 = min (np, n1 + 63)
go to 110
endif
cbugc***DEBUG begins.
cbug 9902 format (/ 'aptpers results:' /
cbug & (i3,' qu,qv,ap= ',1p3e22.14))
cbug write ( 3, 9902) (n, qu(n), qv(n), dap(n), n = 1, np)
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptpers. (+1 line.)
end
UCRL-WEB-209832