subroutine aptradv (ax, ay, az, bx, by, bz, px, py, pz, np, tol,
& qu, qv, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTRADV
c
c call aptradv (ax, ay, az, bx, by, bz, px, py, pz, np, tol,
c & qu, qv, nerr)
c
c Version: aptradv Updated 1990 November 29 15:30.
c aptradv Originated 1990 October 2 10:30.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To do a geometric view factor projection of the np points
c p = (px, py, pz), relative to the differential surface element
c centered at point a = (ax, ay, az) with normal vector
c b = (bx, by, bz). First, each point "p" is projected radially
c onto the unit sphere centered at point "a". Then, each point
c on the unit sphere is projected perpendicularly onto the plane
c of the unit circle centered at point "a" with normal vector "b".
c
c The projection of each point "p" is point q = (qu, qv).
c The orientation of the "u" axis in the projection plane is
c arbitrary. Any points "p" behind the projection plane will
c be assigned very large values of qu and qv.
c
c If the points "p" bound a surface, then the ratio of the
c projected area to the area of the circle is the differential
c geometric view factor of the surface, as seen from the
c differential surface element at point "a" with normal vector
c "b".
c
c Flag nerr indicates any input error.
c
c Input: ax, ay, az, bx, by, bz, px, py, pz, np, tol.
c
c Output: qu, qv, 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 view point "a", in the
c surface with normal vector "b". The center of the
c unit sphere and unit circle of the projection.
c
c bx,by,bz Input The x, y, z components of the vector normal to the
c surface through point "a".
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
c np Input Size of arrays px, py, pz, qu, qv.
c
c px,py,pz Input The x, y, z coordinates of point "p". Size np.
c
c qu, qv Output The u and v coordinates of the geometric view factor
c projection of point "p", relative to the diffential
c surface element at point "a" with normal vector "b".
c The orientation of the "u" axis is arbitrary.
c For points "p" behind point "a", relative to the
c direction "b", point "q" will have very large
c coordinates. Size np.
c
c tol Input Numerical tolerance limit. With 64-bit floating
c point numbers, 1.e-5 to 1.e-11 is recommended.
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 u of point "q".
dimension qu (1)
c---- Coordinate v of point "q".
dimension qv (1)
c.... Local variables.
c---- A very big number.
common /laptradv/ big
c---- Distance between points "a" and "p".
common /laptradv/ dap
c---- A very small number.
common /laptradv/ fuz
c---- Index in arrays.
common /laptradv/ n
c---- First index of subset of data.
common /laptradv/ n1
c---- Last index of subset of data.
common /laptradv/ n2
c---- Index in external array.
common /laptradv/ nn
c---- Size of current subset of data.
common /laptradv/ ns
c---- Modified value of px.
common /laptradv/ ppx (64)
c---- Modified value of py.
common /laptradv/ ppy (64)
c---- Modified value of pz.
common /laptradv/ ppz (64)
c---- Rotation matrix.
common /laptradv/ rotm (3,3)
c---- Component x of unit vector "ub".
common /laptradv/ ubx
c---- Component y of unit vector "ub".
common /laptradv/ uby
c---- Component z of unit vector "ub".
common /laptradv/ ubz
c---- Length of normal vector "b".
common /laptradv/ vlenb
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptradv finding the F-view from point a:' /
cbug & ' ax,ay,az= ',1p3e22.14 /
cbug & ' with normal vector b:' /
cbug & ' bx,by,bz= ',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 & (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 normal 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
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.... 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 (/ 'aptradv 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 "b" to the z axis.
call aptmopv (rotm, 0, 0., 0., 0., ppx, ppy, ppz, ns, tol, nerr)
cbugc***DEBUG begins.
cbug 9802 format (/ 'aptradv 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.... Do a geometric view factor projection.
c---- Loop over subset of data.
do 130 n = 1, ns
dap = sqrt (ppx(n)**2 + ppy(n)**2 + ppz(n)**2)
ppx(n) = ppx(n) / (dap + fuz)
ppy(n) = ppy(n) / (dap + fuz)
nn = n + n1 - 1
if (ppz(n) .ge. 0.0) then
qu(nn) = ppy(n)
qv(nn) = ppx(n)
else
qu(nn) = big
qv(nn) = big
endif
c---- End of loop over subset of data.
130 continue
cbugc***DEBUG begins.
cbug 9803 format (/ 'aptradv temporaries (projected):' /
cbug & (i3,' ppx,ppy=',1p2e22.14))
cbug write ( 3, 9803) (n, ppx(n), ppy(n), n = 1, np)
cbugc***DEBUG ends.
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 (/ 'aptradv results:' /
cbug & (i3,' qu,qv =',1p2e22.14))
cbug write ( 3, 9902) (n, qu(n), qv(n), n = 1, np)
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptradv. (+1 line.)
end
UCRL-WEB-209832