subroutine aptsver (ax, ay, az, bx, by, bz, cx, cy, cz, tol,
& sa, sb, sc, sum, wa, wb, wc, sx, sy, sz, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTSVER
c
c call aptsver (ax, ay, az, bx, by, bz, cx, cy, cz, tol,
c & sa, sb, sc, sum, wa, wb, wc, sx, sy, sz, nerr)
c
c Version: aptsver Updated 2001 March 13 13:45.
c aptsver Originated 2001 March 7 13:40.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: For the triangle with vertices a = (ax, ay, az),
c b = (bx, by, bz), and c = (cx, cy, cz) in 3-D space, to find
c the Steiner vertex s = (sx, sy, sz), from which the sum of the
c lengths sum = sa + sb + sc of the lines to the three triangle
c vertices is minimized. The angles between each pair of those
c lines is 120 degrees, unless the triangle has a vertex angle
c exceeding 120 degrees, in which case the Steiner vertex is at
c that vertex. The relative vertex weights wa, wb and wc are
c also returned.
c Flag nerr indicates any failure to find a true Steiner vertex.
c
c Input: ax, ay, az, bx, by, bz, cx, cy, cz, tol.
c
c Output: sa, sb, sc, wa, wb, wc, sx, sy, sz, nerr.
c
c Calls: aptvang, aptvdis
c
c Glossary:
c
c ax,ay,az Input The x, y, z coordinates of vertex "a" of the triangle.
c
c bx,by,bz Input The x, y, z coordinates of vertex "b" of the triangle.
c
c cx,cy,cz Input The x, y, z coordinates of vertex "c" of the triangle.
c
c ab Output The length of edge "ab" of the triangle.
c
c bc Output The length of edge "bc" of the triangle.
c
c ca Output The length of edge "ca" of the triangle.
c
c sa,b,c Output The distances from the Steiner vertex "s" to the
c triangle vertices "a", "b" and "c".
c
c sum Output The sum of distances, sum = sa + sb + sc.
c
c nerr Output Indicates special case, if not 0.
c 1 if the triangle has two coincident vertices.
c Steiner vertex is put there.
c 2 if the triangle has a vertex angle of 120 degrees
c or more. Steiner vertex is put there.
c
c sx,sy,sz Output The x, y, z coordinates of the Steiner vertex.
c
c tol Input Numerical tolerance limit.
c On computers with 64-bit real numbers, recommend
c 1.e-5 to 1.e-11.
c
c wa,b,c Output Relative eights of triangle vertices "a", "b" and "c",
c relative to Steiner vertex.
c sx = wa * ax + wb * bx + wc * cx
c sy = wa * ay + wb * by + wc * cy
c sz = wa * az + wb * bz + wc * cz.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
c.... Specify variable types.
implicit none
c.... Arguments.
integer nerr ! Special case indicator.
real ax ! Coordinate x of vertex a.
real ay ! Coordinate y of vertex a.
real az ! Coordinate z of vertex a.
real bx ! Coordinate x of vertex b.
real by ! Coordinate y of vertex b.
real bz ! Coordinate z of vertex b.
real cx ! Coordinate x of vertex c.
real cy ! Coordinate y of vertex c.
real cz ! Coordinate z of vertex c.
real sa ! Distance from Steiner vertex to point a.
real sb ! Distance from Steiner vertex to point b.
real sc ! Distance from Steiner vertex to point c.
real sum ! Sub sa + sb + sc.
real sx ! Coordinate x of Steiner vertex.
real sy ! Coordinate y of Steiner vertex.
real sz ! Coordinate z of Steiner vertex.
real tol ! Precision limit.
real wa ! Weight of vertex a.
real wb ! Weight of vertex b.
real wc ! Weight of vertex c.
c.... Local variables.
integer nerra ! Error return from apt subroutine.
real ab ! Length of edge vector ab.
real abx ! Component x of edge vector ab.
real aby ! Component y of edge vector ab.
real abz ! Component z of edge vector ab.
real anga ! Angle at vertex a, radians.
real angb ! Angle at vertex b, radians.
real angc ! Angle at vertex c, radians.
real angm ! Angle of 120 degrees, in radians.
real angra ! Angle angm - anga.
real angrb ! Angle angm - angb.
real angrc ! Angle angm - angc.
real bc ! Length of edge vector bc.
real bcx ! Component x of edge vector bc.
real bcy ! Component y of edge vector bc.
real bcz ! Component z of edge vector bc.
real ca ! Length of edge vector ca.
real cax ! Component x of edge vector ca.
real cay ! Component y of edge vector ca.
real caz ! Component z of edge vector ca.
real costha ! Cosine of angle at vertex a.
real costhb ! Cosine of angle at vertex b.
real costhc ! Cosine of angle at vertex c.
real den ! Denominator sa*sb + sb*sc + sc*sa.
real dena ! Denominator for finding sa.
real denb ! Denominator for finding sb.
real denc ! Denominator for finding sc.
real fact ! Constant fact = 2.0 / sqrt (3.0).
real pi ! Constant pi = 3.14159265358979323...
c.... Intrinsic functions.
intrinsic cos ! Cosine function.
intrinsic sin ! Sine function.
c=======================================================================********
cbugc***DEBUG begins.
cbug 9901 format (/ 'aptsver finding triangle Steiner vertex.',
cbug & ' tol=',1pe22.14 /
cbug & ' ax,ay,az= ',1p3e22.14 /
cbug & ' bx,by,bz= ',1p3e22.14 /
cbug & ' cx,cy,cz= ',1p3e22.14 )
cbug write ( 3, 9901) tol, ax, ay, az, bx, by, bz, cx, cy, cz
cbugc***DEBUG ends.
c=======================================================================********
c.... Initialize.
nerr = 0
c=======================================================================********
c.... Find the edge vectors of the triangle, "ab", "bc" and "ca",
c.... and the edge lengths, ab, bc and ca.
c.... Put the Steiner vertex at any edge of zero length.
call aptvdis (ax, ay, az, bx, by, bz, 1, tol,
& abx, aby, abz, ab, nerra)
call aptvdis (bx, by, bz, cx, cy, cz, 1, tol,
& bcx, bcy, bcz, bc, nerra)
call aptvdis (cx, cy, cz, ax, ay, az, 1, tol,
& cax, cay, caz, ca, nerra)
if (ab .eq. 0.0) then ! Steiner vertex at vertex A = B.
nerr = 1
sa = 0.0
sb = 0.0
sc = ca
wa = 0.5
wb = 0.5
wc = 0.0
go to 150
endif
if (bc .eq. 0.0) then ! Steiner vertex at vertex B = C.
nerr = 1
sa = ab
sb = 0.0
sc = 0.0
wa = 0.0
wb = 0.5
wc = 0.5
go to 150
endif
if (ca .eq. 0.0) then ! Steiner vertex at vertex C = A.
nerr = 1
sa = 0.0
sb = bc
sc = 0.0
wa = 0.5
wb = 0.0
wc = 0.5
go to 150
endif
c=======================================================================********
c.... Find the vertex angles.
c.... Put the Steiner vertex at any triangle vertex of 120 degrees or more.
pi = 3.14159265358979323
angm = 2.0 * pi / 3.0 ! 120 degrees.
call aptvang (cax, cay, caz, abx, aby, abz, 1, tol, costha, nerra)
anga = acos (-costha)
if (abs (anga - angm) .le. tol * angm) anga = angm
if (anga .ge. angm) then ! Steiner vertex at vertex A.
nerr = 2
sa = 0.0
sb = ab
sc = ca
wa = 1.0
wb = 0.0
wc = 0.0
go to 150
endif
call aptvang (abx, aby, abz, bcx, bcy, bcz, 1, tol, costhb, nerra)
angb = acos (-costhb)
if (abs (angb - angm) .le. tol * angm) angb = angm
if (angb .ge. angm) then ! Steiner vertex at vertex B.
nerr = 2
sa = ab
sb = 0.0
sc = bc
wa = 0.0
wb = 1.0
wc = 0.0
go to 150
endif
call aptvang (bcx, bcy, bcz, cax, cay, caz, 1, tol, costhc, nerra)
angc = acos (-costhc)
if (abs (angc - angm) .le. tol * angm) angc = angm
if (angc .ge. angm) then ! Steiner vertex at vertex C.
nerr = 2
sa = ca
sb = bc
sc = 0.0
wa = 0.0
wb = 0.0
wc = 1.0
go to 150
endif
c=======================================================================********
c.... Find the distances from the Steiner vertex to the triangle vertices.
fact = 2.0 / sqrt (3.0)
angra = angm - anga
angrb = angm - angb
angrc = angm - angc
dena = sqrt (ca**2 + ab**2 + 2.0 * ca * ab * cos (angra))
denb = sqrt (ab**2 + bc**2 + 2.0 * ab * bc * cos (angrb))
denc = sqrt (bc**2 + ca**2 + 2.0 * bc * ca * cos (angrc))
sa = fact * ca * ab * sin (angra) / dena
sb = fact * ab * bc * sin (angrb) / denb
sc = fact * bc * ca * sin (angrc) / denc
den = sa * sb + sb * sc + sc * sa
wa = sb * sc / den
wb = sc * sa / den
wc = sa * sb / den
c=======================================================================********
c.... Find the total length and the Steiner vertex.
150 sum = sa + sb + sc
sx = wa * ax + wb * bx + wc * cx
sy = wa * ay + wb * by + wc * cy
sz = wa * az + wb * bz + wc * cz
c=======================================================================********
210 continue
cbugc***DEBUG begins.
cbug 9916 format (/ 'aptsver results. nerr = ',i2 /
cbug & ' ab,bc,ca ',1p3e22.14 /
cbug & '