subroutine aptetru (ax, ay, az, bx, by, bz, iunit, anga, angb,
& tol, ux, uy, uz, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTETRU
c
c call aptetru (ax, ay, az, bx, by, bz, iunit, anga, angb,
c & tol, ux, uy, uz, nerr)
c
c Version: aptetru Updated 2001 March 28 11:20.
c aptetru Originated 2001 March 28 11:20.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: Given vectors a = (ax, ay, az) and b = (bx, by, bz), to find
c unit vector u = (ux, uy, uz) at angle anga from vector "a" and
c at angle angb from vector "b". Angles anga and angb are in
c degrees (iunit = 0) or in radians (iunit = 1). These vectors
c may represent the edge vectors at the vertex of a polyhedron.
c
c No solution is possible if vectors "a" and "b" are parallel,
c or if the angles anga and angb are too close to zero or 180
c degrees (zero or pi radians), relative to the angle between
c vectors "a" and "b".
c
c Otherwise, two solutions are possible, with vector "u" on either
c side of a plane containing vectors "a" and "b". The solution
c returned here is that for which the vectors "a", "b" and "u"
c firm a positive triple, like the thumb, first finger and middle
c finger of the right hand, respectively. To get the other
c solution, exchange vectors "a" and "b" and exchange angles anga
c and angb.
c
c Flag nerr indicates any input error.
c
c Input: anga, angb, ax, ay, az, bx, by, bz, iunit, tol.
c
c Output: ux, uy, uz, nerr.
c
c Calls: aptvdot, aptvuna, aptvunb, aptvxun
c
c Glossary:
c
c anga Input The angle between vector "a" and vector "u",
c in degrees (iunit = 0) or in radians (iunit = 1).
c
c angb Input The angle between vector "b" and vector "u",
c in degrees (iunit = 0) or in radians (iunit = 1).
c
c ax,ay,az Input The x, y, z components of vector "a".
c
c bx,by,bz Input The x, y, z components of vector "b".
c
c nerr Output Indicates an input error, if not 0.
c 1 if vector "a" has zero length.
c 2 if vector "b" has zero length.
c 3 if vectors "a" and "b" are parallel.
c 4 if angles anga and angb allow no solution.
c
c ux,uy,uz Output The x, y, z components of vector "u".
c Vectors "a", "b" and "u" form a positive triple,
c i.e,: "a" cross "b" dot "u" > 0.
c
c tol Input Numerical tolerance limit.
c On computers with 64-bit floating vector numbers,
c recommend 1.e-5 to 1.e-11.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
implicit none
integer iunit ! Unit for angles, 0 for degrees.
integer nerr ! Error indicator.
integer nerra ! Error indicator.
real anga ! Angle between vector "a" and vector "u".
real angb ! Angle between vector "b" and vector "u".
real arg ! Temporary term, square of a sine?
real ax ! Component x of vector "a".
real ay ! Component y of vector "a".
real az ! Component z of vector "a".
real bx ! Component x of vector "b".
real by ! Component y of vector "b".
real bz ! Component z of vector "b".
real cosab ! Cosine of angle between "a" and "b".
real cosa ! Cosine of angle between "a" and "u".
real cosb ! Cosine of angle between "b" and "u".
real ux ! Component x of vector "u".
real uy ! Component y of vector "u".
real uz ! Component z of vector "u".
real dab ! Length of vector ua cross ub.
real dlena ! Length of vector "a".
real dlenb ! Length of vector "b".
real dlenu ! Length of vector "u".
real erra ! Relative error in calculating fa.
real errb ! Relative error in calculating fb.
real fa ! Coefficient of vector ua.
real fab ! Coefficient of vector uab.
real fb ! Coefficient of vector ub.
real pi ! Math constant 3.14159265358979324...
real sinab2 ! Sine**2 of angle between "a" and "b",
real tol ! Precision constant.
real uabx ! Component x of vector uab.
real uaby ! Component y of vector uab.
real uabz ! Component z of vector uab.
real uax ! Component x of vector ua.
real uay ! Component y of vector ua.
real uaz ! Component z of vector ua.
real ubx ! Component x of vector ub.
real uby ! Component y of vector ub.
real ubz ! Component z of vector ub.
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptetru finding a vertex edge. tol = ',1pe22.14 /
cbug & ' ax,ay,az =',1p3e22.14 /
cbug & ' bx,by,bz =',1p3e22.14 /
cbug & ' anga,angb =',1p2e22.14 /
cbug & ' iunit =',i3 )
cbug write (3, 9901) tol, ax, ay, az, bx, by, bz,
cbug & anga, angb, iunit
cbugc***DEBUG ends.
c.... initialize.
pi = 3.14159265358979324
nerr = 0
ux = -1.e99
uy = -1.e99
uz = -1.e99
c.... Find the unit vector parallel to vector "a".
call aptvunb (ax, ay, az, 1, tol, uax, uay, uaz, dlena, nerra)
if (dlena .eq. 0.0) then
nerr = 1
go to 140
endif
c.... Find the unit vector parallel to vector "b".
call aptvunb (bx, by, bz, 1, tol, ubx, uby, ubz, dlenb, nerra)
if (dlenb .eq. 0.0) then
nerr = 2
go to 140
endif
c.... Find the cosine and sine**2 of the angle between side "a" and side "b"
c.... (vectors "a" and "b").
call aptvdot (uax, uay, uaz, ubx, uby, ubz, 1, tol, cosab, nerra)
if (abs (cosab - 1.0) .le. tol) then
nerr = 3
go to 140
endif
sinab2 = 1.0 - cosab**2
c.... Find the unit vector "uab" perpendicular to vectors "a" and "b".
call aptvxun (uax, uay, uaz, ubx, uby, ubz, 1, tol,
& uabx, uaby, uabz, dab, nerra)
c.... Find the cosines of the angles between vector "u" and vectors "a" and "b".
if (iunit .eq. 0) then
cosa = cos (anga * pi / 180.0)
cosb = cos (angb * pi / 180.0)
else
cosa = cos (anga)
cosb = cos (angb)
endif
c.... Find the unit vector "u" at angle anga from "a", angle angb from "b".
c.... Assume unit vector u = fa * ua + fb * ub + fab * uab,
c.... with unknowns fa, fb and fab. Then:
c.... u dot ua = cosa = fa + fb * cosab
c.... u dot ub = cosb = fa * cosab + fb
c.... u dot u = 1 = fa**2 + fb**2 + fab**2 + 2.0 * fa * fb * cosab
c.... Solving for fa, fb, fab, with truncation tests:
fa = (cosa - cosb * cosab) / sinab2
erra = tol * (abs (cosa) + abs (cosb * cosab)) / sinab2
if (abs (fa) .le. erra) fa = 0.0
fb = (cosb - cosa * cosab) / sinab2
errb = tol * (abs (cosb) + abs (cosa * cosab)) / sinab2
if (abs (fb) .le. errb) fb = 0.0
arg = fa**2 + fb**2 + 2.0 * fa * fb * cosab
if (abs (arg - 1.0) .le. tol) arg = 1.0
if (arg .gt. 1.0) then
nerr = 4
go to 140
endif
fab = sqrt (1.0 - arg)
call aptvdot (fa, fb, fab, uax, ubx, uabx, 1, tol, ux, nerra)
call aptvdot (fa, fb, fab, uay, uby, uaby, 1, tol, uy, nerra)
call aptvdot (fa, fb, fab, uaz, ubz, uabz, 1, tol, uz, nerra)
c.... Truncate small components to zero.
call aptvuna (ux, uy, uz, 1, tol, dlenu, nerra)
140 continue
cbugc***DEBUG begins.
cbug 9902 format (/ 'aptetru results. nerr = ',i3 /
cbug & ' ux,uy,uz =',1p3e22.14 )
cbug write ( 3, 9902) nerr, ux, uy, uz
cbugc***DEBUG ends.
210 return
c.... End of subroutine aptetru. (+1 line.)
end
UCRL-WEB-209832