subroutine aptrabc (a, b, c, tol, ncut, acut1, fcut1, dcut1, & acut2, fcut2, dcut2, dcut3, nerr) ccbeg. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c SUBROUTINE APTRABC c c call aptrabc (a, b, c, tol, ncut, acut1, fcut1, dcut1, c & acut2, fcut2, dcut2, dcut3, nerr) c c Version: aptrabc Updated 1999 September 13 14:00. c aptrabc Originated 1999 August 26 16:50. c c Author: Arthur L. Edwards, LLNL, L-298, Telephone (925) 422-4123. c c c Purpose: For the triangle with edges a, b and c, to find the ncut c (1 to 3) straight line(s) that cut the triangle into sections c with equal perimeters and equal areas, and acut1 and acut2, the c names of the edges or vertices on which the endpoints occur, c to find the fractional distances fcut1 and fcut2 of the c endpoints along the edges, and to find the lengths dcut1, dcut2 c and dcut3 of the small triangular section. c Flag nerr indicates any input error. c c Input: a, b, c, tol. c c Output: ncut, acut1, fcut1, dcut1, acut2, fcut2, dcut2, nerr. c c Calls: aptrich, aptqrts c c c Glossary: c c a Input The length of edge a of the triangle. c c acut1 Output The names of the vertices or edges on which the first c endpoint, E, occurs: A, B, C, a, b or c. Size 3. c c acut2 Output The names of the vertices or edges on which the second c endpoint, F, occurs: A, B, C, a, b or c. Size 3. c c b Input The length of edge b of the triangle. c c c Input The length of edge c of the triangle. c c dcut1 Output The distance of the first endput, E, along the edge. c Size 3. c c dcut2 Output The distance of the second endput, F, along the edge. c Size 3. c c dcut3 Output The length of the cutting line, from E to F. Size 3. c c fcut1 Output The fractional distance of the first endpoint, E, c along the edge. Size 3. c c fcut2 Output The fractional distance of the second endpoint, F, c along the edge. Size 3. c c ncut Output The number of cutting lines (1, 2 or 3). c c nerr Output Indicates an input error, if not 0. c 1 if an edge length is negative or too short. c 2 if the triangle is impossible (an edge length is c not less than the sum of the other two edge lengths. c c tol Input Numerical tolerance limit. c On Cray computers, recommend 1.e-5 to 1.e-11. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccend. implicit none c.... Arguments. real a ! Length of edge a of the triangle. real b ! Length of edge b of the triangle. real c ! Length of edge c of the triangle. real dcut1(3) ! Distance of point E along edge. real dcut2(3) ! Distance of point F along edge. real dcut3(3) ! Length of cutting line. real fcut1(3) ! Fractional distance of E along edge. real fcut2(3) ! Fractional distance of F along edge. integer ncut ! Number of equipartition cutting lines. integer nerr ! Error indicator. real tol ! Numerical tolerence limit. character*1 acut1(3) ! Names of edges or vertices for point E. character*1 acut2(3) ! Names of edges or vertices for point F. c.... Local variables. real aa ! Coefficient of t**2 in quadratic eq. real bb ! Coefficient of t in quadratic eq. real cc ! Coefficient of 1 in quadratic eq. real fab ! Fractional distance from A to B. real fac ! Fractional distance from A to C. real fba ! Fractional distance from B to A. real fbc ! Fractional distance from B to C. real fca ! Fractional distance from C to A. real fcb ! Fractional distance from C to B. integer itrun ! Truncation flag from aptqrts. integer kvera ! Cuts through vertex A. integer kverb ! Cuts through vertex B. integer kverc ! Cuts through vertex C. integer n ! Index of cut line. integer nerra ! Error flag from apt subroutine. integer nroots ! # of real roots of quadratic equation. real per ! Perimeter of triangle, c + a + b. real qq ! bb**2 - 4*aa*bb. real root1 ! First real root of quadatic equation. real root2 ! Second real root of quadratic equation. character*1 acut1x ! Names of edges or vertices for point E. character*1 acut2x ! Names of edges or vertices for point F. cbugc***DEBUG begins. cbug 9901 format (/ 'aptrabc finding triangle partition cuts.', cbug & ' tol=',1pe22.14 / cbug & ' a ,b ,c = ',1p3e22.14 ) cbug write ( 3, 9901) tol, a, b, c cbugc***DEBUG ends. c.... Initialize. nerr = 0 ncut = 0 kvera = 0 kverb = 0 kverc = 0 do 110 n = 1, 3 acut1(n) = '?' acut2(n) = '?' fcut1(n) = -1.e99 fcut2(n) = -1.e99 110 continue c.... Test for input errors. call aptrich (a, b, c, tol, nerr) if (nerr .ne. 0) go to 410 c.... Find lines that cut triangle into equal areas, perimeters. per = a + b + c c=======================================================================******** c.... Try vertex C, from edge a to edge b. aa = 2.0 * a bb = -per cc = b call aptqrts (0, aa, bb, cc, qq, tol, & nroots, root1, root2, itrun) if (nroots .ge. 1) then if ((root1 .ge. (0.5 - tol)) .and. & (root1 .le. (1.0 + tol))) then acut1x = 'a' acut2x = 'b' fcb = root1 fca = 0.5 / fcb fbc = 1.0 - fcb if (abs (fbc) .le. tol) then if (kverb .gt. 0) go to 150 acut1x = 'B' kverb = 1 fbc = 0.0 fcb = 1.0 fac = 0.5 fca = 0.5 endif fac = 1.0 - fca if (abs (fac) .le. tol) then if (kvera .gt. 0) go to 150 acut2x = 'A' kvera = 1 fac = 0.0 fca = 1.0 fbc = 0.5 fcb = 0.5 endif ncut = ncut + 1 acut1(ncut) = acut1x acut2(ncut) = acut2x fcut1(ncut) = fcb fcut2(ncut) = fca dcut1(ncut) = fcb * a dcut2(ncut) = fca * b dcut3(ncut) = sqrt (per * (c - 0.25 * per)) endif ! Cut endpoints are on edges. endif ! Cut endpoints are real. 150 if (nroots .eq. 2) then if ((root2 .ge. (0.5 - tol)) .and. & (root2 .le. (1.0 + tol))) then acut1x = 'a' acut2x = 'b' fcb = root2 fca = 0.5 / fcb fbc = 1.0 - fcb if (abs (fbc) .le. tol) then if (kverb .gt. 0) go to 190 kverb = 1 acut1x = 'B' fbc = 0.0 fcb = 1.0 fca = 0.5 fac = 0.5 endif fac = 1.0 - fca if (abs (fac) .le. tol) then if (kvera .gt. 0) go to 190 kvera = 1 acut2x = 'A' fac = 0.0 fca = 1.0 fbc = 0.5 fcb = 0.5 endif ncut = ncut + 1 acut1(ncut) = acut1x acut2(ncut) = acut2x fcut1(ncut) = fcb fcut2(ncut) = fca dcut1(ncut) = fcb * a dcut2(ncut) = fca * b dcut3(ncut) = sqrt (per * (c - 0.25 * per)) endif ! Cut endpoints are on edges. endif ! Cut endpoints are real. 190 continue c=======================================================================******** c.... Try vertex A, from edge b to edge c. aa = 2.0 * b bb = -per cc = c call aptqrts (0, aa, bb, cc, qq, tol, & nroots, root1, root2, itrun) if (nroots .ge. 1) then if ((root1 .ge. (0.5 - tol)) .and. & (root1 .le. (1.0 + tol))) then acut1x = 'b' acut2x = 'c' fac = root1 fab = 0.5 / fac fca = 1.0 - fac if (abs (fca) .le. tol) then if (kverc .gt. 0) go to 250 kverc = 1 acut1x = 'C' fac = 1.0 fca = 0.0 fab = 0.5 fba = 0.5 endif fba = 1.0 - fab if (abs (fba) .le. tol) then if (kverb .gt. 0) go to 250 kverb = 1 acut2x = 'B' fab = 1.0 fba = 0.0 fac = 0.5 fca = 0.5 endif ncut = ncut + 1 acut1(ncut) = acut1x acut2(ncut) = acut2x fcut1(ncut) = fac fcut2(ncut) = fab dcut1(ncut) = fac * b dcut2(ncut) = fab * c dcut3(ncut) = sqrt (per * (a - 0.25 * per)) endif ! Cut endpoints are on edges. endif ! Cut endpoints are real. 250 if (nroots .eq. 2) then if ((root2 .ge. (0.5 - tol)) .and. & (root2 .le. (1.0 + tol))) then acut1x = 'b' acut2x = 'c' fac = root2 fab = 0.5 / fac fca = 1.0 - fac if (abs (fca) .le. tol) then if (kverc .gt. 0) go to 290 kverc = 1 acut1x = 'C' fac = 1.0 fca = 0.0 fab = 0.5 fba = 0.5 endif fba = 1.0 - fab if (abs (fba) .le. tol) then if (kverb .gt. 0) go to 290 kverb = 1 acut2x = 'B' fab = 1.0 fba = 0.0 fac = 0.5 fca = 0.5 endif ncut = ncut + 1 acut1(ncut) = acut1x acut2(ncut) = acut2x fcut1(ncut) = fac fcut2(ncut) = fab dcut1(ncut) = fac * b dcut2(ncut) = fab * c dcut3(ncut) = sqrt (per * (a - 0.25 * per)) endif ! Cut endpoints are on edges. endif ! Cut endpoints are real. 290 continue c=======================================================================******** c.... Try vertex B, from edge c to edge a. aa = 2.0 * c bb = -per cc = a call aptqrts (0, aa, bb, cc, qq, tol, & nroots, root1, root2, itrun) if (nroots .ge. 1) then if ((root1 .ge. (0.5 - tol)) .and. & (root1 .le. (1.0 + tol))) then acut1x = 'c' acut2x = 'a' fba = root1 fbc = 0.5 / fba fab = 1.0 - fba if (abs (fab) .le. tol) then if (kvera .gt. 0) go to 350 kvera = 1 acut1x = 'A' fab = 0.0 fba = 1.0 fbc = 0.5 fcb = 0.5 endif fcb = 1.0 - fbc if (abs (fcb) .le. tol) then if (kverc .gt. 0) go to 350 acut2x = 'C' kverc = 1 fbc = 1.0 fcb = 0.0 fab = 0.5 fba = 0.5 endif ncut = ncut + 1 acut1(ncut) = acut1x acut2(ncut) = acut2x fcut1(ncut) = fba fcut2(ncut) = fbc dcut1(ncut) = fba * c dcut2(ncut) = fbc * a dcut3(ncut) = sqrt (per * (b - 0.25 * per)) endif ! Cut endpoints are on edges. endif ! Cut endpoints are real. 350 if (nroots .eq. 2) then if ((root2 .ge. (0.5 - tol)) .and. & (root2 .le. (1.0 + tol))) then acut1x = 'c' acut2x = 'a' fba = root2 fbc = 0.5 / fba fab = 1.0 - fba if (abs (fab) .le. tol) then if (kvera .gt. 0) go to 390 kvera = 1 acut1x = 'A' fab = 0.0 fba = 1.0 fbc = 0.5 fcb = 0.5 endif fcb = 1.0 - fbc if (abs (fcb) .le. tol) fcb = 0.0 if (abs (fcb) .le. tol) then if (kverc .gt. 0) go to 390 kverc = 1 acut2x = 'C' fbc = 1.0 fcb = 0.0 fab = 0.5 fba = 0.5 endif ncut = ncut + 1 acut1(ncut) = acut1x acut2(ncut) = acut2x fcut1(ncut) = fba fcut2(ncut) = fbc dcut1(ncut) = fba * c dcut2(ncut) = fbc * a dcut3(ncut) = sqrt (per * (b - 0.25 * per)) endif ! Cut endpoints are on edges. endif ! Cut endpoints are real. 390 continue c=======================================================================******** 410 continue cbugc***DEBUG begins. cbug 9916 format (/ 'aptrabc results. nerr = ',i2,' ncut = ',i1 ) cbug write ( 3, 9916) nerr, ncut cbug cbug 9918 format ('ncut = ',i2,' acut1 = ',a1,' acut2 = ',a1 / cbug & 'fcut1, fcut2 =',1p2e22.14 / cbug & 'dcut1,2,3 =',1p3e22.14 ) cbug cbug do 420 n = 1, ncut cbug write ( 3, 9918) n, acut1(n), acut2(n), cbug & fcut1(n), fcut2(n), cbug & dcut1(n), dcut2(n), dcut3(n) cbug 420 continue cbugc***DEBUG ends. return c.... End of subroutine aptrabc. (+1 line.) end UCRL-WEB-209832