program main !*****************************************************************************80 ! !! geompack2_test() tests geompack2(). ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 12 December 2025 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! This version by John Burkardt. ! implicit none call timestamp ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'geompack2_test():' write ( *, '(a)' ) ' FORTRAN90 version' write ( *, '(a)' ) ' Test geompack2().' call test01 ( ) call test02 ( ) call test03 ( 'cmos.in' ) call test04 ( ) call test05 ( 'annulus.in' ) call test06 ( ) call test07 ( ) call test08 ( 'annulus.in' ) call test09 ( ) call test10 ( ) call test11 ( 'annulus.in' ) call test12 ( ) call test13 ( 'ptpg.in' ) call test14 ( ) call test15 ( 'shr1.in' ) call test16 ( ) call test17 ( 'annulus.in' ) call test18 ( ) call test19 ( ) call test20 ( ) call test21 ( ) ! ! Terminate. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'geompack2_test():' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) ' ' call timestamp ( ) stop 0 end subroutine test01 ( ) !*****************************************************************************80 ! !! TEST01 tests ANGLE, AREAPG, AREATR. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 12 December 2025 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! This version by John Burkardt. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) integer, parameter :: n = 9 real ( kind = rk8 ) ang(n) real ( kind = rk8 ) angle real ( kind = rk8 ) area real ( kind = rk8 ) areapg real ( kind = rk8 ) areatr real ( kind = rk8 ) atr(n) integer i real ( kind = rk8 ) sum2 real ( kind = rk8 ), dimension ( n ) :: xc = (/ & 5.0D+00, 5.0D+00, 7.0D+00, 9.0D+00, 5.0D+00, & 1.0D+00, 2.0D+00, 0.0D+00, 3.0D+00 /) real ( kind = rk8 ), dimension ( n ) :: yc = (/ & 0.0D+00, 2.0D+00, 2.0D+00, 4.0D+00, 8.0D+00, & 7.0D+00, 5.0D+00, 3.0D+00, 3.0D+00 /) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST01' write ( *, '(a)' ) ' AREAPG computes the area of a polygon;' write ( *, '(a)' ) ' AREATR computes the area of a triangle;' write ( *, '(a)' ) ' ANGLE computes the polygonal angle at any vertex.' write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) ' Number of polygonal vertices is N = ', n write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' I, XC(I), YC(I)' write ( *, '(a)' ) ' ' do i = 1, n write ( *, '(i5,2d15.7)' ) i, xc(i), yc(i) end do area = areapg ( n, xc, yc ) write ( *, '(a)' ) ' ' write ( *, '(a,g14.6)' ) ' Area computed directly by AREAPG = ', area sum2 = 0.0D+00 do i = 1, n if ( i == 1 ) then atr(i) = areatr ( xc(1), yc(1), xc(i), yc(i), xc(i+1), yc(i+1) ) ang(i) = angle ( xc(n), yc(n), xc(i), yc(i), xc(i+1), yc(i+1) ) else if ( i == n ) then atr(i) = areatr ( xc(1), yc(1), xc(i), yc(i), xc(1), yc(1) ) ang(i) = angle ( xc(i-1), yc(i-1), xc(i), yc(i), xc(1), yc(1) ) else atr(i) = areatr ( xc(1), yc(1), xc(i), yc(i), xc(i+1), yc(i+1) ) ang(i) = angle ( xc(i-1), yc(i-1), xc(i), yc(i), xc(i+1), yc(i+1) ) end if sum2 = sum2 + atr(i) end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I, AREATR(I), ANGLE(I)' write ( *, '(a)' ) ' ' do i = 1, n write ( *, '(i5,2d15.7)' ) i, atr(i), ang(i) end do write ( *, '(a,g14.6)' ) & 'Area computed indirectly by summing AREATR = ', sum2 return end subroutine test02 ( ) !*****************************************************************************80 ! !! TEST02 tests CMCIRC. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 12 December 2025 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! This version by John Burkardt. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) integer cmcirc integer in real ( kind = rk8 ), parameter :: x0 = 3.0D+00 real ( kind = rk8 ), parameter :: y0 = 3.0D+00 real ( kind = rk8 ), parameter :: x1 = 5.0D+00 real ( kind = rk8 ), parameter :: y1 = 0.0D+00 real ( kind = rk8 ), parameter :: x2 = 0.0D+00 real ( kind = rk8 ), parameter :: y2 = 5.0D+00 real ( kind = rk8 ), parameter :: x3 = -5.0D+00 real ( kind = rk8 ), parameter :: y3 = 0.0D+00 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST02' write ( *, '(a)' ) ' CMCIRC determines if a point lies in, on' write ( *, '(a)' ) ' or outside a circle given by 3 points.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The points defining the circle:' write ( *, '(a)' ) ' ' write ( *, '(a,2g14.6)' ) ' X1,Y1 = ', x1, y1 write ( *, '(a,2g14.6)' ) ' X2,Y2 = ', x2, y2 write ( *, '(a,2g14.6)' ) ' X3,Y3 = ', x3, y3 write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The point to be tested:' write ( *, '(a)' ) ' ' write ( *, '(a,2g14.6)' ) ' X0,Y0 = ', x0, y0 in = cmcirc ( x0, y0, x1, y1, x2, y2, x3, y3 ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Test results:' if ( in == 2 ) then write ( *, '(a)' ) ' The three points are collinear.' else if ( in == 1 ) then write ( *, '(a)' ) ' The point is inside the circle.' else if ( in == 0 ) then write ( *, '(a)' ) ' The point is on the circle.' else if ( in == -1 ) then write ( *, '(a)' ) ' The point is outside the circle.' end if return end subroutine test03 ( filename ) !*****************************************************************************80 ! !! TEST03 tests CVDEC2, FNDSEP, INSED2, INSVR2, JNHOLE, MINANG, RESVRT, SPDEC2. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 12 December 2025 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! This version by John Burkardt. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) integer, parameter :: incr = 10000 integer, parameter :: inunit = 1 integer, parameter :: maxed = 101 integer, parameter :: maxhv = 200 integer, parameter :: maxiw = 900 integer, parameter :: maxnc = 30 integer, parameter :: maxnv = 400 integer, parameter :: maxpv = 1000 integer, parameter :: maxvc = 500 integer, parameter :: maxwk = 1500 real ( kind = rk8 ) angmin real ( kind = rk8 ) angspc real ( kind = rk8 ) angtol integer case real ( kind = rk8 ) degrees_to_radians integer edge(4,maxed) character ( len = * ) filename integer ht(0:maxed-1) integer htsiz integer hvl(maxhv) integer i real ( kind = rk8 ) iang(maxpv) integer icur(maxnc) integer ierror integer ivrt(maxnv) integer iwk(maxiw) integer j integer map(maxnc) integer msglvl integer ncur integer nh integer nhola integer nhole integer npolg integer nrfv integer nsc integer nv integer nvbc(maxnc) integer nvc integer nvcin integer nvert real ( kind = rk8 ), parameter :: pi = 3.141592653589793D+00 integer prime integer pvl(4,maxpv) real ( kind = rk8 ) radians_to_degrees integer regnum(maxhv) character ( len = 20 ) rgname real ( kind = rk8 ) tol real ( kind = rk8 ) tolin real ( kind = rk8 ) vcl(2,maxvc) real ( kind = rk8 ) wk(maxwk) tol = 100.0D+00 * epsilon ( tol ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST03' write ( *, '(a)' ) ' Test CVDEC.' write ( *, '(a)' ) ' Test FNDSEP.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Reading input file: "' // trim ( filename ) // '".' open ( unit = inunit, file = filename, form = 'formatted' ) ! ! Read in vertices of general polygonal region. ! ! CASE = 1 : simple polygon or multiply connected polygonal region ! CASE = 2 : general polygonal region with holes and interfaces ! CASE = -1 or -2: similar, but decompose into simple polygons only ! and don't obtain convex polygon decomposition ! write ( *, '(a)' ) ' ' read ( inunit, '(a)' ) rgname write ( *, '(a)' ) trim ( rgname ) read ( inunit, * ) tolin, angspc, angtol write ( *, '(a,d15.7)' ) ' TOLIN = ', tolin write ( *, '(a,d15.7)' ) ' ANGSPC = ', angspc write ( *, '(a,d15.7)' ) ' ANGTOL = ', angtol angspc = degrees_to_radians ( angspc ) angtol = degrees_to_radians ( angtol ) read ( inunit, * ) case, nvc, ncur, msglvl if ( maxvc < nvc ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST03 - Error!' write ( *, '(a)' ) ' MAXVC < NVC.' return end if if ( maxnc < ncur ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST03 - Error!' write ( *, '(a)' ) ' MAXNC < NCUR.' return end if read ( inunit, * ) nvbc(1:ncur) if ( abs(case) == 2 ) then read ( inunit, * ) icur(1:ncur) end if read ( inunit, * ) ( vcl(1,i), vcl(2,i), i = 1, nvc ) ! ! Call DSMCPR or DSPGDC to set data structures in arrays ! REGNUM, HVL, PVL, IANG, HOLV = IWK. ! if ( abs ( case ) == 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Calling DSMCPR to set data structures.' nhole = ncur - 1 call dsmcpr ( nhole, nvbc, vcl, maxhv, maxpv, maxiw, nvc, npolg, nvert, & nhola, regnum, hvl, pvl, iang, iwk, ierror ) else if ( abs(case) == 2 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Calling DSPGDC to set data structures.' nv = sum ( nvbc(1:ncur) ) if ( maxnv < nv ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST03 - Error!' write ( *, '(a)' ) ' MAXNV < NV.' return end if read ( inunit, * ) ivrt(1:nv) nsc = 0 do i = 1, nv if ( ivrt(i) < 0 ) then nsc = nsc + 1 end if end do nsc = nsc / 2 if ( maxed < nsc ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST03 - Fatal error!' write ( *, '(a)' ) ' MAXED < NSC.' return end if htsiz = min ( prime ( nsc / 2 ), maxed ) call dspgdc ( nvc, vcl, incr, ncur, nvbc, icur, ivrt, maxhv, maxpv, maxiw, & npolg, nvert, nhole, nhola, regnum, hvl, pvl, iang, iwk, htsiz, nsc, & ht, edge, map, ierror ) end if close ( unit = inunit ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST03 - Error!' write ( *, '(a)' ) ' IERROR = ', ierror return end if ! ! Print out data structures and measurements after initialization. ! nvcin = nvc nh = nhole * 2 + nhola write ( *, '(a,i6)' ) 'MSGLVL = ', msglvl write ( *, '(a,i6)' ) 'NVC = ', nvc write ( *,630) (i,vcl(1,i),vcl(2,i),i=1,nvc) write ( *,640) npolg,nhole,(hvl(i),i=1,npolg+nhole) write ( *,650) (regnum(i),i=1,npolg) write ( *,660) nvert,(i,(pvl(j,i),j=1,4),iang(i),i=1,nvert) write ( *,640) nhola,nh,(iwk(i),i=1,nh) nrfv = 0 do i = 1, nvert if ( pi + tol < iang(i) ) then nrfv = nrfv + 1 end if end do angmin = iang(1) do i = 2, nvert angmin = min ( angmin, iang(i) ) end do angmin = radians_to_degrees ( angmin ) write (*,710) nvc,npolg,nvert,nhole,nhola,nrfv,angmin ! ! Obtain simple (and convex) polygon decompositions. ! if ( msglvl == 2 ) then write ( *,670) end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Calling SPDEC2' call spdec2 ( angspc, angtol, nvc, npolg, nvert, nhole, nhola, maxvc, & maxhv, maxpv, maxiw-nh, maxwk, iwk, vcl, regnum, hvl, pvl, iang, & iwk(nh+1), wk, ierror ) if ( 0 < case ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Calling CVDEC2' call cvdec2 ( angspc, angtol, nvc, npolg, nvert, maxvc, maxhv, maxpv, & maxiw, maxwk, vcl, regnum, hvl, pvl, iang, iwk, wk, ierror ) end if if ( msglvl == 2 ) then write ( *,670) 0, 0, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00 end if if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST03 - Error!' write ( *, '(a)' ) ' IERROR = ', ierror return end if ! ! Print out data structures and measurements after decomposition. ! write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) ' NVC = ', nvc write ( *,630) (i,vcl(1,i),vcl(2,i),i=nvcin+1,nvc) write ( *,680) npolg, hvl(1:npolg) write ( *,650) regnum(1:npolg) write ( *,660) nvert,(i,(pvl(j,i),j=1,4),iang(i),i=1,nvert) angmin = minval ( iang(1:nvert) ) angmin = radians_to_degrees ( angmin ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DECOMP:' write ( *, '(a,i6)' ) ' NVC = ', nvc write ( *, '(a,i6)' ) ' NPOLG = ', npolg write ( *, '(a,i6)' ) ' NVERT = ', nvert write ( *, '(a,g14.6)' ) ' ANGMIN = ', angmin 630 format ((1x,i7,2f15.7)) 640 format (/1x,2i7/(1x,10i7)) 650 format (/(1x,10i7)) 660 format (/1x,i7/(1x,5i7,f15.7)) 670 format (1x,2i7,4f15.7) 680 format (/1x,i7/(1x,10i7)) 710 format (1x,'initds: nvc=',i7,' npolg=',i7,' nvert=',i7, & ' nhole=',i7/9x,'nhola=',i7,' nrfv=',i7,' angmin=',f9.3 ) return end subroutine test04 ( ) !*****************************************************************************80 ! !! TEST04 tests DIAEDG. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 12 December 2025 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! This version by John Burkardt. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) integer diaedg integer in real ( kind = rk8 ) x0 real ( kind = rk8 ) x1 real ( kind = rk8 ) x2 real ( kind = rk8 ) x3 real ( kind = rk8 ) y0 real ( kind = rk8 ) y1 real ( kind = rk8 ) y2 real ( kind = rk8 ) y3 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST04' write ( *, '(a)' ) ' DIAEDG determines which diagonal of a' write ( *, '(a)' ) ' quadrilateral is to be preferred, based on' write ( *, '(a)' ) ' the circumcircle criterion.' x0 = 0.0D+00 y0 = 0.0D+00 x1 = 5.0D+00 y1 = 0.0D+00 x2 = 6.0D+00 y2 = 1.0D+00 x3 = 1.0D+00 y3 = 1.0D+00 write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The points defining the quadrilateral:' write ( *, '(a)' ) ' ' write ( *, '(a,2g14.6)' ) ' P0: X0,Y0 = ', x0, y0 write ( *, '(a,2g14.6)' ) ' P1: X1,Y1 = ', x1, y1 write ( *, '(a,2g14.6)' ) ' P2: X2,Y2 = ', x2, y2 write ( *, '(a,2g14.6)' ) ' P3: X3,Y3 = ', x3, y3 in = diaedg ( x0, y0, x1, y1, x2, y2, x3, y3 ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' DIAEDG results:' write ( *, '(a)' ) ' ' if ( in == 1 ) then write ( *, '(a)' ) ' Use diagonal P0--P2.' else if ( in == -1 ) then write ( *, '(a)' ) ' Use diagonal P1--P3.' else if ( in == 0 ) then write ( *, '(a)' ) ' All 4 points lie on a circle.' write ( *, '(a)' ) ' Either diagonal can be used.' end if return end subroutine test05 ( filename ) !*****************************************************************************80 ! !! TEST05 tests DSMCPR, DSPGDC, EDGHT, HOLVRT. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 12 December 2025 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! This version by John Burkardt. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) integer, parameter :: incr = 10000 integer, parameter :: maxed = 101 integer, parameter :: maxho = 50 integer, parameter :: maxnc = 30 integer, parameter :: maxnv = 500 integer case integer edge(4,maxed) character ( len = * ) filename integer holv(maxho) integer ht(0:maxed-1) integer htsiz integer hvl(maxnc*2) integer i real ( kind = rk8 ) iang(maxnv*2) integer icur(maxnc) integer ierror integer, parameter :: inunit = 1 integer ivrt(maxnv) integer j integer map(maxnc) integer ncur integer nhola integer nhole integer npolg integer nsc integer nv integer nvbc(maxnc) integer nvc integer nvert integer prime integer pvl(4,maxnv*2) integer regnum(maxnc*2) character ( len = 20 ) rgname real ( kind = rk8 ) tolin real ( kind = rk8 ) vcl(2,maxnv) ! ! Read in the vertices of a general polygonal region. ! CASE = 1 : simple polygon or multiply connected polygonal region ! CASE = 2 : general polygonal region with holes and interfaces ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST05' write ( *, '(a)' ) ' Reading input file:' write ( *, '(a)' ) filename open ( unit = inunit, file = filename, form = 'formatted' ) read ( inunit, '(a)' ) rgname read ( inunit, * ) tolin read ( inunit, * ) case, nvc, ncur if ( maxnv < nvc ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST05 - Error!' write ( *, '(a)' ) ' MAXNV < NVC.' return end if if ( maxnc < ncur ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST05 - Error!' write ( *, '(a)' ) ' MAXNC < NCUR.' return end if read ( inunit, * ) nvbc(1:ncur) if ( case == 2) then read ( inunit, * ) icur(1:ncur) end if read ( inunit, * ) ( vcl(1,i), vcl(2,i), i = 1, nvc ) ! ! Call DSMCPR or DSPGDC to set data structures in arrays ! REGNUM, HVL, PVL, IANG, HOLV. ! if ( case == 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Calling DSMCPR:' nhole = ncur - 1 call dsmcpr ( nhole, nvbc, vcl, maxnc*2, maxnv*2, maxho, nvc, & npolg, nvert, nhola, regnum, hvl, pvl, iang, holv, ierror ) else if ( case == 2 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Calling DSPGDC:' nv = sum ( nvbc(1:ncur) ) if ( maxnv < nv ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST05 - Fatal error!' write ( *, '(a)' ) ' MAXNV < NV.' return end if read ( inunit, * ) ivrt(1:nv) nsc = 0 do i = 1, nv if ( ivrt(i) < 0 ) then nsc = nsc + 1 end if end do nsc = nsc / 2 if ( maxed < nsc ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST05 - Fatal error!' write ( *, '(a)' ) ' MAXED < NSC.' return end if htsiz = min ( prime(nsc/2), maxed ) call dspgdc ( nvc, vcl, incr, ncur, nvbc, icur, ivrt, maxnc*2, maxnv*2, & maxho, npolg, nvert, nhole, nhola, regnum, hvl, pvl, iang, holv, & htsiz, nsc, ht, edge, map, ierror ) end if close ( unit = inunit ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST05 - Error!' write ( *, '(a)' ) ' IERROR = ', ierror return end if ! ! Print out the data structures. ! write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) ' NVC = ', nvc write ( *, '(a)' ) ' ' write ( *, '(i7,2f15.7)' ) (i,vcl(1:2,i),i=1,nvc) write ( *, '(a)' ) ' ' write ( *,640) npolg,nhole,(hvl(i),i=1,npolg+nhole) write ( *, '(a)' ) ' ' write ( *,650) regnum(1:npolg) write ( *, '(a)' ) ' ' write ( *,660) nvert,(i,(pvl(j,i),j=1,4),iang(i),i=1,nvert) write ( *, '(a)' ) ' ' write ( *,640) nhola,nhole*2+nhola,(holv(i),i=1,nhole*2+nhola) 640 format (1x,2i7/(1x,10i7)) 650 format ((1x,10i7)) 660 format (1x,i7/(1x,5i7,f15.7)) return end subroutine test06 ( ) !*****************************************************************************80 ! !! TEST06 tests DTRIS2; ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 12 December 2025 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! This version by John Burkardt. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) integer, parameter :: g_num = 24 integer a integer b integer c integer d integer diaedg character ( len = 255 ) :: file_name = 'test06_triangulation_plot.eps' real ( kind = rk8 ), dimension (2,g_num) :: g_xy = reshape ( (/ & 1.0000000D+00, 0.0000000D+00, & 0.9238795D+00, 0.3826834D+00, & 0.7071068D+00, 0.7071068D+00, & 0.3826834D+00, 0.9238795D+00, & 0.0000000D+00, 1.0000000D+00, & -0.3826834D+00, 0.9238795D+00, & -0.7071068D+00, 0.7071068D+00, & -0.9238795D+00, 0.3826834D+00, & -1.0000000D+00, 0.0000000D+00, & -0.9238795D+00, -0.3826834D+00, & -0.7071068D+00, -0.7071068D+00, & -0.3826834D+00, -0.9238795D+00, & 0.0000000D+00, -1.0000000D+00, & 0.3826834D+00, -0.9238795D+00, & 0.7071068D+00, -0.7071068D+00, & 0.9238795D+00, -0.3826834D+00, & 0.7500000D+00, 0.0000000D+00, & 0.6767767D+00, 0.1767767D+00, & 0.5000000D+00, 0.2500000D+00, & 0.3232233D+00, 0.1767767D+00, & 0.2500000D+00, 0.0000000D+00, & 0.3232233D+00, -0.1767767D+00, & 0.5000000D+00, -0.2500000D+00, & 0.6767767D+00, -0.1767767D+00 /), (/ 2, g_num /) ) integer i integer ierror integer ind(g_num) integer j integer jp1 integer jp2 integer k integer nlo integer nod_tri(3,2*g_num) integer tnbr(3,2*g_num) integer tri_num do i = 1, g_num ind(i) = i end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST06' write ( *, '(a)' ) ' DTRIS2 constructs the Delaunay triangulation' write ( *, '(a)' ) ' of a set of points in the plane.' write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) ' The number of points to triangulate is ', g_num write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The coordinates of the points are:' write ( *, '(a)' ) ' ' do i = 1, g_num write ( *, '(i5,2f15.7)' ) i, g_xy(1,i), g_xy(2,i) end do call dtris2 ( g_num, g_xy, ind, tri_num, nod_tri, tnbr, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST06 - Error!' write ( *, '(a,i6)' ) ' IERROR = ', ierror return end if nlo = 0 do i = 1, tri_num do j = 1, 3 k = tnbr(j,i) if ( i < k ) then jp1 = j + 1 if ( 3 < jp1 ) then jp1 = 1 end if jp2 = jp1 + 1 if ( 3 < jp2 ) then jp2 = 1 end if a = nod_tri(j,i) b = nod_tri(jp1,i) c = nod_tri(jp2,i) if ( nod_tri(1,k) == b ) then d = nod_tri(3,k) else if ( nod_tri(2,k) == b ) then d = nod_tri(1,k) else d = nod_tri(2,k) end if if ( diaedg ( g_xy(1,c), g_xy(2,c), g_xy(1,a), g_xy(2,a), g_xy(1,d), & g_xy(2,d), g_xy(1,b), g_xy(2,b) ) == 1 ) then nlo = nlo + 1 end if end if end do end do write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) ' The number of triangles = ', tri_num write ( *, '(a,i6)' ) ' NLO = ', nlo write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' I, NOD_TRI(1:3,I), TNBR(1:3,I)' write ( *, '(a)' ) ' ' do i = 1, tri_num write ( *, '(7i8)' ) i, nod_tri(1:3,i), tnbr(1:3,i) end do call triangulation_plot_eps ( file_name, g_num, g_xy, tri_num, nod_tri ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' TRIANGULATION_PLOT_EPS has created an EPS file' write ( *, '(a)' ) ' containing an image of the triangulation, ' write ( *, '(a)' ) ' in ' // trim ( file_name ) return end subroutine test07 ( ) !*****************************************************************************80 ! !! TEST07 tests DTRIW2; ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 12 December 2025 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! This version by John Burkardt. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) real ( kind = rk8 ), parameter :: large = 1000.0D+00 integer, parameter :: maxnp = 10000 integer, parameter :: maxst = 100 integer a integer alg integer b real ( kind = rk8 ) binexp integer c integer d integer diaedg integer i integer ierror integer ind(maxnp+3) integer j integer jp1 integer jp2 integer k integer msglvl integer nlo integer npt integer ntri integer stack(maxst) integer til(3,maxnp*2+1) integer tnbr(3,maxnp*2+1) real ( kind = rk8 ) vcl(2,maxnp+3) ! ! ALG = ! 2: DTRIW2; ! 3: DTRIW2 with bounding triangle; ! 4: DTRIW2 with call to BNSRT2 first. ! MSGLVL ! 0: print arrays; ! 4: also print edges as they are created and swapped. ! ! I HAVE NO IDEA HOW TO CHOOSE BINEXP ! msglvl = 0 binexp = 0.5D+00 vcl(1,1) = 1.0000000D+00 vcl(2,1) = 0.0000000D+00 vcl(1,2) = 0.9238795D+00 vcl(2,2) = 0.3826834D+00 vcl(1,3) = 0.7071068D+00 vcl(2,3) = 0.7071068D+00 vcl(1,4) = 0.3826834D+00 vcl(2,4) = 0.9238795D+00 vcl(1,5) = 0.0000000D+00 vcl(2,5) = 1.0000000D+00 vcl(1,6) = - 0.3826834D+00 vcl(2,6) = 0.9238795D+00 vcl(1,7) = - 0.7071068D+00 vcl(2,7) = 0.7071068D+00 vcl(1,8) = - 0.9238795D+00 vcl(2,8) = 0.3826834D+00 vcl(1,9) = - 1.0000000D+00 vcl(2,9) = 0.0000000D+00 vcl(1,10) = - 0.9238795D+00 vcl(2,10) = - 0.3826834D+00 vcl(1,11) = - 0.7071068D+00 vcl(2,11) = - 0.7071068D+00 vcl(1,12) = - 0.3826834D+00 vcl(2,12) = - 0.9238795D+00 vcl(1,13) = 0.0000000D+00 vcl(2,13) = - 1.0000000D+00 vcl(1,14) = 0.3826834D+00 vcl(2,14) = - 0.9238795D+00 vcl(1,15) = 0.7071068D+00 vcl(2,15) = - 0.7071068D+00 vcl(1,16) = 0.9238795D+00 vcl(2,16) = - 0.3826834D+00 vcl(1,17) = 0.7500000D+00 vcl(2,17) = 0.0000000D+00 vcl(1,18) = 0.6767767D+00 vcl(2,18) = 0.1767767D+00 vcl(1,19) = 0.5000000D+00 vcl(2,19) = 0.2500000D+00 vcl(1,20) = 0.3232233D+00 vcl(2,20) = 0.1767767D+00 vcl(1,21) = 0.2500000D+00 vcl(2,21) = 0.0000000D+00 vcl(1,22) = 0.3232233D+00 vcl(2,22) = - 0.1767767D+00 vcl(1,23) = 0.5000000D+00 vcl(2,23) = - 0.2500000D+00 vcl(1,24) = 0.6767767D+00 vcl(2,24) = - 0.1767767D+00 npt = 24 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST07' write ( *, '(a,i6)' ) ' MSGLVL = ', msglvl write ( *, '(a,i6)' ) ' NPT = ', npt write ( *, '(a,g14.6)' ) ' BINEXP = ', binexp if ( maxnp < npt ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST07 - Fatal error!' write ( *, '(a)' ) ' MAXNP < NPT.' return end if write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) ' The number of points to triangulate is ', npt write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The coordinates of the points are:' write ( *, '(a)' ) ' ' do i = 1, npt write ( *, '(i5,2f15.7)' ) i, vcl(1,i), vcl(2,i) end do do alg = 2, 4 npt = 24 write ( *, '(a,i6)' ) 'ALG = ', alg if ( alg /= 3 ) then do i = 1, npt ind(i) = i end do else vcl(1,npt+1) = -large vcl(2,npt+1) = -large vcl(1,npt+2) = large vcl(2,npt+2) = -large vcl(1,npt+3) = 0.0D+00 vcl(2,npt+3) = large ind(1) = npt + 1 ind(2) = npt + 2 ind(3) = npt + 3 do i = 1, npt ind(i+3) = i end do npt = npt + 3 end if if ( alg == 4 ) then call bnsrt2 ( binexp, npt, vcl, ind, til, tnbr ) end if call dtriw2 ( npt, maxst, vcl, ind, ntri, til, tnbr, stack, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST07 - Error!' write ( *, '(a,i6)' ) ' IERROR = ', ierror return end if nlo = 0 do i = 1, ntri do j = 1, 3 k = tnbr(j,i) if ( i < k ) then jp1 = j + 1 if ( 3 < jp1 ) then jp1 = 1 end if jp2 = jp1 + 1 if ( 3 < jp2 ) then jp2 = 1 end if a = til(j,i) b = til(jp1,i) c = til(jp2,i) if ( til(1,k) == b ) then d = til(3,k) else if ( til(2,k) == b ) then d = til(1,k) else d = til(2,k) end if if ( diaedg(vcl(1,c),vcl(2,c),vcl(1,a),vcl(2,a),vcl(1,d), & vcl(2,d),vcl(1,b),vcl(2,b)) == 1 ) then nlo = nlo + 1 end if end if end do end do write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) ' NLO = ', nlo call delaunay_print ( npt, vcl, ntri, til, tnbr ) end do return end subroutine test08 ( filename ) !*****************************************************************************80 ! !! TEST08 tests EQDIS2, INTPG, MFDEC2, MMASEP, SEPMDF, SEPSHP, SFDWMF, SFUPMF, TRISIZ. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 12 December 2025 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! This version by John Burkardt. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) integer, parameter :: incr = 10000 integer, parameter :: inunit = 1 integer, parameter :: maxed = 101 integer, parameter :: maxhv = 350 integer, parameter :: maxiw = 900 integer, parameter :: maxnc = 30 integer, parameter :: maxnv = 400 integer, parameter :: maxpv = 2000 integer, parameter :: maxvc = 800 integer, parameter :: maxwk = 1500 real ( kind = rk8 ) angmin real ( kind = rk8 ) angspc real ( kind = rk8 ) angtol real ( kind = rk8 ) area(maxhv) integer case real ( kind = rk8 ) degrees_to_radians real ( kind = rk8 ) dmin integer edge(4,maxed) character ( len = * ) filename real ( kind = rk8 ) h(maxhv) logical hflag integer ht(0:maxed-1) integer htsiz integer hvl(maxhv) integer i real ( kind = rk8 ) iang(maxpv) integer icur(maxnc) integer ierror integer ivrt(maxnv) integer iwk(maxiw) integer j real ( kind = rk8 ) kappa integer map(maxnc) integer msglvl integer ncur integer nh integer nhola integer nhole integer nmin integer npolg integer nrfv integer nsc integer ntri(maxhv) integer ntrid integer ntrie integer nv integer nvbc(maxnc) integer nvc integer nvcin integer nvert real ( kind = rk8 ), parameter :: pi = 3.141592653589793D+00 integer prime real ( kind = rk8 ) psi(maxhv) integer pvl(4,maxpv) real ( kind = rk8 ) radians_to_degrees integer regnum(maxhv) character ( len = 20 ) rgname real ( kind = rk8 ) tol real ( kind = rk8 ) tolin real ( kind = rk8 ) umdf2 real ( kind = rk8 ) vcl(2,maxvc) real ( kind = rk8 ) wk(maxwk) ! external umdf2 ! tol = 100.0D+00 * epsilon ( tol ) ! ! Read in vertices of general polygonal region. ! CASE = 1 : simple polygon or multiply connected polygonal region ! CASE = 2 : general polygonal region with holes and interfaces ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST08' write ( *, '(a)' ) ' Reading input file:' write ( *, '(a)' ) filename open ( unit = inunit, file = filename, form = 'formatted' ) read ( inunit, '(a)' ) rgname read ( inunit, * ) tolin,angspc,angtol,kappa,dmin,nmin,ntrid write ( *, '(a)' ) ' ' write ( *, '(a,a)' ) ' Region: ', rgname write ( *, 700) tol,angspc,angtol,kappa,dmin,nmin,ntrid angspc = degrees_to_radians ( angspc ) angtol = degrees_to_radians ( angtol ) hflag = ( 0.0D+00 <= kappa .and. kappa <= 1.0D+00 ) read ( inunit, * ) case,nvc,ncur,msglvl if ( maxvc < nvc ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST08 - Error' write ( *, '(a)' ) ' MAXVC < NVC.' return end if if ( maxnc < ncur ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST08 - Error' write ( *, '(a)' ) ' MAXNC < NCUR.' return end if read ( inunit, * ) nvbc(1:ncur) if ( case == 2 ) then read ( inunit, * ) icur(1:ncur) end if read ( inunit, * ) (vcl(1,i),vcl(2,i),i=1,nvc) ! ! Call DSMCPR or DSPGDC to set data structures in arrays ! REGNUM, HVL, PVL, IANG, HOLV = IWK. ! if ( case == 1 ) then nhole = ncur - 1 call dsmcpr ( nhole, nvbc, vcl, maxhv, maxpv, maxiw, nvc, npolg, nvert, & nhola, regnum, hvl, pvl, iang, iwk, ierror ) else if ( case == 2 ) then nv = sum ( nvbc(1:ncur) ) if ( maxnv < nv ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST08 - Error' write ( *, '(a)' ) ' MAXNV < NV.' return end if read ( inunit, * ) ivrt(1:nv) nsc = 0 do i = 1, nv if ( ivrt(i) < 0 ) then nsc = nsc + 1 end if end do nsc = nsc / 2 if ( maxed < nsc ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST08 - Fatal error' write ( *, '(a)' ) ' MAXED < NSC.' return end if htsiz = min ( prime ( nsc / 2 ), maxed ) call dspgdc ( nvc, vcl, incr, ncur, nvbc, icur, ivrt, maxhv, maxpv, maxiw, & npolg, nvert, nhole, nhola, regnum, hvl, pvl, iang, iwk, htsiz, nsc, & ht, edge, map, ierror ) end if close ( unit = inunit ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST08 - Error!' write ( *, '(a,i6)' ) ' IERROR = ', ierror return end if ! ! Print out data structures and measurements after initialization. ! nvcin = nvc nh = nhole*2 + nhola write ( *,670) msglvl write ( *,630) nvc,(i,vcl(1,i),vcl(2,i),i=1,nvc) write ( *,640) npolg,nhole,(hvl(i),i=1,npolg+nhole) write ( *,650) (regnum(i),i=1,npolg) write ( *,660) nvert,(i,(pvl(j,i),j=1,4),iang(i),i=1,nvert) write ( *,640) nhola,nh,(iwk(i),i=1,nh) nrfv = 0 do i = 1, nvert if ( pi + tol < iang(i) ) then nrfv = nrfv + 1 end if end do angmin = minval ( iang(1:nvert) ) angmin = radians_to_degrees ( angmin ) write (*,710) nvc,npolg,nvert,nhole,nhola,nrfv,angmin ! ! Obtain simple and convex polygon decompositions, and print measurements. ! if ( msglvl == 2) then write ( *,670) end if call spdec2 ( angspc,angtol,nvc,npolg,nvert,nhole,nhola,maxvc,maxhv, & maxpv,maxiw-nh,maxwk,iwk,vcl,regnum,hvl,pvl,iang,iwk(nh+1),wk,ierror) call cvdec2 ( angspc, angtol, nvc, npolg, nvert, maxvc, maxhv, maxpv, maxiw, & maxwk, vcl, regnum, hvl, pvl, iang, iwk, wk, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST08 - Error!' write ( *, '(a,i6)' ) ' IERROR = ', ierror return end if angmin = minval ( iang(1:nvert) ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST08:' write ( *, '(a)' ) ' Before call to EQDIS2:' write ( *, '(a,i6)' ) ' NVC = ', nvc write ( *, '(a,i6)' ) ' NPOLG = ', npolg write ( *, '(a,i6)' ) ' NVERT = ', nvert write ( *, '(a,g14.6)' ) ' ANGMIN = ', angmin ! ! Obtain further convex polygon decomposition based on mesh ! distribution function, and triangle sizes for the polygons. ! write ( *, '(a)' ) 'DEBUG: Call EQDIS2' call eqdis2 ( hflag, umdf2, kappa, angspc, angtol, dmin, nmin, ntrid, nvc, & npolg, nvert, maxvc, maxhv, maxpv, maxiw, maxwk, vcl, regnum, hvl, pvl, & iang, area, psi, h, iwk, wk, ierror ) if ( msglvl == 2 ) then write ( *,670) 0, 0, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00 end if if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST08 - Error!' write ( *, '(a,i6)' ) ' IERROR = ', ierror return end if ! ! Print out data structures and measurements after decomposition. ! write ( *,630) nvc,(i,vcl(1,i),vcl(2,i),i=nvcin+1,nvc) write ( *,680) npolg, hvl(1:npolg) write ( *,650) regnum(1:npolg) write ( *,660) nvert,(i,(pvl(j,i),j=1,4),iang(i),i=1,nvert) angmin = minval ( iang(1:nvert) ) angmin = radians_to_degrees ( angmin ) ntrie = 0 do i = 1, npolg ntri(i) = int ( ntrid * psi(i) * area(i) + 0.5D+00 ) ntrie = ntrie + ntri(i) end do write ( *,690) (i,area(i),psi(i),h(i),ntri(i),i=1,npolg) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST08:' write ( *, '(a)' ) ' After call to EQDIS2:' write ( *, '(a,i6)' ) ' NVC = ', nvc write ( *, '(a,i6)' ) ' NPOLG = ', npolg write ( *, '(a,i6)' ) ' NVERT = ', nvert write ( *, '(a,i6)' ) ' NTRIE = ', ntrie write ( *, '(a,g14.6)' ) ' ANGMIN = ', angmin 630 format (/1x,i7/(1x,i7,2f15.7)) 640 format (/1x,2i7/(1x,10i7)) 650 format (/(1x,10i7)) 660 format (/1x,i7/(1x,5i7,f15.7)) 670 format (1x,2i7,4f15.7) 680 format (/1x,i7/(1x,10i7)) 690 format (/(1x,i7,3e15.7,i7)) 700 format ( 'input : tol=',d15.7,' angspc=',f9.3, & ' angtol=',f9.3/9x,'kappa=',f9.3,' dmin=',f9.3, & ' nmin=',i5,' ntrid=',i7) 710 format (1x,'initds: nvc=',i7,' npolg=',i7,' nvert=',i7, & ' nhole=',i7/9x,'nhola=',i7,' nrfv=',i7,' angmin=',f9.3 ) return end subroutine test09 ( ) !*****************************************************************************80 ! !! TEST09 tests LRLINE. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 12 December 2025 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! This version by John Burkardt. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) real ( kind = rk8 ) dv integer lr integer lrline real ( kind = rk8 ) xu real ( kind = rk8 ) xv1 real ( kind = rk8 ) xv2 real ( kind = rk8 ) yu real ( kind = rk8 ) yv1 real ( kind = rk8 ) yv2 xu = 0.0D+00 yu = 0.0D+00 xv1 = 0.0D+00 yv1 = -1.0D+00 xv2 = 1.0D+00 yv2 = 0.0D+00 dv = 0.0D+00 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST09' write ( *, '(a)' ) ' LRLINE determines if a point is to the right,' write ( *, '(a)' ) ' left, or on a directed line that is a directed' write ( *, '(a)' ) ' distance away from a directed line from one' write ( *, '(a)' ) ' point to another.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The directed base line goes from' write ( *, '(2g14.6)' ) xv1, yv1 write ( *, '(a)' ) ' to' write ( *, '(2g14.6)' ) xv2, yv2 write ( *, '(a)' ) ' ' write ( *, '(a,g14.6)' ) ' The directed line distance is ', dv write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The point to be located is' write ( *, '(2g14.6)' ) xu, yu lr = lrline ( xu, yu, xv1, yv1, xv2, yv2, dv ) write ( *, '(a)' ) ' ' if ( lr == -1 ) then write ( *, '(a)' ) ' The point is to the right of the line.' else if ( lr == 0 ) then write ( *, '(a)' ) ' The point is on the line.' else if ( lr == 1 ) then write ( *, '(a)' ) ' The point is to the left of the line.' end if return end subroutine test10 ( ) !*****************************************************************************80 ! !! TEST10 tests LUFAC, LUSOL. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 12 December 2025 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! This version by John Burkardt. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) integer, parameter :: nmax = 20 real ( kind = rk8 ) a(nmax,nmax) real ( kind = rk8 ) b(nmax) real ( kind = rk8 ) emax real ( kind = rk8 ) esum integer i integer ipvt(nmax) integer j integer n integer seed logical singlr real ( kind = rk8 ) t real ( kind = rk8 ) tol real ( kind = rk8 ) urand tol = 100.0D+00 * epsilon ( tol ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST10' write ( *, '(a)' ) ' LUFAC factors a linear system;' write ( *, '(a)' ) ' LUSOL solves a factored linear system.' n = 4 seed = 1952 b(1:n) = 0.0D+00 do j = 1, n do i = 1, n a(i,j) = urand(seed) * 2.0D+00 - 1.0D+00 b(i) = b(i) + a(i,j) end do end do if ( n <= 4 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'a, b' do i = 1, n write ( *, '(5f15.7)' ) (a(i,j),j=1,n),b(i) end do end if call lufac ( a, nmax, n, tol, ipvt, singlr ) if ( singlr ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The matrix is singular' return end if call lusol ( a, nmax, n, ipvt, b ) if ( n <= 4 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ipvt, lu' ipvt(n) = n do i = 1, n write ( *, '(i5,4f15.7)' ) ipvt(i), (a(i,j),j=1,n) end do end if emax = 0.0D+00 esum = 0.0D+00 do i = 1, n t = abs ( b(i) - 1.0D+00 ) emax = max ( emax, t ) esum = esum + t end do write ( *,630) (b(i),i=1,min(4,n)) write ( *,640) emax,esum 630 format (' x = ',4f15.7) 640 format (' emax,esum = ',2e15.7) return end subroutine test11 ( filename ) !*****************************************************************************80 ! !! TEST11 tests DSMDF2, MDF2, PRMDF2. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 12 December 2025 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! This version by John Burkardt. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) integer, parameter :: incr = 10000 integer, parameter :: inunit = 1 integer, parameter :: maxed = 101 integer, parameter :: maxhv = 200 integer, parameter :: maxiw = 900 integer, parameter :: maxnc = 30 integer, parameter :: maxpv = 1000 integer, parameter :: maxvc = 500 integer, parameter :: maxwk = 1500 real ( kind = rk8 ) angspc real ( kind = rk8 ) angtol real ( kind = rk8 ) area(maxhv) integer case real ( kind = rk8 ) degrees_to_radians integer edge(4,maxed) real ( kind = rk8 ) edgval(maxpv) character ( len = * ) filename integer ht(0:maxed-1) integer htsiz integer hvl(maxhv) integer i real ( kind = rk8 ) iang(maxpv) integer icur(maxnc) integer ierror integer ifv integer ivrt(maxpv) integer iwk(maxiw) integer j integer map(maxnc) real ( kind = rk8 ) mdf2 integer ncur integer nev(maxhv) integer nh integer nhola integer nhole integer npolg integer nsc integer nv integer nvbc(maxnc) integer nvc integer nvert integer prime integer pvl(4,maxpv) integer regnum(maxhv) character ( len = 20 ) rgname real ( kind = rk8 ) tol real ( kind = rk8 ) tolin real ( kind = rk8 ) val(maxhv) real ( kind = rk8 ) vcl(2,maxvc) real ( kind = rk8 ) vrtval(maxvc) real ( kind = rk8 ) widsq(maxhv) real ( kind = rk8 ) wk(maxwk) real ( kind = rk8 ) x integer xivrt(maxhv+1) real ( kind = rk8 ) y tol = 100.0D+00 * epsilon ( tol ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST11' write ( *, '(a)' ) ' Reading input file:' write ( *, '(a)' ) filename open ( unit = inunit, file = filename, form = 'formatted' ) ! ! Read in vertices of general polygonal region. ! CASE = 1 : simple polygon or multiply connected polygonal region ! CASE = 2 : general polygonal region with holes and interfaces ! read ( inunit, '(a)' ) rgname read ( inunit, * ) tolin,angspc,angtol angspc = degrees_to_radians ( angspc ) angtol = degrees_to_radians ( angtol ) read ( inunit, * ) case,nvc,ncur if ( maxvc < nvc ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST11 - Error' write ( *, '(a)' ) ' MAXVC < NVC.' return end if if ( maxnc < ncur ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST11 - Error' write ( *, '(a)' ) ' MAXNC < NCUR.' return end if read ( inunit, * ) (nvbc(i),i=1,ncur) if ( case == 2 ) then read ( inunit, * ) (icur(i),i=1,ncur) end if read ( inunit, * ) (vcl(1,i),vcl(2,i),i=1,nvc) ! ! Call DSMCPR or DSPGDC to set data structures in arrays ! REGNUM, HVL, PVL, IANG, HOLV = IWK. ! if ( case == 1 ) then nhole = ncur - 1 call dsmcpr(nhole,nvbc,vcl,maxhv,maxpv,maxiw,nvc,npolg,nvert, & nhola,regnum,hvl,pvl,iang,iwk, ierror ) else if ( case == 2 ) then nv = 0 do i = 1, ncur nv = nv + nvbc(i) end do if ( maxpv < nv ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST11 - Error' write ( *, '(a)' ) ' MAXPC < NV.' return end if read ( inunit, * ) ivrt(1:nv) nsc = 0 do i = 1, nv if ( ivrt(i) < 0 ) then nsc = nsc + 1 end if end do nsc = nsc / 2 if ( maxed < nsc ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST11 - Error' write ( *, '(a)' ) ' MAXED < NSC.' return end if htsiz = min ( prime(nsc/2), maxed ) call dspgdc(nvc,vcl,incr,ncur,nvbc,icur,ivrt,maxhv,maxpv,maxiw, & npolg,nvert,nhole,nhola,regnum,hvl,pvl,iang,iwk,htsiz,nsc, & ht, edge, map, ierror ) end if if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST11 - Error!' write ( *, '(a,i6)' ) ' IERROR = ', ierror return end if ! ! Obtain simple and convex polygon decompositions. ! nh = nhole*2 + nhola call spdec2(angspc,angtol,nvc,npolg,nvert,nhole,nhola,maxvc,maxhv, & maxpv,maxiw-nh,maxwk,iwk,vcl,regnum,hvl,pvl,iang,iwk(nh+1),wk, ierror ) call cvdec2 ( angspc, angtol, nvc, npolg, nvert, maxvc, maxhv, maxpv, & maxiw, maxwk, vcl, regnum, hvl, pvl, iang, iwk, wk, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST11 - Error!' write ( *, '(a,i6)' ) ' IERROR = ', ierror return end if ! ! Initialize data structure for heuristic mesh distribution function ! and evaluate function at centroid of each convex polygon. ! call dsmdf2 ( .true., nvc, npolg, maxwk, vcl, hvl, pvl, iang, & ivrt, xivrt, widsq, edgval, vrtval, area, wk, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST11 - Error!' write ( *, '(a,i6)' ) ' IERROR = ', ierror return end if do i = 1, npolg call prmdf2(i,widsq(i),ivrt,xivrt,edgval,vrtval,nev(i),ifv,iwk) if ( nev(i) == 0 ) then val(i) = widsq(i) else nv = xivrt(i+1) - xivrt(i) x = 0.0D+00 y = 0.0D+00 do j = xivrt(i), xivrt(i+1)-1 x = x + vcl(1,ivrt(j)) y = y + vcl(2,ivrt(j)) end do val(i) = 1.0D+00 / mdf2(x/nv,y/nv,widsq(i),nev(i),ifv,iwk,ivrt, & edgval,vrtval,vcl) end if end do close ( unit = inunit ) ! ! Print arrays from calls to 3 mdf routines. ! write ( *,630) npolg,(i,xivrt(i),area(i),widsq(i),val(i),nev(i),i=1,npolg) write ( *,640) nvert,nvc,(i,ivrt(i),edgval(i),vrtval(i),i=1,nvc) write ( *,650) (i,ivrt(i),edgval(i),i=nvc+1,nvert) 630 format (/1x,i7/(1x,2i7,3f15.7,i7)) 640 format (/1x,2i7/(1x,2i7,2f15.7)) 650 format (1x,2i7,f15.7) return end subroutine test12 ( ) !*********************************************************************** ! !! TEST12 tests PRIME. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 12 December 2025 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! This version by John Burkardt. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) integer i integer prime write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST12' write ( *, '(a)' ) ' PRIME finds the smallest prime bigger than' write ( *, '(a)' ) ' a given value.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' I, PRIME(I)' write ( *, '(a)' ) ' ' do i = 100, 500, 100 write ( *, '(2i10)' ) i, prime ( i ) end do return end subroutine test13 ( filename ) !*****************************************************************************80 ! !! TEST13 tests PTPOLG. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 12 December 2025 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! This version by John Burkardt. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) integer, parameter :: inunit = 1 integer, parameter :: maxn = 100 real ( kind = rk8 ) a real ( kind = rk8 ) b real ( kind = rk8 ) c integer dim real ( kind = rk8 ) dtol character ( len = * ) filename integer i integer inout integer j integer n real ( kind = rk8 ) nrml(3) integer pgind(0:maxn) real ( kind = rk8 ) pt(3) real ( kind = rk8 ) tol real ( kind = rk8 ) vcl(3,maxn) tol = 100.0D+00 * epsilon ( tol ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST13' write ( *, '(a)' ) ' Reading input file:' write ( *, '(a)' ) filename open ( unit = inunit, file = filename, form = 'formatted' ) dtol = 10.0D+00 * tol read ( inunit, * ) n,dim,a,b if ( n < 3 ) then return end if if ( maxn < n ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DRPTPG - Error' write ( *, '(a)' ) ' MAXN < N.' return end if write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) 'DIM = ', dim write ( *, '(a,g14.6)' ) ' A = ', a write ( *, '(a,g14.6)' ) ' B = ', b write ( *, '(a)' ) ' ' read ( inunit, * ) (vcl(1,i),vcl(2,i),i=1,n) pgind(0) = n do i = 1, n pgind(i) = i if ( dim == 3 ) then vcl(3,i) = a * vcl(1,i) + b * vcl(2,i) end if write ( *, '(i5,3f15.7)' ) i,(vcl(j,i),j=1,dim) end do if ( dim == 3 ) then c = sqrt ( 1.0D+00 + a**2 + b**2) nrml(1) = - a / c nrml(2) = - b / c nrml(3) = 1.0D+00 / c end if write ( *, '(a)' ) ' ' 20 continue read ( inunit, *, end=30 ) pt(1),pt(2) if ( dim == 3 ) then pt(3) = a * pt(1) + b * pt(2) end if call ptpolg ( dim, 3, n, 1, pgind, vcl, pt, nrml, dtol, inout ) write ( *,630) inout,(pt(i),i=1,dim) go to 20 30 continue close ( unit = inunit ) 630 format (1x,'inout=',i3,3x,'pt=',3f15.7) return end subroutine test14 ( ) !*****************************************************************************80 ! !! TEST14 tests ROTIAR. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 12 December 2025 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! This version by John Burkardt. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) integer, parameter :: maxn = 50 integer a(maxn) integer i integer n integer shift n = 10 shift = 3 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST14' write ( *, '(a)' ) ' ROTIAR "rotates" an array.' write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) ' Using N = ', n write ( *, '(a,i6)' ) ' SHIFT = ', shift do i = 1, n a(i) = i end do write ( *, '(10i6)' ) a(1:n) write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Shifted array:' write ( *, '(a)' ) ' ' call rotiar ( n, a, shift ) write ( *, '(10i6)' ) a(1:n) return end subroutine test15 ( filename ) !*****************************************************************************80 ! !! TEST15 tests DIAM2, SHRNK2, WIDTH2. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 12 December 2025 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! This version by John Burkardt. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) integer, parameter :: inunit = 1 integer, parameter :: maxn = 100 real ( kind = rk8 ) diamsq character ( len = * ) filename integer i integer i1 integer i2 integer iedge(0:maxn) integer ierror integer n integer nshr real ( kind = rk8 ) sdist(0:maxn-1) real ( kind = rk8 ) widsq real ( kind = rk8 ) xc(0:maxn) real ( kind = rk8 ) xs(0:maxn) real ( kind = rk8 ) yc(0:maxn) real ( kind = rk8 ) ys(0:maxn) ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST15' write ( *, '(a)' ) ' Reading data file: "' // trim ( filename ) // '".' open ( unit = inunit, file = filename, form = 'formatted' ) read ( inunit, * ) n if ( maxn < n ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST15 - Error!' write ( *, '(a,i6)' ) ' N = ', n write ( *, '(a,i6)' ) ' Maximum legal N is ', maxn return end if read ( inunit, * ) (xc(i),yc(i),sdist(i),i=0,n-1) close ( unit = inunit ) xc(n) = xc(0) yc(n) = yc(0) write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) ' N = ', n write ( *, '(2f15.7)' ) (xc(i),yc(i),i=0,n) xs(0) = 0.0D+00 ys(0) = 0.0D+00 write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Calling SHRNK2' call shrnk2 ( n, xc, yc, sdist, nshr, xs, ys, iedge, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST15 - Error!' write ( *, '(a,i6)' ) ' SHRNK2 returns IERROR = ', ierror return end if write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) ' NSHR = ', n write ( *, '(2f15.7)' ) (xs(i),ys(i),i=0,nshr) call diam2 ( n, xc(1), yc(1), i1, i2, diamsq, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST15 - Error!' write ( *, '(a,i6)' ) ' DIAM2 returns IERROR = ', ierror return end if write ( *, '(a)' ) ' ' write ( *, '(2i5,f18.7)' ) i1,i2,diamsq call width2 ( n, xc(1), yc(1), i1, i2, widsq, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST15 - Error!' write ( *, '(a,i6)' ) ' WIDTH2 returns IERROR = ', ierror return end if write ( *, '(a)' ) ' ' write ( *, '(2i5,f18.7)' ) i1,i2,widsq if ( nshr < 3 ) then return end if call diam2 ( nshr, xs(1), ys(1), i1, i2, diamsq, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST15 - Error!' write ( *, '(a,i6)' ) ' DIAM2 returns IERROR = ', ierror return end if write ( *, '(a)' ) ' ' write ( *, '(2i5,f18.7)' ) i1,i2,diamsq call width2 ( nshr, xs(1), ys(1), i1, i2, widsq, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST15 - Error!' write ( *, '(a,i6)' ) ' WIDTH2 returns IERROR = ', ierror return end if write ( *, '(a)' ) ' ' write ( *, '(2i5,f18.7)' ) i1,i2,widsq return end subroutine test16 ( ) !*****************************************************************************80 ! !! TEST16 tests DHPSRT, IHPSRT, RANDPT. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 12 December 2025 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! This version by John Burkardt. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) integer, parameter :: maxk = 4 integer, parameter :: maxn = 100 integer axis real ( kind = rk8 ) da(maxk,maxn) integer i integer ia(maxk,maxn) logical iflag integer j integer k integer map(maxn) integer n integer nptav real ( kind = rk8 ), dimension ( maxk ) :: scale = (/ & 1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00 /) integer seed real ( kind = rk8 ), dimension ( maxk ) :: trans = (/ & 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00 /) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST16' k = 2 n = 10 seed = 1952 axis = 0 nptav = 0 iflag = (k < 0) k = abs(k) if ( k < 1 .or. maxk < k ) then return end if if ( n < 1 .or. maxn < n ) then return end if do j = 1, n map(j) = j end do call randpt ( k, n, seed, axis, nptav, scale, trans, maxk, da ) if ( iflag ) then do j = 1, n do i = 1, k ia(i,j) = int ( n * da(i,j) ) end do end do do j = 1, n write ( *, '(5i5)' ) j, (ia(i,j),i=1,k) end do call ihpsrt ( k, n, maxk, ia, map ) write ( *, '(a)' ) ' ' do j = 1, n write ( *, '(5i5)' ) map(j), (ia(i,map(j)),i=1,k) end do else do j = 1, n write ( *, '(i5,4f15.7)' ) j, (da(i,j),i=1,k) end do call dhpsrt ( k, n, maxk, da, map ) write ( *,'(a)' ) ' ' do j = 1, n write ( *, '(i5,4f15.7)' ) map(j), (da(i,map(j)),i=1,k) end do end if return end subroutine test17 ( filename ) !*****************************************************************************80 ! !! TEST17 tests BEDGMV, CVDTRI, FNDTRI, INTTRI, LOP, MTREDG, ROTPG, TMERGE, TRINBR, TRIPR2, TRPOLG. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 12 December 2025 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! This version by John Burkardt. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) integer, parameter :: incr = 10000 integer, parameter :: inunit = 1 integer, parameter :: maxed = 1201 integer, parameter :: maxhv = 350 integer, parameter :: maxiw = 900 integer, parameter :: maxnc = 30 integer, parameter :: maxnv = 400 integer, parameter :: maxpv = 2000 integer, parameter :: maxti = 8000 integer, parameter :: maxvc = 5000 integer, parameter :: maxwk = 1500 integer, parameter :: nfreq = 12 integer a integer afreq(0:nfreq-1) real ( kind = rk8 ) anga real ( kind = rk8 ) angb real ( kind = rk8 ) angc real ( kind = rk8 ) angle real ( kind = rk8 ) angmax real ( kind = rk8 ) angmin real ( kind = rk8 ) angspc real ( kind = rk8 ) angtol real ( kind = rk8 ) area(maxhv) integer b integer c integer case real ( kind = rk8 ) degrees_to_radians real ( kind = rk8 ) delta real ( kind = rk8 ) dmin integer edge(4,maxed) character ( len = * ) filename real ( kind = rk8 ) h(maxhv) logical hflag integer ht(0:maxed-1) integer htsiz integer hvl(maxhv) integer i real ( kind = rk8 ) iang(maxpv) integer icur(maxnc) integer ierror integer ivrt(maxnv) integer iwk(maxiw) integer j real ( kind = rk8 ) kappa integer map(maxnc) integer msglvl integer ncur integer nh integer nhola integer nhole integer nmin integer npolg integer nrfv integer nsc integer ntri integer ntrid integer nv integer nvbc(maxnc) integer nvc integer nvcin integer nvert real ( kind = rk8 ), parameter :: pi = 3.141592653589793D+00 integer prime real ( kind = rk8 ) psi(maxhv) integer pvl(4,maxpv) real ( kind = rk8 ) radians_to_degrees integer regnum(maxhv) character ( len = 20 ) rgname integer til(3,maxti) integer tnbr(3,maxti) real ( kind = rk8 ) tol real ( kind = rk8 ) tolin integer tstart(maxhv) real ( kind = rk8 ), external :: umdf2 real ( kind = rk8 ) vcl(2,maxvc) integer vnum(maxpv) integer vstart(maxpv) real ( kind = rk8 ) wk(maxwk) ! tol = 100.0D+00 * epsilon ( tol ) ! ! Read in vertices of general polygonal region. ! CASE = 1 : simple polygon or multiply connected polygonal region ! CASE = 2 : general polygonal region with holes and interfaces ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST17' write ( *, '(a)' ) ' Reading input file: "' // trim ( filename ) // '".' open ( unit = inunit, file = filename, form = 'formatted' ) read ( inunit, '(a)' ) rgname read ( inunit, * ) tolin, angspc, angtol, kappa, dmin, nmin, ntrid write ( *, '(a)' ) ' ' write ( *, 710 ) rgname, tol, angspc, angtol, kappa, dmin, nmin, ntrid angspc = degrees_to_radians ( angspc ) angtol = degrees_to_radians ( angtol ) hflag = ( 0.0D+00 <= kappa .and. kappa <= 1.0D+00 ) read ( inunit, * ) case, nvc, ncur, msglvl if ( maxvc < nvc ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST17 - Error' write ( *, '(a)' ) ' MAXVC < NVC.' return end if if ( maxnc < ncur ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST17 - Error' write ( *, '(a)' ) ' MAXNC < NCUR.' return end if read ( inunit, * ) nvbc(1:ncur) if ( case == 2 ) then read ( inunit, * ) icur(1:ncur) end if read ( inunit, * ) (vcl(1,i),vcl(2,i),i=1,nvc) ! ! Call DSMCPR or DSPGDC to set the data structures in arrays ! REGNUM, HVL, PVL, IANG, HOLV = IWK. ! if ( case == 1 ) then nhole = ncur - 1 call dsmcpr ( nhole, nvbc, vcl, maxhv, maxpv, maxiw, nvc, npolg, nvert, & nhola, regnum, hvl, pvl, iang, iwk, ierror ) else if ( case == 2 ) then nv = sum ( nvbc(1:ncur) ) if ( maxnv < nv ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST17 - Error' write ( *, '(a)' ) ' MAXNV < NV.' return end if read ( inunit, * ) ivrt(1:nv) nsc = 0 do i = 1, nv if ( ivrt(i) < 0 ) then nsc = nsc + 1 end if end do nsc = nsc / 2 if ( maxed < nsc ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST17 - Error' write ( *, '(a)' ) ' MAXED < NSC.' return end if htsiz = min ( prime(nsc/2), maxed ) call dspgdc ( nvc, vcl, incr, ncur, nvbc, icur, ivrt, maxhv, maxpv, maxiw, & npolg, nvert, nhole, nhola, regnum, hvl, pvl, iang, iwk, htsiz, nsc, & ht, edge, map, ierror ) end if close ( unit = inunit ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST17 - Error!' write ( *, '(a,i6)' ) ' IERROR = ', ierror return end if ! ! Print out the data structures and measurements after initialization. ! nvcin = nvc nh = nhole * 2 + nhola write ( *, '(a,i6)' ) 'MSGLVL = ', msglvl write ( *,630) nvc,(i,vcl(1,i),vcl(2,i),i=1,nvc) write ( *,640) npolg,nhole,(hvl(i),i=1,npolg+nhole) write ( *,650) regnum(1:npolg) write ( *,660) nvert, (i,(pvl(j,i),j=1,4), iang(i), i = 1,nvert) write ( *,640) nhola, nh, iwk(1:nh) nrfv = 0 do i = 1, nvert if ( pi + tol < iang(i) ) then nrfv = nrfv + 1 end if end do angmin = minval ( iang(1:nvert) ) angmin = radians_to_degrees ( angmin ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST17: INITDS:' write ( *, '(a,i6)' ) ' NVC = ', nvc write ( *, '(a,i6)' ) ' NPOLG = ', npolg write ( *, '(a,i6)' ) ' NVERT = ', nvert write ( *, '(a,i6)' ) ' NHOLE = ', nhole write ( *, '(a,i6)' ) ' NHOLA = ', nhola write ( *, '(a,i6)' ) ' NRFV = ', nrfv write ( *, '(a,g14.6)' ) ' ANGMIN = ', angmin ! ! Obtain simple and convex polygon decompositions, and print measurements. ! if ( msglvl == 2 ) then write ( *,670) end if call spdec2 ( angspc, angtol, nvc, npolg, nvert, nhole, nhola, maxvc, maxhv, & maxpv, maxiw-nh, maxwk, iwk, vcl, regnum, hvl, pvl, iang, iwk(nh+1), wk, & ierror ) call cvdec2 ( angspc, angtol, nvc, npolg, nvert, maxvc, maxhv, maxpv, maxiw, & maxwk, vcl, regnum, hvl, pvl, iang, iwk, wk, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST17 - Error!' write ( *, '(a,i6)' ) ' IERROR = ', ierror return end if angmin = minval ( iang(1:nvert) ) angmin = radians_to_degrees ( angmin ) write (*,730) nvc, npolg, nvert, angmin ! ! Obtain further convex polygon decomposition based on mesh ! distribution function, and triangle sizes for the polygons. ! Then print measurements. ! call eqdis2 ( hflag, umdf2, kappa, angspc, angtol, dmin, nmin, ntrid, nvc, & npolg, nvert, maxvc, maxhv, maxpv, maxiw, maxwk, vcl, regnum, hvl, pvl, & iang, area, psi, h, iwk, wk, ierror ) if ( msglvl == 2 ) then write ( *,670) 0, 0, 0.0D+00, 0.0D+00, 0.0D+00, 0.0D+00 end if if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST17 - Error!' write ( *, '(a,i6)' ) ' IERROR = ', ierror return end if write ( *,630) nvc,(i,vcl(1,i),vcl(2,i),i=nvcin+1,nvc) write ( *,680) npolg, hvl(1:npolg) write ( *,650) regnum(1:npolg) angmin = radians_to_degrees ( minval ( iang(1:nvert) ) ) write (*,740) nvc, npolg, nvert, angmin htsiz = min ( prime(nvcin*2), maxed ) nvcin = nvc ! ! Triangulate each convex polygon in the decomposition, according to ! the mesh spacings in the H array. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Calling TRIPR2' call tripr2 ( nvc, npolg, nvert, maxvc, maxti, maxiw, maxwk, h, vcl, hvl, & pvl, iang, ntri, til, vstart, vnum, tstart, iwk, wk, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST17 - Error!' write ( *, '(a,i6)' ) ' IERROR = ', ierror return end if ! ! Compute TNBR array. Then print arrays and measurements. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Calling TRINBR' call trinbr ( nvc, ntri, til, tnbr, htsiz, maxed, ht, edge, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST17 - Error!' write ( *, '(a,i6)' ) ' IERROR = ', ierror return end if write ( *,690) nvert,(i,(pvl(j,i),j=1,4),iang(i),vstart(i),vnum(i),i=1,nvert) write ( *,630) nvc,(i,vcl(1,i),vcl(2,i),i=nvcin+1,nvc) write ( *,650) tstart(1:npolg) write ( *,700) ntri,(i,(til(j,i),j=1,3),(tnbr(j,i),j=1,3),i=1,ntri) angmin = pi angmax = 0.0D+00 delta = pi / real ( nfreq, kind = rk8 ) afreq(0:nfreq-1) = 0 do i = 1, ntri a = til(1,i) b = til(2,i) c = til(3,i) anga = angle(vcl(1,c),vcl(2,c),vcl(1,a),vcl(2,a),vcl(1,b),vcl(2,b)) angb = angle(vcl(1,a),vcl(2,a),vcl(1,b),vcl(2,b),vcl(1,c),vcl(2,c)) angc = pi - anga - angb angmin = min ( angmin, anga, angb, angc ) angmax = max ( angmax, anga, angb, angc ) if ( abs ( anga - 0.5D+00 * pi ) <= tol ) then anga = 0.5D+00 * pi - 2.0D+00 * tol end if j = min ( nfreq - 1, int ( ( anga + tol ) / delta ) ) afreq(j) = afreq(j) + 1 if ( abs ( angb - 0.5D+00 * pi ) <= tol ) then angb = 0.5D+00 * pi - 2.0D+00 * tol end if j = min ( nfreq - 1, int ( ( angb + tol ) / delta ) ) afreq(j) = afreq(j) + 1 if ( abs ( angc - 0.5D+00 * pi ) <= tol ) then angc = 0.5D+00 * pi - 2.0D+00 * tol end if j = min ( nfreq - 1, int ( ( angc + tol ) / delta ) ) ! ! Having some odd cases where J comes out NEGATIVE. ! if ( j < 0 ) then j = 0 end if afreq(j) = afreq(j) + 1 end do angmin = radians_to_degrees ( angmin ) angmax = radians_to_degrees ( angmax ) write (*,750) nvc,ntri,angmin,angmax write (*,760) ( i*180.0D+00 / real ( nfreq, kind = rk8 ), & afreq(i) / real ( 3*ntri, kind = rk8 ),i=0,nfreq-1) 630 format (/1x,i7/(1x,i7,2f15.7)) 640 format (/1x,2i7/(1x,10i7)) 650 format (/(1x,10i7)) 660 format (/1x,i7/(1x,5i7,f15.7)) 670 format (1x,2i7,4f15.7) 680 format (/1x,i7/(1x,10i7)) 690 format (/1x,i7/(1x,5i7,f15.7,2i7)) 700 format (/1x,i7/(1x,7i7)) 710 format (1x,a20/1x,'input : tol=',d15.7,' angspc=',f9.3, & ' angtol=',f9.3/9x,'kappa=',f9.3,' dmin=',f9.3, & ' nmin=',i5,' ntrid=',i7) 730 format (1x,'decomp: nvc=',i7,' npolg=',i7,' nvert=',i7/ & 9x,'angmin=',f9.3 ) 740 format (1x,'eqdist: nvc=',i7,' npolg=',i7,' nvert=',i7/ & 9x,'angmin=',f9.3 ) 750 format (1x,'triang: nvc=',i7,' ntri=',i7 / & 9x,'angmin=',f9.3,' angmax=',f9.3 ) 760 format (1x,'relative frequency of triangle angles'/4(f8.1,f10.5)) return end subroutine test18 ( ) !*****************************************************************************80 ! !! TEST18 tests VISVRT, VORNBR. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 12 December 2025 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! This version by John Burkardt. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) integer, parameter :: maxn = 200 real ( kind = rk8 ) angle real ( kind = rk8 ) angspc real ( kind = rk8 ) degrees_to_radians integer i integer ierror integer ivert integer ivis(0:maxn) integer ivor(0:maxn) integer j integer maxnv integer n integer nvis integer nvor integer nvrt integer nvsvrt real ( kind = rk8 ) phi real ( kind = rk8 ) theta(0:maxn) real ( kind = rk8 ) x(maxn) real ( kind = rk8 ) xc(0:maxn) real ( kind = rk8 ) xeye real ( kind = rk8 ) xvor(0:maxn) real ( kind = rk8 ) y(maxn) real ( kind = rk8 ) yc(0:maxn) real ( kind = rk8 ) yeye real ( kind = rk8 ) yvor(0:maxn) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST18' write ( *, '(a)' ) ' Test being SKIPPED FOR NOW.' return read ( *, * ) n if ( n < 3 ) then return end if if ( maxn < n ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST18 - Error' write ( *, '(a)' ) ' MAXN < N.' return end if read ( *, * ) ( x(i), y(i), i = 1, n ) read ( *, * ) ivert read ( *, * ) angspc ! ! 1 <= IVERT <= N is index of polygon vertex for eyepoint. ! ANGSPC is angle spacing parameter in degrees. ! angspc = degrees_to_radians ( angspc ) xeye = x(ivert) yeye = y(ivert) nvrt = n - 2 j = -1 do i = ivert+1, n j = j + 1 xc(j) = x(i) yc(j) = y(i) end do do i = 1, ivert-1 j = j + 1 xc(j) = x(i) yc(j) = y(i) end do write ( *, '(a,2g14.6)' ) '(XEYE,YEYE)= ', xeye, yeye write ( *,630) nvrt,(xc(i),yc(i),i=0,nvrt) call vispol ( xeye, yeye, nvrt, xc, yc, nvis, ivis ) write ( *,640) nvis,(xc(i),yc(i),ivis(i),i=0,nvis) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST18 - Fatal error!' write ( *, '(a,i6)' ) ' IERROR = ', ierror return end if phi = angle ( xc(nvis), yc(nvis), xeye, yeye, xc(0), yc(0) ) maxnv = nvis + int ( phi / angspc ) if ( maxn < maxnv ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST18 - Fatal error!' write ( *, '(a)' ) ' MAXN < MAXNV.' return end if call visvrt ( angspc, xeye, yeye, nvis, xc, yc, ivis, maxnv, nvsvrt, theta ) write ( *,650) nvsvrt,(xc(i),yc(i),ivis(i),theta(i),i=0,nvsvrt) nvor = -1 call vornbr ( xeye, yeye, nvsvrt, xc, yc, nvor, ivor, xvor, yvor, ierror ) write ( *,640) nvor,(xvor(i),yvor(i),ivor(i),i=0,nvor) 630 format (/1x,i5/(1x,2f15.7)) 640 format (/1x,i5/(1x,2f15.7,i5)) 650 format (/1x,i5/(1x,2f15.7,i5,f15.7)) return end subroutine test19 ( ) !*****************************************************************************80 ! !! TEST19 tests ROTIPG, VISPOL. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 12 December 2025 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! This version by John Burkardt. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) integer, parameter :: maxn = 200 integer i integer ivert integer ivis(0:maxn+1) integer j integer n integer nvis integer nvrt integer vptype real ( kind = rk8 ) x(maxn) real ( kind = rk8 ) xc(0:maxn+1) real ( kind = rk8 ) xeye real ( kind = rk8 ) y(maxn) real ( kind = rk8 ) yc(0:maxn+1) real ( kind = rk8 ) yeye write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST19' write ( *, '(a)' ) ' TEST BEING SKIPPED FOR NOW.' return read ( *, * ) n if ( n < 3 ) then return end if if ( maxn < n ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST19 - Error' write ( *, '(a)' ) ' MAXN < N.' return end if read ( *,*) (x(i),y(i),i=1,n) read ( *,*) ivert,vptype,xeye,yeye ! ! 1 <= IVERT <= N is index of polygon vertex or IVERT <= 0 if ! ROTIPG is to be called. BVTYPE = 0, 1, or 2 for boundary, ! interior, or blocked exterior viewpoint. (XEYE,YEYE) is needed ! for non-boundary viewpoints only, and must be visible from ! (X(IVERT),Y(IVERT)) if 0 < IVERT. ! if ( vptype == 0 ) then xeye = x(ivert) yeye = y(ivert) nvrt = n - 2 j = -1 do i = ivert+1, n j = j + 1 xc(j) = x(i) yc(j) = y(i) end do do i = 1, ivert-1 j = j + 1 xc(j) = x(i) yc(j) = y(i) end do else if ( vptype == 1 ) then nvrt = n if ( ivert <= 0 ) then do i = 1, n xc(i) = x(i) yc(i) = y(i) end do xc(0) = xc(n) yc(0) = yc(n) call rotipg ( xeye, yeye, nvrt, xc, yc ) else j = -1 do i = ivert, n j = j + 1 xc(j) = x(i) yc(j) = y(i) end do do i = 1, ivert j = j + 1 xc(j) = x(i) yc(j) = y(i) end do end if else nvrt = n if ( ivert <= 0 ) then do i = 1, n xc(n-i) = x(i) yc(n-i) = y(i) end do xc(n) = xc(0) yc(n) = yc(0) call rotipg ( xeye, yeye, nvrt, xc, yc ) else j = -1 do i = ivert, 1, -1 j = j + 1 xc(j) = x(i) yc(j) = y(i) end do do i = n, ivert, -1 j = j + 1 xc(j) = x(i) yc(j) = y(i) end do end if end if write ( *, '(a)' ) ' ' write ( *, '(a,2g14.6)' ) '(XEYE,YEYE)=', xeye,yeye write ( *,630) nvrt,(xc(i),yc(i),i=0,nvrt) call vispol ( xeye, yeye, nvrt, xc, yc, nvis, ivis ) write ( *,640) nvis,(xc(i),yc(i),ivis(i),i=0,nvis) 630 format (/1x,i5/(1x,2f15.7)) 640 format (/1x,i5/(1x,2f15.7,i5)) return end subroutine pcdec ( ) !*****************************************************************************80 ! !! PCDEC generates PICTEX commands from DRDEC or DREQD output. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 12 December 2025 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! This version by John Burkardt. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) integer, parameter :: edgv = 4 integer, parameter :: loc = 1 integer, parameter :: maxho = 50 integer, parameter :: maxnc = 30 integer, parameter :: maxnv = 500 integer, parameter :: succ = 3 integer holv(maxho) integer hvl(maxnc*2) integer i real ( kind = rk8 ) iang(maxnv*2) integer j integer j1 integer l integer l1 logical mid integer msglvl integer nh integer nhola integer nhole integer npolg integer nvc integer nvert integer pvl(4,maxnv*2) integer regnum(maxnc*2) integer t real ( kind = rk8 ) vcl(2,maxnv) real ( kind = rk8 ) x1 real ( kind = rk8 ) x2 real ( kind = rk8 ) xmax real ( kind = rk8 ) xmin real ( kind = rk8 ) y1 real ( kind = rk8 ) y2 real ( kind = rk8 ) ymax real ( kind = rk8 ) ymin read ( *,*) msglvl if ( msglvl /= 2) then return end if read ( *,*) nvc read ( *,*) (t,vcl(1,i),vcl(2,i),i=1,nvc) read ( *,*) npolg,nhole read ( *,*) (hvl(i),i=1,npolg+nhole) read ( *,*) (regnum(i),i=1,npolg) read ( *,*) nvert read ( *,*) (t,(pvl(j,i),j=1,4),iang(i),i=1,nvert) read ( *,*) nhola,nh read ( *,*) (holv(i),i=1,nh) xmin = vcl(1,1) ymin = vcl(2,1) xmax = xmin ymax = ymin do i = 2, nvc xmin = min(xmin,vcl(1,i)) xmax = max(xmax,vcl(1,i)) ymin = min(ymin,vcl(2,i)) ymax = max(ymax,vcl(2,i)) end do write ( *, '(a)' ) '\\begin{figure}[htbp]' write ( *, '(a)' ) '\\begin{center}' write ( *, '(a)' ) '\\mbox{\\beginpicture' write ( *, '(a)' ) '\\setcoordinatesystem units <0.2cm,0.2cm>' write ( *,610) xmin,xmax,ymin,ymax write ( *, '(a)' ) '\\setlinear' do i = 1, npolg+nhole mid = .false. j = hvl(i) l = pvl(loc,j) 20 continue j1 = pvl(succ,j) l1 = pvl(loc,j1) if ( pvl(edgv,j) < j ) then if ( .not. mid ) then write ( *, '(a)' ) '\\plot' mid = .true. write ( *, '(2f15.7)' ) vcl(1,l),vcl(2,l) end if write ( *, '(2f15.7)' ) vcl(1,l1),vcl(2,l1) else if ( mid ) then write ( *, '(a)' ) '/' mid = .false. end if j = j1 l = l1 if ( j /= hvl(i) ) then go to 20 end if if ( mid ) then write ( *, '(a)' ) '/' end if end do 40 continue read ( *,*) i,j,x1,y1,x2,y2 if ( i == 0 ) then go to 50 end if write ( *,630) x1,y1,x2,y2 go to 40 50 continue write ( *, '(a)' ) '\\end{picture}' write ( *, '(a)' ) '\\end{center}' write ( *, '(a)' ) '\\end{figure}' 610 format ('\\setplotarea x from ',f9.4,' to ',f9.4,', y from ',f9.4, & ' to ',f9.4) 630 format ('\\plot ',4f15.7,' /') return end subroutine pcds ( ) !*****************************************************************************80 ! !! PCDS generates PICTEX commands from DRDS output. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 12 December 2025 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! This version by John Burkardt. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) integer, parameter :: edgv = 4 integer, parameter :: loc = 1 integer, parameter :: maxho = 50 integer, parameter :: maxnc = 30 integer, parameter :: maxnv = 500 integer, parameter :: succ = 3 integer holv(maxho) integer hvl(maxnc*2) integer i real ( kind = rk8 ) iang(maxnv*2) integer j integer j1 integer l integer l1 logical mid integer nh integer nhola integer nhole integer npolg integer nvc integer nvert integer pvl(4,maxnv*2) integer regnum(maxnc*2) integer t real ( kind = rk8 ) vcl(2,maxnv) real ( kind = rk8 ) xmax real ( kind = rk8 ) xmin real ( kind = rk8 ) ymax real ( kind = rk8 ) ymin ! read ( *,*) nvc read ( *,*) (t,vcl(1,i),vcl(2,i),i=1,nvc) read ( *,*) npolg,nhole read ( *,*) (hvl(i),i=1,npolg+nhole) read ( *,*) (regnum(i),i=1,npolg) read ( *,*) nvert read ( *,*) (t,(pvl(j,i),j=1,4),iang(i),i=1,nvert) read ( *,*) nhola,nh read ( *,*) (holv(i),i=1,nh) xmin = vcl(1,1) ymin = vcl(2,1) xmax = xmin ymax = ymin do i = 2, nvc xmin = min(xmin,vcl(1,i)) xmax = max(xmax,vcl(1,i)) ymin = min(ymin,vcl(2,i)) ymax = max(ymax,vcl(2,i)) end do write ( *, '(a)' ) '\\begin{figure}[htbp]' write ( *, '(a)' ) '\\begin{center}' write ( *, '(a)' ) '\\mbox{\\beginpicture' write ( *, '(a)' ) '\\setcoordinatesystem units <0.2cm,0.2cm>' write ( *, 610) xmin,xmax,ymin,ymax write ( *, '(a)' ) '\\setlinear' do i = 1, npolg+nhole mid = .false. j = hvl(i) l = pvl(loc,j) 20 continue j1 = pvl(succ,j) l1 = pvl(loc,j1) if ( pvl(edgv,j) < j ) then if ( .not. mid ) then write ( *, '(a)' ) '\\plot' mid = .true. write ( *, '(2f15.7)' ) vcl(1,l),vcl(2,l) end if write ( *, '(2f15.7)' ) vcl(1,l1),vcl(2,l1) else if ( mid ) then write ( *, '(a)' ) '/' mid = .false. end if j = j1 l = l1 if ( j /= hvl(i) ) then go to 20 end if if ( mid ) then write ( *, '(a)' ) '/' end if end do write ( *, '(a)' ) '\\endpicture}' write ( *, '(a)' ) '\\end{center}' write ( *, '(a)' ) '\\end{figure}' 610 format ('\\setplotarea x from ',f9.4,' to ',f9.4,', y from ',f9.4, & ' to ',f9.4) return end subroutine pctri ( ) !*****************************************************************************80 ! !! PCTRI generates PICTEX commands from DRTRI output. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 12 December 2025 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! This version by John Burkardt. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) integer, parameter :: maxti = 8000 integer, parameter :: maxvc = 5000 integer a integer b integer i integer j integer jp1 integer msglvl integer nh integer nhola integer nhole integer npolg integer ntri integer nvc integer nvcin integer nvert integer t integer til(3,maxti) integer tnbr(3,maxti) real ( kind = rk8 ) vcl(2,maxvc) real ( kind = rk8 ) x real ( kind = rk8 ) xmax real ( kind = rk8 ) xmin real ( kind = rk8 ) ymax real ( kind = rk8 ) ymin ! read ( *,*,end=40) xmin,xmax,ymin,ymax read ( *,*) msglvl if ( msglvl /= 0) then return end if read ( *,*) nvc read ( *,*) (t,vcl(1,i),vcl(2,i),i=1,nvc) read ( *,*) npolg,nhole read ( *,*) (t,i=1,npolg+nhole) read ( *,*) (t,i=1,npolg) read ( *,*) nvert read ( *,*) ((t,j=1,5),x,i=1,nvert) read ( *,*) nhola,nh read ( *,*) (t,i=1,nh) nvcin = nvc read ( *,*) nvc read ( *,*) (t,vcl(1,i),vcl(2,i),i=nvcin+1,nvc) read ( *,*) npolg read ( *,*) (t,i=1,npolg) read ( *,*) (t,i=1,npolg) read ( *,*) nvert read ( *,*) ((t,j=1,5),x,t,t,i=1,nvert) nvcin = nvc read ( *,*) nvc read ( *,*) (t,vcl(1,i),vcl(2,i),i=nvcin+1,nvc) read ( *,*) (t,i=1,npolg) read ( *,*) ntri read ( *,*) (t,(til(j,i),j=1,3),(tnbr(j,i),j=1,3),i=1,ntri) ! write ( *, '(a)' ) '\\begin{figure}[htbp]' write ( *, '(a)' ) '\\begin{center}' write ( *, '(a)' ) '\\mbox{\\beginpicture' write ( *, '(a)' ) '\\setcoordinatesystem units <1.0cm,1.0cm>' write ( *,610) xmin,xmax,ymin,ymax write ( *, '(a)' ) '\\setlinear' do i = 1, ntri do j = 1, 3 if ( tnbr(j,i) < i ) then if ( j <= 2 ) then jp1 = j + 1 else jp1 = 1 end if a = til(j,i) b = til(jp1,i) if ( xmin <= vcl(1,a) .and. vcl(1,a) <= xmax .and. & ymin <= vcl(2,a) .and. vcl(2,a) <= ymax .and. & xmin <= vcl(1,b) .and. vcl(1,b) <= xmax .and. & ymin <= vcl(2,b) .and. vcl(2,b) <= ymax) then write ( *,620) vcl(1,a),vcl(2,a),vcl(1,b),vcl(2,b) end if end if end do end do write ( *, '(a)' ) '\\endpicture}' write ( *, '(a)' ) '\\end{center}' write ( *, '(a)' ) '\\end{figure}' 40 continue 610 format ('\\setplotarea x from ',f9.4,' to ',f9.4,', y from ',f9.4, & ' to ',f9.4) 620 format ('\\plot ',4f15.7,' /') return end subroutine pitri ( ) !*****************************************************************************80 ! !! PITRI generates PIC commands from DRTRI output. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 12 December 2025 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! This version by John Burkardt. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) integer, parameter :: maxti = 8000 integer, parameter :: maxvc = 5000 integer a integer b integer i integer j integer jp1 integer msglvl integer nh integer nhola integer nhole integer npolg integer ntri integer nvc integer nvcin integer nvert integer t integer til(3,maxti) integer tnbr(3,maxti) real ( kind = rk8 ) vcl(2,maxvc) real ( kind = rk8 ) x real ( kind = rk8 ) xmax real ( kind = rk8 ) xmin real ( kind = rk8 ) ymax real ( kind = rk8 ) ymin read ( *,*,end=40) xmin,xmax,ymin,ymax read ( *,*) msglvl if ( msglvl /= 0 ) then return end if read ( *, * ) nvc read ( *, * ) (t,vcl(1,i),vcl(2,i),i=1,nvc) read ( *, * ) npolg,nhole read ( *, * ) (t,i=1,npolg+nhole) read ( *, * ) (t,i=1,npolg) read ( *, * ) nvert read ( *, * ) ((t,j=1,5),x,i=1,nvert) read ( *, * ) nhola,nh read ( *, * ) (t,i=1,nh) nvcin = nvc read ( *, * ) nvc read ( *, * ) (t,vcl(1,i),vcl(2,i),i=nvcin+1,nvc) read ( *, * ) npolg read ( *, * ) (t,i=1,npolg) read ( *, * ) (t,i=1,npolg) read ( *, * ) nvert read ( *, * ) ((t,j=1,5),x,t,t,i=1,nvert) nvcin = nvc read ( *, * ) nvc read ( *, * ) (t,vcl(1,i),vcl(2,i),i=nvcin+1,nvc) read ( *, * ) (t,i=1,npolg) read ( *, * ) ntri read ( *, * ) (t,(til(j,i),j=1,3),(tnbr(j,i),j=1,3),i=1,ntri) write ( *, '(a)' ) '.ps 5.5i' do i = 1, ntri do j = 1, 3 if ( tnbr(j,i) < i ) then if ( j <= 2 ) then jp1 = j + 1 else jp1 = 1 end if a = til(j,i) b = til(jp1,i) if ( xmin <= vcl(1,a) .and. vcl(1,a) <= xmax .and. & ymin <= vcl(2,a) .and. vcl(2,a) <= ymax .and. & xmin <= vcl(1,b) .and. vcl(1,b) <= xmax .and. & ymin <= vcl(2,b) .and. vcl(2,b) <= ymax) then write ( *,620) -vcl(2,a),vcl(1,a),-vcl(2,b),vcl(1,b) end if ! ! Rotate by 90 degrees in above line. ! end if end do end do write ( *, '(a)' ) '.pe' 40 continue 620 format ('line from ',f11.7,', ',f11.7,' to ', f11.7,', ',f11.7) return end subroutine test20 ( ) !*****************************************************************************80 ! !! TEST20 tests XEDGE. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 12 December 2025 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! This version by John Burkardt. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) logical intsct integer mode real ( kind = rk8 ) xu real ( kind = rk8 ) xv1 real ( kind = rk8 ) xv2 real ( kind = rk8 ) xw1 real ( kind = rk8 ) xw2 real ( kind = rk8 ) yu real ( kind = rk8 ) yv1 real ( kind = rk8 ) yv2 real ( kind = rk8 ) yw1 real ( kind = rk8 ) yw2 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST20' write ( *, '(a)' ) ' XEDGE determines whether two edges,' write ( *, '(a)' ) ' or an edge and a ray, intersect.' write ( *, '(a)' ) ' (An edge is a finite line segment.)' write ( *, '(a)' ) ' (A ray is a semi-infinite line segment.)' mode = 0 xv1 = 3.0D+00 yv1 = 0.0D+00 xv2 = 3.0D+00 yv2 = 2.0D+00 xw1 = 0.0D+00 yw1 = 0.0D+00 xw2 = 6.0D+00 yw2 = 2.0D+00 write ( *, '(a)' ) ' ' if ( mode == 0 ) then write ( *, '(a)' ) ' Edge 1 is from' write ( *, * ) ' (', xv1, ',', yv1, ') to' write ( *, * ) ' (', xv2, ',', yv2, ').' else write ( *, '(a)' ) ' Ray 1 is from' write ( *, * ) ' (', xv1, ',', yv1, ') through' write ( *, * ) ' (', xv2, ',', yv2, ').' end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Edge 2 is from' write ( *, * ) ' (', xw1, ',', yw1, ') to' write ( *, * ) ' (', xw2, ',', yw2, ').' call xedge ( mode, xv1, yv1, xv2, yv2, xw1, yw1, xw2, yw2, xu, yu, intsct ) if ( .not. intsct ) then write ( *, '(a)' ) 'The objects either do not intersect' write ( *, '(a)' ) 'or are parallel, or coincide.' else write ( *, '(a,2g14.6)' ) 'The objects intersect at ', xu, yu end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) '(Expecting the answer (3,1) ).' return end subroutine test21 ( ) !*****************************************************************************80 ! !! TEST21 tests XLINE. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 12 December 2025 ! ! Author: ! ! Original FORTRAN77 version by Barry Joe. ! This version by John Burkardt. ! implicit none integer, parameter :: rk8 = kind ( 1.0D+00 ) real ( kind = rk8 ) dv real ( kind = rk8 ) dw logical parall real ( kind = rk8 ) xu real ( kind = rk8 ) xv1 real ( kind = rk8 ) xv2 real ( kind = rk8 ) xw1 real ( kind = rk8 ) xw2 real ( kind = rk8 ) yu real ( kind = rk8 ) yv1 real ( kind = rk8 ) yv2 real ( kind = rk8 ) yw1 real ( kind = rk8 ) yw2 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST21' write ( *, '(a)' ) ' XLINE finds the intersection of two lines.' write ( *, '(a)' ) ' Each line is defined as the line a given' write ( *, '(a)' ) ' distance to the left of a line through two' write ( *, '(a)' ) ' points.' xv1 = 0.0D+00 yv1 = 0.0D+00 xv2 = 0.0D+00 yv2 = 1.0D+00 dv = - 6.0D+00 xw1 = 0.0D+00 yw1 = 0.0D+00 xw2 = 3.0D+00 yw2 = 1.0D+00 dw = 0.0D+00 write ( *, '(a)' ) ' ' write ( *, * ) ' Line 1 is ', dv, ' units left of the line' write ( *, * ) ' through (', xv1, ',', yv1, ') and' write ( *, * ) ' (', xv2, ',', yv2, ').' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' ' write ( *, * ) ' Line 2 is ', dw, ' units left of the line' write ( *, * ) ' through (', xw1, ',', yw1, ') and' write ( *, * ) ' (', xw2, ',', yw2, ').' call xline ( xv1, yv1, xv2, yv2, xw1, yw1, xw2, yw2, dv, dw, xu, yu, parall ) write ( *, '(a)' ) ' ' if ( parall ) then write ( *, '(a)' ) ' The lines are parallel or coincide.' else write ( *, '(a,2g14.6)' ) ' The lines intersect at ', xu, yu end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' (Expecting the answer (6,2) ).' return end subroutine timestamp ( ) !*****************************************************************************80 ! !! timestamp() prints the current YMDHMS date as a time stamp. ! ! Example: ! ! 31 May 2001 9:45:54.872 AM ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 12 December 2025 ! ! Author: ! ! John Burkardt ! implicit none character ( len = 8 ) ampm integer d integer h integer m integer mm character ( len = 9 ), parameter, dimension(12) :: month = (/ & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' /) integer n integer s integer values(8) integer y call date_and_time ( values = values ) y = values(1) m = values(2) d = values(3) h = values(5) n = values(6) s = values(7) mm = values(8) if ( h < 12 ) then ampm = 'AM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h < 12 ) then ampm = 'PM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Midnight' else ampm = 'AM' end if end if end if write ( *, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm ) return end