subroutine aptripr (ab, bc, ca, tol, elen, da, db, dc, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTRIPR
c
c call aptripr (ab, bc, ca, tol, elen, da, db, dc, nerr)
c
c Version: aptripr Updated 2001 April 10 14:45.
c aptripr Originated 2001 April 2 15:30.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c Source: On the Compass Cluster EAST machine, in ~edwards/work/apt.
c
c Purpose: Given the edge lengths ab, bc and ca of triangle abc,
c to find the edge length elen of an equilateral triangle with
c vertices projected from vertices a, b and c, in the direction
c normal to the plane of triangle abc, and to find the distances
c da, db and dc of those vertices from vertices a, b and c.
c Flag nerr indicates any input errors.
c
c Input: ab, bc, ca, tol.
c
c Output: elen, da, db, dc, nerr.
c
c Calls: aptrich
c
c Glossary:
c
c ab Input The length of edge ab of triangle abc.
c
c bc Input The length of edge bc of triangle abc.
c
c ca Input The length of edge ca of triangle abc.
c
c da Output The distance from vertex A to the corresponding vertex
c of the equilateral triangle.
c
c db Output The distance from vertex B to the corresponding vertex
c of the equilateral triangle
c (arbitrarily set at zero).
c
c dc Output The distance from vertex C to the corresponding vertex
c of the equilateral triangle.
c
c elen Output The length of the edges of the equilateral triangle.
c
c nerr Output Indicates an input error, if not 0.
c 1 if an edge length is not greater than
c emin = tol * (abs (ab) + abs (bc) + abs (ca)).
c One or more edges are too short, or negative.
c 2 if an edge length is not less than
c emax = 0.5 * (ab + bc + ca) * (1.0 - tol).
c One edge is too long.
c
c tol Input Numerical tolerance limit.
c On computers with 64-bit floating point numbers,
c recommend 1.e-5 to 1.e-11.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Specifications.
implicit none
intrinsic max
intrinsic min
c.... Arguments.
integer nerr ! Error indicator.
real ab, bc, ca ! Length of edges ab, bc, ca of triangle.
real da, db, dc ! Distance vertices a, b, c projected.
real elen ! Length of edges of equilateral triangle.
real tol ! Numerical tolerence limit.
c.... Local variables.
real arg1, arg2 ! Intermediate.
real dab, dbc, dca ! Equal db - da, dc - db, da - cd.
real elen2 ! Square of elen.
real err2 ! Estimated error in arg2.
real edab, edbc, edca ! Estimated errors in dab, dbc, dca.
real shift ! Equals min (da, db, dc).
real sum1 ! Equals dab + dbc - dca.
real sum2 ! Equals dab - dbc + dca.
real sum3 ! Equals dab - dbc - dca.
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptripr finding equilateral projection triangle.',
cbug & ' tol=',1pe22.14 /
cbug & ' ab,bc,ca = ',1p3e22.14 )
cbug write ( 3, 9901) tol, ab, bc, ca
cbugc***DEBUG ends.
c.... Initialize.
nerr = 0
elen = -1.e99
da = -1.e99
db = -1.e99
dc = -1.e99
c.... Test for impossible triangle (an edge too short or long).
call aptrich (ab, bc, ca, tol, nerr)
if (nerr .ne. 0) go to 410
c.... Find the equilateral triangle edge length elen.
arg1 = ab**2 + bc**2 + ca**2
arg2 = ab**4 + bc**4 + ca**4 -
& ab**2 * bc**2 - bc**2 * ca**2 - ca**2 * ab**2
err2 = tol * (ab**4 + bc**4 + ca**4 +
& ab**2 * bc**2 + bc**2 * ca**2 + ca**2 * ab**2)
if (abs (arg2) .le. err2) arg2 = 0.0
elen2 = (arg1 + 2.0 * sqrt (arg2)) / 3.0
elen = sqrt (elen2)
c.... Find the distances between the original and projected vertices.
dab = elen2 - ab**2
edab = tol * (elen2 + ab**2)
if (abs (dab) .le. edab) dab = 0.0
dab = sqrt (dab)
dbc = elen2 - bc**2
edbc = tol * (elen2 + bc**2)
if (abs (dbc) .le. edbc) dbc = 0.0
dbc = sqrt (dbc)
dca = elen2 - ca**2
edca = tol * (elen2 + ca**2)
if (abs (dca) .le. edca) dca = 0.0
dca = sqrt (dca)
sum1 = abs (dab + dbc - dca)
sum2 = abs (dab - dbc + dca)
sum3 = abs (dab - dbc - dca)
cbugc***DEBUG begins.
cbug 9914 format (' sum1,2,3 = '1p3e22.14 /
cbug & ' dab,bc,ca= '1p3e22.14 )
cbug write ( 3, 9914) sum1, sum2, sum3, dab, dbc, dca
cbugc***DEBUG ends.
c.... Find the correct signs, so that dca + dab + dbc = 0.
if ((sum1 .le. sum2) .and. (sum1 .le. sum3)) then
dca = -dca
elseif ((sum2 .le. sum3) .and. (sum2 .le. sum1)) then
dbc = -dbc
elseif ((sum3 .le. sum1) .and. (sum3 .le. sum2)) then
dbc = -dbc
dca = -dca
endif
cbugc***DEBUG begins.
cbug write ( 3, 9914) sum1, sum2, sum3, dab, dbc, dca
cbugc***DEBUG ends.
da = 0.0
db = dab
dc = -dca
c.... Make all distances non-negative, with one zero.
shift = min (da, db, dc)
da = da - shift
db = db - shift
dc = dc - shift
410 continue
cbugc***DEBUG begins.
cbug 9916 format (/ 'aptripr results. nerr = ',i2 /
cbug & ' elen = '1pe22.14 /
cbug & ' da,db,dc = '1p3e22.14 )
cbug write ( 3, 9916) nerr, elen, da, db, dc
cbugc***DEBUG ends.
return
c.... End of subroutine aptripr. (+1 line.)
end
UCRL-WEB-209832