subroutine apttrig (ax, ay, az, bx, by, bz, cx, cy, cz, tol,
& abx, aby, abz, dab, bcx, bcy, bcz, dbc,
& cax, cay, caz, dca, anga, angb, angc,
& vnx, vny, vnz, areatr,
& vmax, vmay, vmaz, dma, pmax, pmay, pmaz,
& vmbx, vmby, vmbz, dmb, pmbx, pmby, pmbz,
& vmcx, vmcy, vmcz, dmc, pmcx, pmcy, pmcz,
& cgx, cgy, cgz,
& viax, viay, viaz, dvia, anax, anay, anaz,
& vibx, viby, vibz, dvib, anbx, anby, anbz,
& vicx, vicy, vicz, dvic, ancx, ancy, ancz,
& cix, ciy, ciz, radinsc,
& vqabx, vqaby, vqabz, dqab,
& vqbcx, vqbcy, vqbcz, dqbc,
& vqcax, vqcay, vqcaz, dqca,
& ccx, ccy, ccz, radcirc,
& vhax, vhay, vhaz, dha, alax, alay, alaz,
& vhbx, vhby, vhbz, dhb, albx, alby, albz,
& vhcx, vhcy, vhcz, dhc, alcx, alcy, alcz,
& cox, coy, coz, nerr)
ccbeg.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c SUBROUTINE APTTRIG
c
c call apttrig (ax, ay, az, bx, by, bz, cx, cy, cz, tol,
c & abx, aby, abz, dab, bcx, bcy, bcz, dbc,
c & cax, cay, caz, dca, anga, angb, angc,
c & vnx, vny, vnz, areatr,
c & vmax, vmay, vmaz, dma, pmax, pmay, pmaz,
c & vmbx, vmby, vmbz, dmb, pmbx, pmby, pmbz,
c & vmcx, vmcy, vmcz, dmc, pmcx, pmcy, pmcz,
c & cgx, cgy, cgz,
c & viax, viay, viaz, dvia, anax, anay, anaz,
c & vibx, viby, vibz, dvib, anbx, anby, anbz,
c & vicx, vicy, vicz, dvic, ancx, ancy, ancz,
c & cix, ciy, ciz, radinsc,
c & vqabx, vqaby, vqabz, dqab,
c & vqbcx, vqbcy, vqbcz, dqbc,
c & vqcax, vqcay, vqcaz, dqca,
c & ccx, ccy, ccz, radcirc,
c & vhax, vhay, vhaz, dha, alax, alay, alaz,
c & vhbx, vhby, vhbz, dhb, albx, alby, albz,
c & vhcx, vhcy, vhcz, dhc, alcx, alcy, alcz,
c & cox, coy, coz, nerr)
c
c Version: apttrig Updated 1994 March 24 14:00.
c apttrig Originated 1994 March 24 14:00.
c
c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123.
c
c
c Purpose: To find geometric data for the triangle with vertices
c a = (ax, ay, az), b = (bx, by, bz), and c = (cx, cy, cz) in
c 3-D space. The data includes the edge vectors and their
c lengths, the vertex angles, the vector normal to the plane of
c the triangle, the area of the triangle, the median vectors,
c their lengths, their intercepts (the midpoints of the edges),
c the centroid or center of gravity "cg", the angle bisector
c vectors, their lengths, their intercepts, the center of the
c inscribed circle "ci" and its radius, the perpendicular
c bisector vectors, the center of the circumscribed circle "cc"
c and its radius, the altitude vectors, their lengths, their
c intercepts, and the orthocenter "co".
c Flag nerr indicates any input error.
c
c Input: 10 arguments:
c ax, ay, az, bx, by, bz, cx, cy, cz, tol Vertices, tolerance
c
c Output: 109 arguments:
c GENERAL (19):
c abx, aby, abz, dab Vector "ab", length
c bcx, bcy, bcz, dbc Vector "bc", length
c cax, cay, caz, dca Vector "ca", length
c anga, angb, angc Vertex angles (deg)
c vnx, vny, vnz, areatr Normal vector, area
c RELATED TO THE MEDIANS (24):
c vmax, vmay, vmaz, dma, pmax, pmay, pmaz Median at "a", ...
c vmbx, vmby, vmbz, dmb, pmbx, pmby, pmbz Median at "b", ...
c vmcx, vmcy, vmcz, dmc, pmcx, pmcy, pmcz Median at "c", ...
c cgx, cgy, cgz Centroid
c RELATED TO THE ANGLE BISECTORS (25):
c viax, viay, viaz, dvia, anax, anay, anaz Angle bisector at "a"
c vibx, viby, vibz, dvib, anbx, anby, anbz Angle bisector at "b"
c vicx, vicy, vicz, dvic, ancx, ancy, ancz Angle bisector at "c"
c cix, ciy, ciz, Center inscr circle
c radinsc Radius inscr circle
c RELATED TO THE PERPENDICULAR BISECTORS (16):
c vqabx, vqaby, vqabz, dqab Perp bisector of "ab"
c vqbcx, vqbcy, vqbcz, dqbc Perp bisector of "bc"
c vqcax, vqcay, vqcaz, dqca Perp bisector of "ca"
c ccx, ccy, ccz Center circum circle
c radcirc Radius circum circle
c RELATED TO THE ALTITUDES (24):
c vhax, vhay, vhaz, dha, alax, alay, alaz Altitude at "a", ...
c vhbx, vhby, vhbz, dhb, albx, alby, albz Altitude at "b", ...
c vhcx, vhcy, vhcz, dhc, alcx, alcy, alcz Altitude at "c", ...
c cox, coy, coz Orthocenter
c ERROR FLAG (1):
c nerr 1 if zero area
c
c Calls: aptcirp, aptlnln, apttrid, aptvang, aptvaxb, aptvdis, aptvsum,
c aptvdot
c
c
c Glossary:
c
c abx,y,z Output The x,y,z components of the edge vector "ab" from
c vertex "a" to vertex "b".
c
c alax,y,z Output The x,y,z coordinates of the intercept of the altitude
c vector "vha" from vertex "a" on edge "bc".
c
c albx,y,z Output The x,y,z coordinates of the intercept of the altitude
c vector "vhb" from vertex "b" on edge "ca".
c
c alcx,y,z Output The x,y,z coordinates of the intercept of the altitude
c vector "vhc" from vertex "c" on edge "ab".
c
c anax,y,z Output The x,y,z coordinates of the intercept of the angle
c bisector vector "via" from vertex "a" on edge "bc".
c
c anbx,y,z Output The x,y,z coordinates of the intercept of the angle
c bisector vector "vib" from vertex "b" on edge "ca".
c
c ancx,y,z Output The x,y,z coordinates of the intercept of the angle
c bisector vector "vic" from vertex "c" on edge "ab".
c
c anga Output The interior angle at vertex "a", degrees.
c Set to -1.E99 if the triangle has zero area.
c
c angb Output The interior angle at vertex "b", degrees.
c Set to -1.E99 if the triangle has zero area.
c
c angb Output The interior angle at vertex "c", degrees.
c Set to -1.E99 if the triangle has zero area.
c
c bcx,y,z Output The x,y,z components of the edge vector "bc" from
c vertex "b" to vertex "c".
c
c cax,y,z Output The x,y,z components of the edge vector "ca" from
c vertex "c" to vertex "a".
c
c cox,y,z Output The x,y,z coordinates of the orthocenter (intersection
c of altitudes).
c
c areatr Output Area of the triangle.
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 ccx,y,z Output The x, y, z coordinates of the center of the
c circumscribed circle (at the intersection of the
c perpendicular bisectors of the edges, "vqab", "vqbc"
c and "vqca"). The radius is radcirc.
c Very large if the triangle has zero area (nerr = 1).
c
c cgx,y,z Output The x, y, z coordinates of the center of gravity or
c centroid of the triangle (at the intersection of the
c medians "vma", "vmb" and "vmc").
c
c cix,y,z Output The x, y, z coordinates of the center of the
c inscribed circle (at the intersection of the vertex
c angle bisectors "via", "vib" and "vic").
c The radius is radinsc.
c
c cx,cy,cz Input The x, y, z coordinates of vertex "c" of the triangle.
c
c dab Output The length of edge "ab" of the triangle.
c
c dbc Output The length of edge "bc" of the triangle.
c
c dca Output The length of edge "ca" of the triangle.
c
c dha,b,c Output The lengths of the altitude vectors "vha", "vhb", and
c "vhc", from the vertices "a", "b" and "c" to the
c intercepts "ala", "alb" and "alc" on the opposite
c edges.
c
c dma,b,c Output The lengths of the median vectors "vma", "vmb", "vmc",
c from vertices "a", "b" and "c" to the midpoints
c "pma", "pmb" and "pmc" of the opposite edges.
c
c dqab Output The distance from the midpoint "pmc" of edge "ab" to
c the center of the circumscribed circle, "cc".
c
c dqbc Output The distance from the midpoint "pma" of edge "bc" to
c the center of the circumscribed circle, "cc".
c
c dqca Output The distance from the midpoint "pmb" of edge "ca" to
c the center of the circumscribed circle, "cc".
c
c dvia,b,c Output The lengths of the angle bisector vectors "via", "vib"
c and "vic", from vertices "a", "b" and "c" to the
c opposite edges.
c
c nerr Output Indicates an input error, if not 0.
c 1 if the triangle has zero area (vertices are
c coincident or colinear).
c
c pmax,y,z Output The x,y,z coordinates of the midpoint of edge "bc".
c The intercept of median "vma" on edge "bc".
c
c pmbx,y,z Output The x,y,z coordinates of the midpoint of edge "ca".
c The intercept of median "vmb" on edge "ca".
c
c pmcx,y,z Output The x,y,z coordinates of the midpoint of edge "ab".
c The intercept of median "vmc" on edge "ab".
c
c radcirc Output The radius of the circumscribed circle. Centered at
c "cc". Very large if the triangle has zero area.
c
c radinsc Output The radius of the inscribed circle. Centered at "ci".
c
c tol Input Numerical tolerance limit.
c On Cray computers, recommend 1.e-5 to 1.e-11.
c
c vhax,y,z Output The x,y,z components of the altitude vector from vertex
c "a" to edge "bc". Length is dha.
c Intersects the other altitude vectors at the
c orthocenter "co".
c
c vhbx,y,z Output The x,y,z components of the altitude vector from vertex
c "b" to edge "ca". Length is dhb.
c Intersects the other altitude vectors at the
c orthocenter "co".
c
c vhcx,y,z Output The x,y,z components of the altitude vector from vertex
c "c" to edge "ab". Length is dhc.
c Intersects the other altitude vectors at the
c orthocenter "co".
c
c viax,y,z Output The x,y,z components of the angle bisector vector from
c vertex "a" to edge "bc". Length is dvia.
c Intersects the other angle bisector vectors at the
c center of the inscribed circle "ci".
c
c vibx,y,z Output The x,y,z components of the angle bisector vector from
c vertex "b" to edge "ca". Length is dvib.
c Intersects the other angle bisector vectors at the
c center of the inscribed circle "ci".
c
c vicx,y,z Output The x,y,z components of the angle bisector vector from
c vertex "c" to edge "ab". Length is dvic.
c Intersects the other angle bisector vectors at the
c center of the inscribed circle "ci".
c
c vmax,y,z Output The x,y,z components of the median vector "vma" from
c vertex "a" to the center of edge "bc". Intersects
c the other median vectors at the centroid, "cg".
c Length is dma.
c
c vmbx,y,z Output The x,y,z components of the median vector "vmb" from
c vertex "b" to the center of edge "ca". Intersects
c the other median vectors at the centroid, "cg".
c Length is dmb.
c
c vmcx,y,z Output The x,y,z components of the median vector "vmc" from
c vertex "c" to the center of edge "ab". Intersects
c the other median vectors at the centroid, "cg".
c Length is dmc.
c
c vnx,y,z Output The x,y,z components of the vector normal to the plane
c of the triangle. Length is 2.0 * areatr.
c
c vqabx,y,z Output The x,y,z components of the perpendicular bisector
c vector from the midpoint of edge "ab" to the center
c of the circumscribed circle, "cc". Length is dqab.
c
c vqbcx,y,z Output The x,y,z components of the perpendicular bisector
c vector from the midpoint of edge "bc" to the center
c of the circumscribed circle, "cc". Length is dqbc.
c
c vqcax,y,z Output The x,y,z components of the perpendicular bisector
c vector from the midpoint of edge "ca" to the center
c of the circumscribed circle, "cc". Length is dqca.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
ccend.
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c Local glossary.
c
c dota Local The dot product of the edge vectors "ca" and "ab".
c
c dotb Local The dot product of the edge vectors "ab" and "bc".
c
c dotc Local The dot product of the edge vectors "bc" and "ca".
c
c vpabx,y,z Local The x,y,z components of a vector perpendicular to
c edge "ab", in the plane of the triangle.
c
c vpbcx,y,z Local The x,y,z components of a vector perpendicular to
c edge "bc", in the plane of the triangle.
c
c vpcax,y,z Local The x,y,z components of a vector perpendicular to
c edge "ca", in the plane of the triangle.
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c.... Parameters.
parameter (small = 1.e-99)
parameter (big = 1.e+99)
parameter (pi = 3.14159265358979323)
c.... Local variables.
cbugc***DEBUG begins.
cbug 9901 format (/ 'apttrig finding triangle data:' /
cbug & ' ax,ay,az= ',1p3e22.14 /
cbug & ' bx,by,bz= ',1p3e22.14 /
cbug & ' cx,cy,cz= ',1p3e22.14 )
cbug write ( 3, 9901) ax, ay, az, bx, by, bz, cx, cy, cz
cbugc***DEBUG ends.
c.... Initialize.
nerr = 0
c=======================================================================********
c.... Find the edge vectors of the triangle, "ab", "bc" and "ca",
c.... and the edge lengths, dab, dbc and dca.
call aptvdis (ax, ay, az, bx, by, bz, 1, tol,
& abx, aby, abz, dab, nerra)
call aptvdis (bx, by, bz, cx, cy, cz, 1, tol,
& bcx, bcy, bcz, dbc, nerra)
call aptvdis (cx, cy, cz, ax, ay, az, 1, tol,
& cax, cay, caz, dca, nerra)
if ((dab .eq. 0.0) .or.
& (dbc .eq. 0.0) .or.
& (dca .eq. 0.0)) then
nerr = 1
endif
cbugc***DEBUG begins.
cbug 9902 format (/
cbug & 'apttrig results (edge vectors and lengths):' /
cbug & ' abx,aby,abz= ',1p3e22.14 /
cbug & ' bcx,bcy,bcz= ',1p3e22.14 /
cbug & ' cax,cay,caz= ',1p3e22.14 /
cbug & ' dab,dbc,dca= ',1p3e22.14 )
cbug write ( 3, 9902) abx, aby, abz, bcx, bcy, bcz, cax, cay, caz,
cbug & dab, dbc, dca
cbugc***DEBUG ends.
c=======================================================================********
c.... Find the vertex angles.
call aptvang (cax, cay, caz, abx, aby, abz, 1, tol, costha, nerra)
call aptvang (abx, aby, abz, bcx, bcy, bcz, 1, tol, costhb, nerra)
call aptvang (bcx, bcy, bcz, cax, cay, caz, 1, tol, costhc, nerra)
anga = (180.0 / pi) * acos (-costha)
angb = (180.0 / pi) * acos (-costhb)
angc = (180.0 / pi) * acos (-costhc)
if (dab .eq. 0.0) then
anga = -big
angb = -big
endif
if (dbc .eq. 0.0) then
angb = -big
angc = -big
endif
if (dca .eq. 0.0) then
angc = -big
anga = -big
endif
cbugc***DEBUG begins.
cbug 9903 format (/
cbug & 'apttrig results (vertex angles, degrees):' /
cbug & ' anga,b,c= ',1p3e22.14 )
cbug write ( 3, 9903) anga, angb, angc
cbugc***DEBUG ends.
c=======================================================================********
c.... Find the normal vector, and the area of the triangle.
call aptvaxb (abx, aby, abz, bcx, bcy, bcz, 1, tol,
& vnx, vny, vnz, dnorm, nerra)
if (dnorm .eq. 0.0) then
nerr = 1
endif
areatr = 0.5 * dnorm
cbugc***DEBUG begins.
cbug 9904 format (/
cbug & 'apttrig results (normal vector and length, triangle area):' /
cbug & ' vnx,vny,vnz= ',1p3e22.14 /
cbug & ' dnorm= ',1pe22.14 /
cbug & ' area= ',1pe22.14 )
cbug write ( 3, 9904) vnx, vny, vnz, dnorm, areatr
cbugc***DEBUG ends.
c=======================================================================********
c.... Find the area of the triangle, and
c.... the coordinates of the center of gravity.
call apttrid (ax, ay, az, bx, by, bz, cx, cy, cz, 1, tol,
& dab, dbc, dca, areatr, cgx, cgy, cgz, nerra)
if (areatr .eq. 0.0) then
nerr = 1
endif
cbugc***DEBUG begins.
cbug 9905 format (/
cbug & 'apttrig results (edge lengths, area, center of gravity):' /
cbug & ' dab,dbc,dca= ',1p3e22.14 /
cbug & ' area= ',1pe22.14 /
cbug & ' cgx,cgy,cgz= ',1p3e22.14 )
cbug write ( 3, 9905) dab, dbc, dca, areatr, cgx, cgy, cgz
cbugc***DEBUG ends.
c=======================================================================********
c.... Find the vectors perpendicular to the edges, in the plane of the triangle.
call aptvaxb (vnx, vny, vnz, abx, aby, abz, 1, tol,
& vpabx, vpaby, vpabz, vlenab, nerra)
call aptvaxb (vnx, vny, vnz, bcx, bcy, bcz, 1, tol,
& vpbcx, vpbcy, vpbcz, vlenbc, nerra)
call aptvaxb (vnx, vny, vnz, cax, cay, caz, 1, tol,
& vpcax, vpcay, vpcaz, vlenca, nerra)
if ((vlenab .eq. 0.0) .or.
& (vlenbc .eq. 0.0) .or.
& (vlenca .eq. 0.0)) then
nerr = 1
endif
cbugc***DEBUG begins.
cbug 9906 format (/
cbug & 'apttrig results (vectors perpendicular to edges, lengths):' /
cbug & ' vpabx,y,z= ',1p3e22.14 /
cbug & ' vpbcx,y,z= ',1p3e22.14 /
cbug & ' vpcax,y,z= ',1p3e22.14 /
cbug & ' lengths= ',1p3e22.14 )
cbug write ( 3, 9906) vpabx, vpaby, vpabz,
cbug & vpbcx, vpbcy, vpbcz,
cbug & vpcax, vpcay, vpcaz, vlenab, vlenbc, vlenca
cbugc***DEBUG ends.
c=======================================================================********
c.... Find the midpoints of the edges of the triangle.
call aptvsum (0, 0.5, ax, ay, az, 0.5, bx, by, bz, 1, tol,
& pmcx, pmcy, pmcz, dpmc, nerra)
call aptvsum (0, 0.5, bx, by, bz, 0.5, cx, cy, cz, 1, tol,
& pmax, pmay, pmaz, dpma, nerra)
call aptvsum (0, 0.5, cx, cy, cz, 0.5, ax, ay, az, 1, tol,
& pmbx, pmby, pmbz, dpmb, nerra)
cbugc***DEBUG begins.
cbug 9907 format (/
cbug & 'apttrig results (midpoints of edges):' /
cbug & ' pmax,y,z= ',1p3e22.14 /
cbug & ' pmbx,y,z= ',1p3e22.14 /
cbug & ' pmcx,y,z= ',1p3e22.14 )
cbug write ( 3, 9907) pmax, pmay, pmaz, pmbx, pmby, pmbz,
cbug & pmcx, pmcy, pmcz
cbugc***DEBUG ends.
c=======================================================================********
c.... Find the vectors from the vertices to the median of the opposite edge.
c.... These meet at the centroid, or center of gravity, of the triangle.
call aptvsum (0, 0.5, abx, aby, abz, -0.5, cax, cay, caz, 1, tol,
& vmax, vmay, vmaz, dma, nerra)
call aptvsum (0, 0.5, bcx, bcy, bcz, -0.5, abx, aby, abz, 1, tol,
& vmbx, vmby, vmbz, dmb, nerra)
call aptvsum (0, 0.5, cax, cay, caz, -0.5, bcx, bcy, bcz, 1, tol,
& vmcx, vmcy, vmcz, dmc, nerra)
cbugc***DEBUG begins.
cbug 9908 format (/
cbug a 'apttrig results (median vectors, lengths):' /
cbug & ' vmax,y,z= ',1p3e22.14 /
cbug & ' vmbx,y,z= ',1p3e22.14 /
cbug & ' vmcx,y,z= ',1p3e22.14 /
cbug & ' dma,b,c ',1p3e22.14 )
cbug write ( 3, 9908) vmax, vmay, vmaz,
cbug & vmbx, vmby, vmbz,
cbug & vmcx, vmcy, vmcz, dma, dmb, dmc
cbugc***DEBUG ends.
c=======================================================================********
c.... Find the altitude vectors, and their heights.
call aptvaxb (vnx, vny, vnz, bcx, bcy, bcz, 1, tol,
& vhax, vhay, vhaz, vlen, nerra)
vhax = vpbcx / (dbc**2 + small)
vhay = vpbcy / (dbc**2 + small)
vhaz = vpbcz / (dbc**2 + small)
dha = 2.0 * areatr / (dbc + small)
call aptvaxb (vnx, vny, vnz, cax, cay, caz, 1, tol,
& vhbx, vhby, vhbz, vlen, nerra)
vhbx = vpcax / (dca**2 + small)
vhby = vpcay / (dca**2 + small)
vhbz = vpcaz / (dca**2 + small)
dhb = 2.0 * areatr / (dca + small)
call aptvaxb (vnx, vny, vnz, abx, aby, abz, 1, tol,
& vhcx, vhcy, vhcz, vlen, nerra)
vhcx = vpabx / (dab**2 + small)
vhcy = vpaby / (dab**2 + small)
vhcz = vpabz / (dab**2 + small)
dhc = 2.0 * areatr / (dab + small)
cbugc***DEBUG begins.
cbug 9909 format (/
cbug & 'apttrig results (altitude vectors and their heights):' /
cbug & ' vhax,y,z= ',1p3e22.14 /
cbug & ' vhbx,y,z= ',1p3e22.14 /
cbug & ' vhcx,y,z= ',1p3e22.14 /
cbug & ' heights= ',1p3e22.14 )
cbug write ( 3, 9909) vhax, vhay, vhaz,
cbug & vhbx, vhby, vhbz,
cbug & vhcx, vhcy, vhcz,
cbug & dha, dhb, dhc
cbugc***DEBUG ends.
c=======================================================================********
c.... Find the intercepts of the altitudes on the edges.
call aptvsum (0, 1.0, ax, ay, az, -1.0, vhax, vhay, vhaz,
& 1, tol, alax, alay, alaz, dala, nerra)
call aptvsum (0, 1.0, bx, by, bz, -1.0, vhbx, vhby, vhbz,
& 1, tol, albx, alby, albz, dalb, nerra)
call aptvsum (0, 1.0, cx, cy, cz, -1.0, vhcx, vhcy, vhcz,
& 1, tol, alcx, alcy, alcz, dalc, nerra)
cbugc***DEBUG begins.
cbug 9910 format (/
cbug & 'apttrig results (altitude intercepts):' /
cbug & ' alax,y,z= ',1p3e22.14 /
cbug & ' albx,y,z= ',1p3e22.14 /
cbug & ' alcx,y,z= ',1p3e22.14 )
cbug write ( 3, 9910) alax, alay, alaz,
cbug & albx, alby, albz,
cbug & alcx, alcy, alcz
cbugc***DEBUG ends.
c=======================================================================********
c.... Find the orthocenter (triple intersections of altitudes).
c.... Intersect line "a" to "ala" and "b" to "alb" and "c" to "alc".
call aptlnln (ax, ay, az, alax, alay, alaz,
& bx, by, bz, albx, alby, albz, 1, tol,
& dpmin, fracaa, fracbb,
& cox, coy, coz, cxx, cyy, czz, itrun, ipar, nerra)
cbugc***DEBUG begins.
cbug 9920 format (/
cbug & 'apttrig results (altitude intersection = orthocenter):' /
cbug & ' cox,coy,coz= ',1p3e22.14 )
cbug write ( 3, 9920) cox, coy, coz
cbugc***DEBUG ends.
c=======================================================================********
c.... Find the dot products of the edge vectors.
call aptvdot (cax, cay, caz, abx, aby, abz, 1, tol, dota, nerra)
call aptvdot (abx, aby, abz, bcx, bcy, bcz, 1, tol, dotb, nerra)
call aptvdot (bcx, bcy, bcz, cax, cay, caz, 1, tol, dotc, nerra)
cbugc***DEBUG begins.
cbug 9911 format (/
cbug & 'apttrig results (dot products of edge vectors):' /
cbug & ' dota,b,c= ',1p3e22.14 )
cbug write ( 3, 9911) dota, dotb, dotc
cbugc***DEBUG ends.
c=======================================================================********
c.... Find the perpendicular bisectors of the edges (to the center of the
c.... circumscribed circle), and the distances from the edges to the center.
vqabx = -0.5 * dotc * vpabx / (dnorm**2 + small)
vqaby = -0.5 * dotc * vpaby / (dnorm**2 + small)
vqabz = -0.5 * dotc * vpabz / (dnorm**2 + small)
dqab = 0.5 * abs (dotc) * dab / (dnorm + small)
vqbcx = -0.5 * dota * vpbcx / (dnorm**2 + small)
vqbcy = -0.5 * dota * vpbcy / (dnorm**2 + small)
vqbcz = -0.5 * dota * vpbcz / (dnorm**2 + small)
dqbc = 0.5 * abs (dota) * dbc / (dnorm + small)
vqcax = -0.5 * dotb * vpcax / (dnorm**2 + small)
vqcay = -0.5 * dotb * vpcay / (dnorm**2 + small)
vqcaz = -0.5 * dotb * vpcaz / (dnorm**2 + small)
dqca = 0.5 * abs (dotb) * dca / (dnorm + small)
cbugc***DEBUG begins.
cbug 9912 format (/
cbug & 'apttrig results (perpendicular bisectors and distances to',
cbug & ' intersection):' /
cbug & ' vqabx,y,z= ',1p3e22.14 /
cbug & ' dqab= ',1pe22.14 /
cbug & ' vqbcx,y,z= ',1p3e22.14 /
cbug & ' dqbc= ',1pe22.14 /
cbug & ' vqcax,y,z= ',1p3e22.14 /
cbug & ' dqca= ',1pe22.14 )
cbug write ( 3, 9912) vqabx, vqaby, vqabz, dqab,
cbug & vqbcx, vqbcy, vqbcz, dqbc,
cbug & vqcax, vqcay, vqcaz, dqca
cbugc***DEBUG ends.
c=======================================================================********
c.... Find the center and radius of the circumscribed circle.
call aptcirp (ax, ay, az, bx, by, bz, cx, cy, cz, 1, tol,
& ccx, ccy, ccz, vnx, vny, vnz, radcirc, nerra)
cbugc***DEBUG begins.
cbug 9922 format (/
cbug & 'apttrig results (circumscribed circle, center and radius):' /
cbug & ' ccx,y,z= ',1p3e22.14 /
cbug & ' radcirc= ',1pe22.14 /
cbug & ' vnx,vny,vnz= ',1p3e22.14)
cbug write ( 3, 9922) ccx, ccy, ccz, radcirc, vnx, vny, vnz
cbugc***DEBUG ends.
c=======================================================================********
c.... Find the angle bisectors, from the vertex to the opposite edge.
wt1 = dca / (dca + dab + small)
wt2 = -dab / (dca + dab + small)
call aptvsum (0, wt1, abx, aby, abz, wt2, cax, cay, caz, 1, tol,
& viax, viay, viaz, dvia, nerra)
wt1 = dab / (dab + dbc + small)
wt2 = -dbc / (dab + dbc + small)
call aptvsum (0, wt1, bcx, bcy, bcz, wt2, abx, aby, abz, 1, tol,
& vibx, viby, vibz, dvib, nerra)
wt1 = dbc / (dbc + dca + small)
wt2 = -dca / (dbc + dca + small)
call aptvsum (0, wt1, cax, cay, caz, wt2, bcx, bcy, bcz, 1, tol,
& vicx, vicy, vicz, dvic, nerra)
cbugc***DEBUG begins.
cbug 9913 format (/
cbug & 'apttrig results (angle bisectors, lengths):' /
cbug & ' viax,y,z= ',1p3e22.14 /
cbug & ' vibx,y,z= ',1p3e22.14 /
cbug & ' vicx,y,z= ',1p3e22.14 /
cbug & ' lengths= ',1p3e22.14 )
cbug write ( 3, 9913) viax, viay, viaz,
cbug & vibx, viby, vibz,
cbug & vicx, vicy, vicz, dvia, dvib, dvic
cbugc***DEBUG ends.
c=======================================================================********
c.... Find the intercepts of the angle bisectors on the opposite edge.
wt1 = dca / (dca + dab + small)
wt2 = dab / (dca + dab + small)
call aptvsum(0, wt1, bx, by, bz, wt2, cx, cy, cz, 1, tol,
& anax, anay, anaz, dana, nerra)
wt1 = dab / (dab + dbc + small)
wt2 = dbc / (dab + dbc + small)
call aptvsum(0, wt1, cx, cy, cz, wt2, ax, ay, az, 1, tol,
& anbx, anby, anbz, danb, nerra)
wt1 = dbc / (dbc + dca + small)
wt2 = dca / (dbc + dca + small)
call aptvsum(0, wt1, ax, ay, az, wt2, bx, by, bz, 1, tol,
& ancx, ancy, ancz, danc, nerra)
cbugc***DEBUG begins.
cbug 9914 format (/
cbug & 'apttrig results (angle bisector intercepts):' /
cbug & ' anax,y,z= ',1p3e22.14 /
cbug & ' anbx,y,z= ',1p3e22.14 /
cbug & ' ancx,y,z= ',1p3e22.14 )
cbug write ( 3, 9914) anax, anay, anaz,
cbug & anbx, anby, anbz,
cbug & ancx, ancy, ancz
cbugc***DEBUG ends.
c=======================================================================********
c.... Find the center and radius of the inscribed circle.
cix = dbc * ax + dca * bx + dab * cx
ciy = dbc * ay + dca * by + dab * cy
ciz = dbc * az + dca * bz + dab * cz
derrx = tol * (dbc * abs (ax) + dca * abs (bx) + dab * abs (cx))
derry = tol * (dbc * abs (ay) + dca * abs (by) + dab * abs (cy))
derrz = tol * (dbc * abs (az) + dca * abs (bz) + dab * abs (cz))
if (abs (cix) .le. derrx) cix = 0.0
if (abs (ciy) .le. derry) ciy = 0.0
if (abs (ciz) .le. derrz) ciz = 0.0
cix = cix / (dbc + dca + dab + small)
ciy = ciy / (dbc + dca + dab + small)
ciz = ciz / (dbc + dca + dab + small)
radinsc = 2.0 * areatr / (dbc + dca + dab + small)
cbugc***DEBUG begins.
cbug 9915 format (/
cbug & 'apttrig results (center and radius of inscribed circle):' /
cbug & ' cix,y,z= ',1p3e22.14 /
cbug & ' radinsc= ',1pe22.14 )
cbug write ( 3, 9915) cix, ciy, ciz, radinsc
cbugc***DEBUG ends.
c=======================================================================********
210 continue
cbugc***DEBUG begins.
cbug 9916 format (/
cbug & 'apttrig results (ERROR FLAG): nerr = ',i2 )
cbug if (nerr .ne. 0) then
cbug write ( 3, 9916) nerr
cbug endif
cbugc***DEBUG ends.
return
c.... End of subroutine apttrig. (+1 line.)
end
UCRL-WEB-209832