subroutine aptbang (ax, ay, az, bx, by, bz, cx, cy, cz, & np, tol, bdx, bdy, bdz, dx, dy, dz, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTBANG c c call aptbang (ax, ay, az, bx, by, bz, cx, cy, cz, c & np, tol, bdx, bdy, bdz, dx, dy, dz, nerr) c c Version: aptbang Updated 1990 November 27 14:00. c aptbang Originated 1990 March 8 17:40. c c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123. c c c Purpose: To find, for each of np sets of input data, the bisector c bd = (bdx, bdy, bdz) of the angle "abc" formed by the points c a = (ax, ay, az), b = (bx, by, bz), and c = (cx, cy, cz), c and point d = (dx, dy, dz), the intercept of the bisector on c the line "ca". If points "a", "b" and "c" are colinear, c vector "bd" will be zero, and point "d" will be point "b". c Flag nerr indicates any input error, if not zero. c c History: 1990 March 30. Fixed array index error which affected problems c with np .gt. 64. c c Input: ax, ay, az, bx, by, bz, cx, cy, cz, np, tol. c c Output: bdx, bdy, bdz, dx, dy, dz, nerr. c c Calls: aptvdis, aptvuna c c c Glossary: c c ax,ay,az Input The x, y, z coordinates of point "a". Size np. c Size np. c c bdx,y,z Output The x, y, z components of the vector "bd" which c bisects angle "abc", and connects points "b" and "d". c Size np. c c bx,by,bz Input The x, y, z coordinates of point "b". Size np. c c cx,cy,cz Input The x, y, z coordinates of point "c". Size np. c c dx,dy,dz Output The x, y, z coordinates of point "d" on line "ca". c The intercept of bisector "bd" on line "ca". c Size np. c c nerr Output Indicates an input error, if not 0. c 1 if np is not positive. c c np Input Size of arrays ax, ay, az, bx, by, bz, cx, cy, cz, c bdx, bdy, bdz, dx, dy, dz. 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 "a". dimension ax (1) c---- Coordinate y of point "a". dimension ay (1) c---- Coordinate z of point "a". dimension az (1) c---- Component x of vector "bd". dimension bdx (1) c---- Component y of vector "bd". dimension bdy (1) c---- Component z of vector "bd". dimension bdz (1) c---- Coordinate x of point "b". dimension bx (1) c---- Coordinate y of point "b". dimension by (1) c---- Coordinate z of point "b". dimension bz (1) c---- Coordinate x of point "c". dimension cx (1) c---- Coordinate y of point "c". dimension cy (1) c---- Coordinate z of point "c". dimension cz (1) c---- Coordinate x of point "d". dimension dx (1) c---- Coordinate y of point "d". dimension dy (1) c---- Coordinate z of point "d". dimension dz (1) c.... Local variables. c---- Distance from "a" to "b". common /laptbang/ distab (64) c---- Distance from "b" to "c". common /laptbang/ distbc (64) c---- Estimated error in dx. common /laptbang/ dxerr c---- Estimated error in dy. common /laptbang/ dyerr c---- Estimated error in dz. common /laptbang/ dzerr c---- Temporary factor. common /laptbang/ fact c---- A very small number. common /laptbang/ fuz c---- Index in arrays. common /laptbang/ n c---- First index of subset of data. common /laptbang/ n1 c---- Last index of subset of data. common /laptbang/ n2 c---- Index in external array. common /laptbang/ nn c---- Size of current subset of data. common /laptbang/ ns c---- Component x of unit vector along "ab". common /laptbang/ uabx (64) c---- Component y of unit vector along "ab". common /laptbang/ uaby (64) c---- Component z of unit vector along "ab". common /laptbang/ uabz (64) c---- Component x of unit vector along "bc". common /laptbang/ ubcx (64) c---- Component y of unit vector along "bc". common /laptbang/ ubcy (64) c---- Component z of unit vector along "bc". common /laptbang/ ubcz (64) c---- Magnitude of a vector. common /laptbang/ vlen (64) c---- Component x of "uab" - "ubc". common /laptbang/ vx (64) c---- Component y of "uab" - "ubc". common /laptbang/ vy (64) c---- Component z of "uab" - "ubc". common /laptbang/ vz (64) cbugc***DEBUG begins. cbug 9901 format (/ 'aptbang finding bisector of angle abc. tol=',1pe22.14) cbug 9902 format (i3,' ax,ay,az=',1p3e22.14 / cbug & ' bx,by,bz=',1p3e22.14 / cbug & ' cx,cy,cz=',1p3e22.14) cbug write ( 3, 9901) tol cbug write ( 3, 9902) (n, ax(n), ay(n), az(n), bx(n), by(n), bz(n), cbug & cx(n), cy(n), cz(n), n = 1, np) cbugc***DEBUG ends. c.... Initialize. c---- A very small number. fuz = 1.e-99 nerr = 0 c.... Test for input errors. if (np .le. 0) then nerr = 1 go to 210 endif c.... Set up the indices of the first subset of data. n1 = 1 n2 = min (np, 64) c.... Loop over subsets of data. 110 ns = n2 - n1 + 1 c.... Find the unit vector "uab" along line "ab", distance distab. call aptvdis (ax(n1), ay(n1), az(n1), & bx(n1), by(n1), bz(n1), ns, tol, & uabx, uaby, uabz, distab, nerr) call aptvuna (uabx, uaby, uabz, ns, 0., vlen, nerr) c.... Find the unit vector "ubc" along line "bc", distance distbc. call aptvdis (bx(n1), by(n1), bz(n1), & cx(n1), cy(n1), cz(n1), ns, tol, & ubcx, ubcy, ubcz, distbc, nerr) call aptvuna (ubcx, ubcy, ubcz, ns, 0., vlen, nerr) c.... Find the difference between vectors "uab" and "ubc". call aptvdis (uabx, uaby, uabz, & ubcx, ubcy, ubcz, ns, tol, & vx, vy, vz, vlen, nerr) c.... Find the angle bisector "bd". c---- Loop over subset of data. do 120 n = 1, ns nn = n + n1 - 1 fact = distab(n) * distbc(n) / & (distab(n) + distbc(n) + fuz) bdx(nn) = fact * vx(n) bdy(nn) = fact * vy(n) bdz(nn) = fact * vz(n) c---- End of loop over subset of data. 120 continue c.... Find the intercept point "d". c---- Loop over subset of data. do 130 n = 1, ns nn = n + n1 - 1 fact = distab(n) + distbc(n) + fuz dx(nn) = (distab(n) * cx(nn) + distbc(n) * ax(nn)) / fact dy(nn) = (distab(n) * cy(nn) + distbc(n) * ay(nn)) / fact dz(nn) = (distab(n) * cz(nn) + distbc(n) * az(nn)) / fact c---- End of loop over subset of data. 130 continue c.... See if truncation to zero is allowed. c---- Truncation is allowed. if (tol .gt. 0.0) then c---- Loop over subset of data. do 140 n = 1, ns nn = n + n1 - 1 fact = distab(n) + distbc(n) + fuz dxerr = tol * (distab(n) * abs (cx(nn)) + & distbc(n) * abs (ax(nn))) / fact dyerr = tol * (distab(n) * abs (cy(nn)) + & distbc(n) * abs (ay(nn))) / fact dzerr = tol * (distab(n) * abs (cz(nn)) + & distbc(n) * abs (az(nn))) / fact if (abs (dx(nn)) .lt. dxerr) then dx(nn) = 0.0 endif if (abs (dy(nn)) .lt. dyerr) then dy(nn) = 0.0 endif if (abs (dz(nn)) .lt. dzerr) then dz(nn) = 0.0 endif c---- End of loop over subset of data. 140 continue c---- Tested tol. endif 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 9903 format (/ 'aptbang results:' / cbug & (i3,' bdx,y,z=',1p3e22.14 / cbug & ' dx,dy,dz=',1p3e22.14)) cbug write ( 3, 9903) (n, bdx(n), bdy(n), bdz(n), cbug & dx(n), dy(n), dz(n), n = 1, np) cbugc***DEBUG ends. 210 return c.... End of subroutine aptbang. (+1 line.) end UCRL-WEB-209832