subroutine aptripv (ax, ay, az, bx, by, bz, cx, cy, cz, np, tol,
& dx, dy, dz, dlength, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTRIPV
c
c call aptripv (ax, ay, az, bx, by, bz, cx, cy, cz, np, tol,
c & dx, dy, dz, dlength, nerr)
c
c Version: aptripv Updated 1995 November 28 14:30.
c aptripv Originated 1995 November 28 14:30.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: For each of np vector triples a = (ax, ay, az), b = (bx, by, bz)
c and c = (cx, cy, cz), to find the vector triple product
c d = (dx, dy, dz) and its length, dlength, where
c d = a cross (b cross c) = (c cross b) cross a =
c (a dot c) b - (a dot b) c.
c Vector "d" is perpendicular to "a", in the plane of "b" and "c".
c The values of dx, dy and dz will be truncated to zero, if less
c than the estimated error in their calculation, based on tol.
c Flag nerr indicates any input error.
c
c Input: ax, ay, az, bx, by, bz, cx, cy, cz, np, tol.
c
c Output: dx, dy, dz, dlength, nerr.
c
c Glossary:
c
c ax,ay,az Input The x, y, z components of input vector "a". Size np.
c
c bx,by,bz Input The x, y, z components of input vector "b". Size np.
c
c cx,cy,cz Input The x, y, z components of input vector "c". Size np.
c
c dx,dy,dz Output The vector triple product of vectors "a", "b" and "c",
c a cross (b cross c) = (c cross b) cross a. 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 dx, dy, dz, dlength.
c
c tol Input Numerical tolerance limit.
c For 64-bit floating point arithmetic, recommend
c 1.e-5 to 1.e-11.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Dimensioned arguments.
dimension ax (1) ! Component x of input vector "a".
dimension ay (1) ! Component y of input vector "a".
dimension az (1) ! Component z of input vector "a".
dimension bx (1) ! Component x of input vector "b".
dimension by (1) ! Component y of input vector "b".
dimension bz (1) ! Component z of input vector "b".
dimension cx (1) ! Component x of input vector "c".
dimension cy (1) ! Component y of input vector "c".
dimension cz (1) ! Component z of input vector "c".
dimension dx (1) ! Component x of output vector "d".
dimension dy (1) ! Component y of output vector "d".
dimension dz (1) ! Component z of output vector "d".
dimension dlength (1) ! Length of output vector "d".
c.... Local variables.
common /laptripv/ n ! Index, 1 to np.
common /laptripv/ dxerr ! Estimated error in dx.
common /laptripv/ dyerr ! Estimated error in dy.
common /laptripv/ dzerr ! Estimated error in dz.
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptripv finding vector triple product of vectors:' /
cbug & ' nopt=',i2 /
cbug & (i3,' ax,ay,az=',1p3e22.14 /
cbug & ' bx,by,bz=',1p3e22.14 /
cbug & ' cx,cy,cz=',1p3e22.14 ))
cbug write ( 3, 9901) nopt, (n, ax(n), ay(n), az(n),
cbug & bx(n), by(n), bz(n), cx(n), cy(n), cz(n), n = 1, np)
cbugc***DEBUG ends.
c.... Initialize.
nerr = 0
c.... Test for input errors.
if (np .le. 0) then
nerr = 1
go to 210
endif
c.... Find the vector triple products.
do 110 n = 1, np ! Loop over vectors.
dx(n) = ay(n) * (bx(n) * cy(n) - cx(n) * by(n)) -
& az(n) * (bz(n) * cx(n) - cz(n) * bx(n))
dy(n) = az(n) * (by(n) * cz(n) - cy(n) * bz(n)) -
& ax(n) * (bx(n) * cy(n) - cx(n) * by(n))
dz(n) = ax(n) * (bz(n) * cx(n) - cz(n) * bx(n)) -
& ay(n) * (by(n) * cz(n) - cy(n) * bz(n))
dlength(n) = sqrt (dx(n)**2 + dy(n)**2 + dz(n)**2)
110 continue ! End of loop over vectors.
c.... See if result should be tested for truncation to zero.
if (tol .gt. 0.0) then ! Truncate small result to zero.
do 120 n = 1, np ! Loop over vectors.
dxerr = 3.0 * tol *
& (abs (ay(n)) * (abs (bx(n) * cy(n)) +
& abs (cx(n) * by(n))) +
& abs (az(n)) * (abs (bz(n) * cx(n)) +
& abs (cz(n) * bx(n))))
dyerr = 3.0 * tol *
& (abs (az(n)) * (abs (by(n) * cz(n)) +
& abs (cy(n) * bz(n))) +
& abs (ax(n)) * (abs (bx(n) * cy(n)) +
& abs (cx(n) * by(n))))
dzerr = 3.0 * tol *
& (abs (ax(n)) * (abs (bz(n) * cx(n)) +
& abs (cz(n) * bx(n))) +
& abs (ay(n)) * (abs (by(n) * cz(n)) +
& abs (cy(n) * bz(n))))
if (abs (dx(n)) .le. dxerr) dx(n) = 0.0
if (abs (dy(n)) .le. dyerr) dy(n) = 0.0
if (abs (dz(n)) .le. dzerr) dz(n) = 0.0
dlength(n) = sqrt (dx(n)**2 + dy(n)**2 + dz(n)**2)
120 continue ! End of loop over vectors.
endif ! Tested tol.
cbugc***DEBUG begins.
cbug 9902 format (/ 'aptripv results:' /
cbug & (i3,' dx,dy,dz=',1p3e22.14 /
cbug & ' dlength= ',1pe22.14 ))
cbug write ( 3, 9902) (n, dx(n), dy(n), dz(n), dlength(n),
cbug & n = 1, np)
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptripv. (+1 line.)
end
UCRL-WEB-209832