program geompack_prb ! !*********************************************************************** ! !! GEOMPACK_PRB runs the GEOMPACK tests. ! implicit none ! call timestamp ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'GEOMPACK_PRB:' write ( *, '(a)' ) ' GEOMPACK tests.' 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 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'GEOMPACK_PRB:' write ( *, '(a)' ) ' Normal end of execution.' stop end subroutine test01 ! !*********************************************************************** ! !! TEST01 tests ANGLE; !! TEST01 tests AREAPG; !! TEST01 tests AREATR. ! implicit none ! integer, parameter :: n = 9 ! double precision ang(n) double precision angle double precision area double precision areapg double precision areatr double precision atr(n) integer i double precision sum2 double precision, 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 /) double precision, 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 ! !*********************************************************************** ! !! TEST02 tests CMCIRC. ! implicit none ! integer cmcirc integer in double precision, parameter :: x0 = 3.0D+00 double precision, parameter :: y0 = 3.0D+00 double precision, parameter :: x1 = 5.0D+00 double precision, parameter :: y1 = 0.0D+00 double precision, parameter :: x2 = 0.0D+00 double precision, parameter :: y2 = 5.0D+00 double precision, parameter :: x3 = -5.0D+00 double precision, 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 ) ! !*********************************************************************** ! !! TEST03 tests CVDEC2; !! TEST03 tests FNDSEP; !! TEST03 tests INSED2; !! TEST03 tests INSVR2; !! TEST03 tests JNHOLE; !! TEST03 tests MINANG; !! TEST03 tests RESVRT; !! TEST03 tests SPDEC2. ! implicit none ! 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 ! double precision angmin double precision angspc double precision angtol integer case double precision degrees_to_radians integer edge(4,maxed) character ( len = * ) filename integer ht(0:maxed-1) integer htsiz integer hvl(maxhv) integer i double precision iang(maxpv) integer icur(maxnc) integer ierror integer iprt 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 double precision d_pi integer prime integer pvl(4,maxpv) double precision radians_to_degrees integer regnum(maxhv) character ( len = 20 ) rgname double precision tol double precision tolin double precision vcl(2,maxvc) double precision wk(maxwk) ! tol = 100.0D+00 * epsilon ( tol ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST03' 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 ! CASE = -1 or -2: similar, but decompose into simple polygons only ! and don't obtain convex polygon decomposition ! read ( inunit, '(a)' ) rgname read ( inunit, * ) tolin, angspc, angtol write ( *, '(a)' ) ' ' write ( *, '(a)' ) rgname write ( *, 700 ) tol,angspc,angtol angspc = degrees_to_radians ( angspc ) angtol = degrees_to_radians ( angtol ) read ( inunit, * ) case, nvc, ncur, msglvl if ( nvc > maxvc ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DRDEC - Error!' write ( *, '(a)' ) ' NVC > MAXVC.' return end if if ( ncur > maxnc ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DRDEC - Error!' write ( *, '(a)' ) ' NCUR > MAXNC.' return end if read ( inunit, * ) ( nvbc(i), i = 1, ncur ) if ( abs(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 ( abs ( 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 ( abs(case) == 2 ) then nv = 0 do i = 1, ncur nv = nv + nvbc(i) end do if ( nv > maxnv ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DRDEC - Error!' write ( *, '(a)' ) ' NV > MAXNV.' return end if read ( inunit, * ) (ivrt(i),i=1,nv) nsc = 0 do i = 1, nv if ( ivrt(i) < 0) then nsc = nsc + 1 end if end do nsc = nsc / 2 if ( nsc > maxed ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DRDEC - Error!' write ( *, '(a)' ) ' NSC > MAXED.' 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)' ) 'DRDEC - 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 ( iang(i) > d_pi ( ) + tol ) 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 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 ) if ( case > 0 ) then 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)' ) 'DRDEC - 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)) 700 format ('input : tol=',d15.7,' angspc=',f9.3,' angtol=',f9.3) 710 format (1x,'initds: nvc=',i7,' npolg=',i7,' nvert=',i7, & ' nhole=',i7/9x,'nhola=',i7,' nrfv=',i7,' angmin=',f9.3 ) return end subroutine test04 ! !*********************************************************************** ! !! TEST04 tests DIAEDG. ! implicit none ! integer diaedg integer in double precision x0 double precision x1 double precision x2 double precision x3 double precision y0 double precision y1 double precision y2 double precision 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 ) ! !*********************************************************************** ! !! TEST05 tests DSMCPR; !! TEST05 tests DSPGDC; !! TEST05 tests EDGHT; !! TEST05 tests HOLVRT. ! implicit none ! integer, parameter :: incr = 10000 integer, parameter :: inunit = 1 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 double precision iang(maxnv*2) integer icur(maxnc) integer ierror 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 double precision tolin double precision vcl(2,maxnv) ! ! 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)' ) '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 ( nvc > maxnv ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST05 - Error!' write ( *, '(a)' ) ' NVC > MAXNV.' return end if if ( ncur > maxnc ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST05 - Error!' write ( *, '(a)' ) ' NCUR > MAXNC.' 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 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 nv = 0 do i = 1, ncur nv = nv + nvbc(i) end do if ( nv > maxnv ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST05 - Error!' write ( *, '(a)' ) ' NV > MAXNV.' return end if read ( inunit, * ) ( ivrt(i), i = 1, nv ) nsc = 0 do i = 1, nv if ( ivrt(i) < 0 ) then nsc = nsc + 1 end if end do nsc = nsc / 2 if ( nsc > maxed ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST05 - Error!' write ( *, '(a)' ) ' NSC > MAXED.' 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 data structures. ! write ( *, '(a)' ) ' ' 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,nhole*2+nhola,(holv(i),i=1,nhole*2+nhola) 630 format ((1x,i7,2f15.7)) 640 format (/1x,2i7/(1x,10i7)) 650 format (/(1x,10i7)) 660 format (/1x,i7/(1x,5i7,f15.7)) return end subroutine test06 ! !*********************************************************************** ! !! TEST06 tests DTRIS2; ! implicit none ! integer, parameter :: maxnp = 10000 integer, parameter :: maxst = 100 ! integer a integer b double precision 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) double precision vcl(2,maxnp+3) ! 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 do i = 1, npt ind(i) = i end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST06' write ( *, '(a,i6)' ) ' MSGLVL = ', msglvl write ( *, '(a,g14.6)' ) ' BINEXP = ', binexp if ( npt > maxnp ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST06 - Error!' write ( *, '(a)' ) ' NPT > MAXNP.' 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 call dtris2 ( npt, maxst, vcl, ind, ntri, til, tnbr, stack, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST06 - Error!' write ( *, '(a,i6)' ) ' IERROR = ', ierror return end if nlo = 0 do i = 1, ntri do j = 1, 3 k = tnbr(j,i) if ( k > i ) then jp1 = j + 1 if ( jp1 > 3) then jp1 = 1 end if jp2 = jp1 + 1 if ( jp2 > 3) 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)' ) 'The number of triangles in the triangulation, NTRI = ', ntri write ( *, '(a,i6)' ) 'NLO = ', nlo write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I, TIL(1:3,I), TNBR(1:3,I)' write ( *, '(a)' ) ' ' do i = 1, ntri write ( *, '(7i8)' ) i, (til(j,i),j=1,3), (tnbr(j,i),j=1,3) end do return end subroutine test07 ! !*********************************************************************** ! !! TEST07 tests DTRIW2; ! implicit none ! double precision, parameter :: large = 1000.0D+00 integer, parameter :: maxnp = 10000 integer, parameter :: maxst = 100 ! integer a integer alg integer b double precision 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) double precision 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 ( npt > maxnp ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST07 - Error!' write ( *, '(a)' ) ' NPT > MAXNP.' 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 ( k > i ) then jp1 = j + 1 if ( jp1 > 3) then jp1 = 1 end if jp2 = jp1 + 1 if ( jp2 > 3) 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 ) ! !*********************************************************************** ! !! TEST08 tests EQDIS2; !! TEST08 tests INTPG; !! TEST08 tests MFDEC2; !! TEST08 tests MMASEP; !! TEST08 tests SEPMDF; !! TEST08 tests SEPSHP; !! TEST08 tests SFDWMF; !! TEST08 tests SFUPMF; !! TEST08 tests TRISIZ. ! implicit none ! 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 ! double precision angmin double precision angspc double precision angtol double precision area(maxhv) integer case double precision degrees_to_radians double precision dmin integer edge(4,maxed) character ( len = * ) filename double precision h(maxhv) logical hflag integer ht(0:maxed-1) integer htsiz integer hvl(maxhv) integer i double precision iang(maxpv) integer icur(maxnc) integer ierror integer ivrt(maxnv) integer iwk(maxiw) integer j double precision 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 double precision d_pi integer prime double precision psi(maxhv) integer pvl(4,maxpv) double precision radians_to_degrees integer regnum(maxhv) character ( len = 20 ) rgname double precision tol double precision tolin double precision umdf2 double precision vcl(2,maxvc) double precision 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 = ( kappa >= 0.0D+00 .and. kappa <= 1.0D+00 ) read ( inunit, * ) case,nvc,ncur,msglvl if ( nvc > maxvc ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST08 - Error' write ( *, '(a)' ) ' NVC > MAXVC.' return end if if ( ncur > maxnc ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST08 - Error' write ( *, '(a)' ) ' NCUR > MAXNC.' 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 ( nv > maxnv ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST08 - Error' write ( *, '(a)' ) ' NV > MAXNV.' 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 ( nsc > maxed ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST08 - Error' write ( *, '(a)' ) ' NSC > MAXED.' 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 ( iang(i) > d_pi ( ) + tol ) 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) 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 ! !*********************************************************************** ! !! TEST09 tests LRLINE. ! implicit none ! double precision dv integer lr integer lrline double precision xu double precision xv1 double precision xv2 double precision yu double precision yv1 double precision 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 ! !*********************************************************************** ! !! TEST10 tests LUFAC; !! TEST10 tests LUSOL. ! implicit none ! integer, parameter :: nmax = 20 ! double precision a(nmax,nmax) double precision b(nmax) double precision emax double precision esum integer i integer ipvt(nmax) integer j integer n integer seed logical singlr double precision t double precision tol real 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) = dble ( 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 ) ! !*********************************************************************** ! !! TEST11 tests DSMDF2; !! TEST11 tests MDF2; !! TEST11 tests PRMDF2. ! implicit none ! 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 ! double precision angspc double precision angtol double precision area(maxhv) integer case double precision degrees_to_radians integer edge(4,maxed) double precision edgval(maxpv) character ( len = * ) filename integer ht(0:maxed-1) integer htsiz integer hvl(maxhv) integer i double precision iang(maxpv) integer icur(maxnc) integer ierror integer ifv integer ivrt(maxpv) integer iwk(maxiw) integer j integer map(maxnc) double precision 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 double precision tol double precision tolin double precision val(maxhv) double precision vcl(2,maxvc) double precision vrtval(maxvc) double precision widsq(maxhv) double precision wk(maxwk) double precision x integer xivrt(maxhv+1) double precision 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 ( nvc > maxvc ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST11 - Error' write ( *, '(a)' ) ' NVC > MAXVC.' return end if if ( ncur > maxnc ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST11 - Error' write ( *, '(a)' ) ' NCUR > MAXNC.' 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 ( nv > maxpv ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST11 - Error' write ( *, '(a)' ) ' NV > MAXPC.' return end if read ( inunit, * ) (ivrt(i),i=1,nv) nsc = 0 do i = 1, nv if ( ivrt(i) < 0 ) then nsc = nsc + 1 end if end do nsc = nsc/2 if ( nsc > maxed ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST11 - Error' write ( *, '(a)' ) ' NSC > MAXED.' 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) 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. ! implicit none ! 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 ) ! !*********************************************************************** ! !! TEST13 tests PTPOLG. ! implicit none ! integer, parameter :: inunit = 1 integer, parameter :: maxn = 100 ! double precision a double precision b double precision c integer dim double precision dtol character ( len = * ) filename integer i integer inout integer j integer n double precision nrml(3) integer pgind(0:maxn) double precision pt(3) double precision tol double precision 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 ( n > maxn ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DRPTPG - Error' write ( *, '(a)' ) ' N > MAXN.' 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 ! !*********************************************************************** ! !! TEST14 tests ROTIAR. ! implicit none ! 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 ) ! !*********************************************************************** ! !! TEST15 tests DIAM2; !! TEST15 tests SHRNK2; !! TEST15 tests WIDTH2. ! implicit none ! integer, parameter :: inunit = 1 integer, parameter :: maxn = 100 ! double precision diamsq character ( len = * ) filename integer i integer i1 integer i2 integer iedge(0:maxn) integer ierror integer n integer nshr double precision sdist(0:maxn-1) double precision widsq double precision xc(0:maxn) double precision xs(0:maxn) double precision yc(0:maxn) double precision ys(0:maxn) ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST15' write ( *, '(a)' ) ' Reading data file:' write ( *, '(a)' ) filename open ( unit = inunit, file = filename, form = 'formatted' ) read ( inunit, * ) n if ( n > maxn ) 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 call shrnk2 ( n, xc, yc, sdist, nshr, xs, ys, iedge ) 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 ! !*********************************************************************** ! !! TEST16 tests DHPSRT; !! TEST16 tests IHPSRT; !! TEST16 tests RANDPT. ! implicit none ! integer, parameter :: maxk = 4 integer, parameter :: maxn = 100 ! integer axis double precision da(maxk,maxn) integer i integer ia(maxk,maxn) logical iflag integer j integer k integer map(maxn) integer n integer nptav double precision, dimension ( maxk ) :: scale = (/ & 1.0D+00, 1.0D+00, 1.0D+00, 1.0D+00 /) integer seed double precision, 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. k > maxk ) then return end if if ( n < 1 .or. n > maxn ) 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 ) ! !*********************************************************************** ! !! TEST17 tests BEDGMV; !! TEST17 tests CVDTRI; !! TEST17 tests FNDTRI; !! TEST17 tests INTTRI; !! TEST17 tests LOP; !! TEST17 tests MTREDG; !! TEST17 tests ROTPG; !! TEST17 tests TMERGE; !! TEST17 tests TRINBR; !! TEST17 tests TRIPR2; !! TEST17 tests TRPOLG. ! implicit none ! 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) double precision anga double precision angb double precision angc double precision angle double precision angmax double precision angmin double precision angspc double precision angtol double precision area(maxhv) integer b integer c integer case double precision degrees_to_radians double precision delta double precision dmin integer edge(4,maxed) character ( len = * ) filename double precision h(maxhv) logical hflag integer ht(0:maxed-1) integer htsiz integer hvl(maxhv) integer i double precision iang(maxpv) integer icur(maxnc) integer ierror integer ivrt(maxnv) integer iwk(maxiw) integer j double precision 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 double precision d_pi integer prime double precision psi(maxhv) integer pvl(4,maxpv) double precision radians_to_degrees integer regnum(maxhv) character ( len = 20 ) rgname integer til(3,maxti) integer tnbr(3,maxti) double precision tol double precision tolin integer tstart(maxhv) double precision, external :: umdf2 double precision vcl(2,maxvc) integer vnum(maxpv) integer vstart(maxpv) double precision 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:' 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 ( *, 710 ) rgname, tol, angspc, angtol, kappa, dmin, nmin, ntrid angspc = degrees_to_radians ( angspc ) angtol = degrees_to_radians ( angtol ) hflag = ( kappa >= 0.0D+00 .and. kappa <= 1.0D+00 ) read ( inunit, * ) case, nvc, ncur, msglvl if ( nvc > maxvc ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST17 - Error' write ( *, '(a)' ) ' NVC > MAXVC.' return end if if ( ncur > maxnc ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST17 - Error' write ( *, '(a)' ) ' NCUR > MAXNC.' 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 ( nv > maxnv ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST17 - Error' write ( *, '(a)' ) ' NV > MAXNV.' 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 ( nsc > maxed ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST17 - Error' write ( *, '(a)' ) ' NSC > MAXED.' 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 ( iang(i) > d_pi ( ) + tol ) 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 ) 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. ! 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. ! call trinbr ( nvc, ntri, til, tnbr, htsiz, maxed, ht, edge ) 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 = d_pi ( ) angmax = 0.0D+00 delta = d_pi ( ) / dble(nfreq) 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 = d_pi ( ) - anga - angb angmin = min ( angmin, anga, angb, angc ) angmax = max ( angmax, anga, angb, angc ) if ( abs ( anga - 0.5D+00 * d_pi ( ) ) <= tol ) then anga = 0.5D+00 * d_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 * d_pi ( ) ) <= tol ) then angb = 0.5D+00 * d_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 * d_pi ( ) ) <= tol ) then angc = 0.5D+00 * d_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/dble(nfreq),afreq(i)/dble(3*ntri),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 ! !*********************************************************************** ! !! TEST18 tests VISVRT; !! TEST18 tests VORNBR. ! implicit none ! integer, parameter :: maxn = 200 ! double precision angle double precision angspc double precision 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 double precision phi double precision theta(0:maxn) double precision x(maxn) double precision xc(0:maxn) double precision xeye double precision xvor(0:maxn) double precision y(maxn) double precision yc(0:maxn) double precision yeye double precision 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 ( n > maxn ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST18 - Error' write ( *, '(a)' ) ' N > MAXN.' 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 ( maxnv > maxn ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST18 - Error' write ( *, '(a)' ) ' MAXNV > MAXN.' 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 ! !*********************************************************************** ! !! TEST19 tests ROTIPG; !! TEST19 tests VISPOL. ! implicit none ! integer, parameter :: maxn = 200 ! integer i integer ierror integer ivert integer ivis(0:maxn+1) integer j integer n integer nvis integer nvrt integer vptype double precision x(maxn) double precision xc(0:maxn+1) double precision xeye double precision y(maxn) double precision yc(0:maxn+1) double precision yeye ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST19' write ( *, '(a)' ) ' TEST BEING SKIPPED FOR NOW.' return read ( *, * ) n if ( n < 3 ) then return end if if ( n > maxn ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST19 - Error' write ( *, '(a)' ) ' N > MAXN.' 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 IVERT > 0. ! 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 if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST19 - Error.' write ( *, '(a,i6)' ) ' IERROR = ', ierror return 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 ! !*********************************************************************** ! !! PCDEC generates PICTEX commands from DRDEC or DREQD output. ! implicit none ! 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 double precision 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 double precision vcl(2,maxnv) double precision x1 double precision x2 double precision xmax double precision xmin double precision y1 double precision y2 double precision ymax double precision 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)' ) '\\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) 630 format ('\\plot ',4f15.7,' /') return end subroutine pcds ! !*********************************************************************** ! !! PCDS generates PICTEX commands from DRDS output. ! implicit none ! 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 double precision 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 double precision vcl(2,maxnv) double precision xmax double precision xmin double precision ymax double precision 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 ! !*********************************************************************** ! !! PCTRI generates PICTEX commands from DRTRI output. ! implicit none ! 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) double precision vcl(2,maxvc) double precision x double precision xmax double precision xmin double precision ymax double precision 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 ( vcl(1,a) >= xmin .and. vcl(1,a) <= xmax .and. & vcl(2,a) >= ymin .and. vcl(2,a) <= ymax .and. & vcl(1,b) >= xmin .and. vcl(1,b) <= xmax .and. & vcl(2,b) >= ymin .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 ! !*********************************************************************** ! !! PITRI generates PIC commands from DRTRI output. ! implicit none ! 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) double precision vcl(2,maxvc) double precision x double precision xmax double precision xmin double precision ymax double precision 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 ( vcl(1,a) >= xmin .and. vcl(1,a) <= xmax .and. & vcl(2,a) >= ymin .and. vcl(2,a) <= ymax .and. & vcl(1,b) >= xmin .and. vcl(1,b) <= xmax .and. & vcl(2,b) >= ymin .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 ! !*********************************************************************** ! !! TEST20 tests XEDGE. ! implicit none ! logical intsct integer mode double precision xu double precision xv1 double precision xv2 double precision xw1 double precision xw2 double precision yu double precision yv1 double precision yv2 double precision yw1 double precision 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 ! !*********************************************************************** ! !! TEST21 tests XLINE. ! implicit none ! double precision dv double precision dw logical parall double precision xu double precision xv1 double precision xv2 double precision xw1 double precision xw2 double precision yu double precision yv1 double precision yv2 double precision yw1 double precision 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