c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12
*MAIN
c
c     main program to -
c
c     1. read/compute a Delaunay triangulation of a set of points;
c        if a Delaunay triangulation is computed, its dual is a
c        Voronoi diagram of the set of points; if a Delaunay
c        triangulation is read that is not constrained, i. e.,
c        no line segments were inserted into it, its dual is a
c        Voronoi diagram of the set of points;
c
c     2. insert a set of points and/or a set of line segments into
c        current Delaunay triangulation and obtain a Delaunay
c        triangulation that contains the inserted points and/or
c        line segments; if only points are inserted and the
c        resulting Delaunay triangulation is not constrained,
c        its dual is a Voronoi diagram;
c
c     3. delete a set of points that are vertices in the current
c        Delaunay triangulation and obtain a Delaunay triangulation
c        that does not contain the deleted points; if the resulting
c        Delaunay triangulation is not constrained, its dual is
c        a Voronoi diagram;
c
c     4. compute circumcenters (Voronoi vertices) of triangles in
c        final Delaunay triangulation; if the Delaunay triangulation
c        is not constrained, its dual is a Voronoi diagram and these
c        are its vertices, one per triangle; 
c
c     5. insert edges of a rectangle into Voronoi diagram of set of
c        points (rectangle must contain set of points in its interior)
c        so that parts of Voronoi cells outside rectangle (in
c        particular unbounded cells) are clipped by these edges; the
c        result is the Voronoi diagram in the interior of the rectangle
c        together with the boundary of the rectangle; it is assumed
c        that the corresponding Delaunay triangulation is not
c        constrained so that the Voronoi diagram is well defined and
c        it is also assumed that the circumcenters (Voronoi vertices)
c        of triangles in Delaunay triangulation are computed.
c        
c
c     Author: J. Bernal
c
c     Disclaimer:
c
c     This software was developed at the National Institute of Standards
c     and Technology by employees of the Federal Government in the
c     course of their official duties. Pursuant to title 17 Section 105
c     of the United States Code this software is not subject to
c     copyright protection and is in the public domain. This software is
c     experimental. NIST assumes no responsibility whatsoever for its
c     use by other parties, and makes no guarantees, expressed or
c     implied, about its quality, reliability, or any other
c     characteristic. We would appreciate acknowledgement if the
c     software is used.
c
c**********************************************************************
c
c     Description of some of the variables and parameters used
c     --------------------------------------------------------
c    
c     nv     Number of points or vertices in triangulation including
c            duplicate and deleted points or vertices.
c
c     ivnxt  Number of triangles in triangulation.
c
c     nvmz   Set in program to the value -nvmax-1 (nvmax is described
c            below). See description of arrays is, icon for its use.
c
c     nvmv   Used when program reads information about an existing
c            triangulation from a file. Value is the value that nvmz
c            had during the execution of this program that created
c            the triangulation about which information is being read.
c
c     np     Number of points or vertices in xy input files that will
c            not be processed when computing initial triangulation
c            from scratch and/or inserting an additional set of
c            vertices into the current triangulation (see description
c            of array ip below).
c
c     nl     Number of line segments to be inserted.
c
c     lflag  Integer flag indicating the presence of inserted line
c            segments in triangulation. Equals 1 if any. Zero if none.
c
c     lflug  Integer flag indicating Voronoi cells have been clipped
c            by the edges of a rectangle that contains the set of
c            points for which Voronoi diagram was computed. Equals 1
c            if this is the case. Zero otherwise.
c
c     nd     Number of vertices to be deleted from triangulation.
c
c     nmax   Parameter which is the dimension of arrays x, y, is
c            (described below). Value is the maximum allowed number
c            of points or vertices in triangulation. Value at least
c            that of nv.
c
c     nvmax  Parameter which is the second dimension of array icon
c            (described below) and the dimension of arrays ifl, id, vx,
c            vy (described below). Value is the maximum allowed number
c            of triangles in triangulation. Value at least that of
c            ivnxt. Value at least two times that of nmax for most runs.
c
c     nlmax  Parameter which is the dimension of arrays xrg, yrg,
c            xds, yds (described below). Value is the maximum allowed
c            number of line segments to be inserted into triangulation.
c            Value at least that of nl.
c
c     ndmax  Parameter which is the dimension of arrays xdl, ydl
c            (described below). Value is the maximum allowed number
c            of vertices to be deleted from triangulation. Value at
c            least that of nd.
c
c     nfmax  Parameter which is the dimension of several arrays used
c            internally during the insertion of line segments and
c            deletion of vertices (see below). Value at some multiple
c            of the square root of nmax (described above) seems to work.
c
c     nhmax  Parameter which is the dimension of array ia used
c            internally (see below). Value at least the number of
c            triangles that are affected in triangulation every time
c            a line segment is inserted or a vertex is deleted.
c            Value at some multiple of the square root of nmax
c            (described above) seems to work.
c
c     nkmax  Parameter which is the dimension of array iaze used
c            internally (see below). Value equal to the number of
c            vertices that are deleted plus 1 works.
c
c     njmax  Parameter which is the dimension of array iave used
c            internally (see below). Value at two times that of nkmax
c            seems to work.
c
c     npmax  Parameter which is the dimension of array ip (described
c            below). Value at least equal to the number of points in xy
c            input files that will not be processed when computing the
c            initial triangulation from scratch and/or inserting an
c            additional set of points into the current triangulation.
c
c     epf    Tolerance used by the program to switch from floating
c            point arithmetic to exact arithmetic by testing against
c            this tolerance whether certain quantities are too close
c            to zero; setting it equal to numbers such as 0.1, 0.01
c            has worked well so far.
c
c     x      Array that contains x-coordinates of points or vertices 1
c            through nv. Values obtained by program from input files.
c            Each line of these input files corresponds to one point.
c            Two numbers are found in each line. The first number is
c            the x-coordinate of the point that corresponds to the
c            line. The second number is the y-coordinate of the point
c            (see description of array y below). nmax is the dimension
c            for this array.
c
c     y      Array that contains y-coordinates of points or vertices 1
c            through nv. Values obtained by program from input files
c            (see description of array x above). nmax is the dimension
c            for this array.
c
c     xrg    Array that contains x-coordinates of origin endpoints of
c            line segments 1 through nl. Values obtained by program
c            from input file. Each line of this input file corresponds
c            to one line segment. Four numbers are found in each line.
c            The 1st number is the x-coordinate of the origin endpoint
c            of the line segment that corresponds to the line. The 2nd
c            number is the y-coordinate of the origin endpoint of the
c            line segment, the 3rd number is the x-coordinate of the
c            destination endpoint of the line segment, and the 4th
c            number is the y-coordinate of the destination endpoint of
c            the line segment (see description of arrays yrg, xds, yds
c            below). nlmax is the dimension for this array.
c
c     yrg    Array that contains y-coordinates of origin endpoints of
c            line segments 1 through nl. Values obtained by program
c            from input file (see description of array xrg above).
c            nlmax is the dimension for this array.
c
c     xds    Array that contains x-coordinates of destination endpoints
c            of line segments 1 through nl. Values obtained by program
c            from input file (see description of array xrg above).
c            nlmax is the dimension for this array.
c
c     yds    Array that contains y-coordinates of destination endpoints
c            of line segments 1 through nl. Values obtained by program
c            from input file (see description of array xrg above).
c            nlmax is the dimension for this array.
c
c     xdl    Array that contains x-coordinates of points or vertices 1
c            through nd that are to be deleted from triangulation.
c            Values obtained by program from input file. Each line of
c            this input file corresponds to one point. Two numbers are
c            found in each line. The first number is the x-coordinate
c            of the point that corresponds to the line. The second
c            number is the y-coordinate of the point (see description
c            of array ydl below). ndmax is the dimension for this array.
c
c     ydl    Array that contains y-coordinates of points or vertices 1
c            through nd that are to be deleted from triangulation.
c            Values obtained by program from input file (see description
c            of array xdl above). ndmax is the dimension for this array.
c
c     vx     Array that contains on output x-coordinates of
c            circumcenters (Voronoi vertices) of triangles 1 through
c            ivnxt in final Delaunay triangulation. nvmax is the
c            dimension for this array.
c
c     vy     Array that contains on output y-coordinates of
c            circumcenters (Voronoi vertices) of triangles 1 through
c            ivnxt in final Delaunay triangulation. nvmax is the
c            dimension for this array.
c
c     is     Array with dimension nmax. From i = 1 to nv, is(i) is the
c            triangle number or index for a triangle in triangulation
c            with vertex i as one of its vertices. is(i) equals nvmz
c            if vertex i was not processed when it was its turn to be
c            inserted into current triangulation. Equals 0 if it was
c            deleted from the triangulation. Equals a negative integer
c            larger than nvmz if vertex i is a duplication of vertex
c            -is(i) which was a vertex in the triangulation when it was
c            the turn of vertex i to be inserted in the triangulation.
c            If vertex i is on the boundary of the convex hull of the
c            vertices in the triangulation then is(i) is the triangle
c            number or index for the triangle in the triangulation that
c            is the first triangle found in a counterclockwise
c            direction in the "fan" of triangles with vertex i as a
c            vertex.
c
c     icon   Two-dimensional array with 6 as its first dimension
c            and nvmax as its second. For j = 1, ..., ivnxt, icon(i,j),
c            i = 1, ..., 6, provides information about the j-th
c            triangle in the triangulation. For j = 1, ..., ivnxt,
c            icon(i,j), i = 1, 2, 3, are the triangle numbers or
c            indices for the three triangles adjacent to triangle j in
c            the order in which they appear in a counterclockwise
c            direction around triangle j. Given i, i=1, 2, or 3, if
c            icon(i,j) is negative and not equal to nvmz, then
c            -icon(i,j) is the actual index for the ith adjacent
c            triangle to triangle j and the negativity of icon(i,j)
c            simply implies that the edge shared by triangle j and
c            triangle -icon(i,j) is part of an inserted line segment.
c            Given i, i=1, 2, or 3, if icon(i,j) equals either zero or
c            nvmz then the ith adjcent triangle to triangle j does not
c            exist and the edge of triangle j that triangle j would
c            share with that triangle if it existed is part of the
c            boundary of the convex hull of the vertices in the
c            triangulation. If in addition icon(i,j) equals nvmz then
c            the aforementioned edge is part of an inserted line
c            segment. For j = 1, ..., ivnxt, icon(i,j), i = 4, 5, 6,
c            are the vertex numbers or indices for the three vertices
c            of triangle j in the order in which they appear in a
c            counterclockwise direction around triangle j. In addition,
c            assuming without any loss of generality that for i=1, 2,
c            3, icon(i,j) is positive, for i= 4, 5, 6, icon(i,j) is the
c            vertex number for the vertex of triangle j that is not a
c            vertex for triangle icon(i-3,j).
c
c     ifl    Array with dimension nvmax. For i = 1,..., ivnxt, ifl(i)
c            will equal 1 if the size of Voronoi vertex i is not
c            considered to be too large. If the size of vertex i is
c            considered to be too large ifl(i) will equal 0.
c
c     id     Array with dimension nvmax used internally.
c
c     ix
c     iy
c     ix2
c     iy2    Arrays of dimension nmax used internally.
c
c     ia     Array of dimension nhmax used internally.
c
c     iave   Array of dimension njmax used internally.
c
c     iaze   Array of dimension nkmax used internally.
c
c     ifun
c     irun
c     isun
c     izun
c     igun
c     iwun   Arrays of dimension nfmax used internally.
c
c     ip     Array that contains indices of points or vertices in xy
c            input files that will not be processed when computing
c            initial triangulation from scratch and/or inserting an
c            additional set of points into the current triangulation.
c            For i = 1, ..., np, ip(i) is the i-th point not to be
c            processed. Values obtained by program from an input file.
c            npmax is the dimension of this array.
c
c***********************************************************************
c
      program main
c
      integer nmax, nvmax, nlmax, ndmax
      integer nfmax, nhmax, nkmax, njmax, npmax
      parameter (nmax = 160000, nvmax = 2.1*nmax)
      parameter (nlmax = 10000, ndmax = 30000)
      parameter (nfmax = 5 * 400, nhmax = 5 * 400)
      parameter (nkmax = ndmax+1, njmax = 2 * nkmax, npmax = 10000)
c
      double precision x(nmax), y(nmax), vx(nvmax), vy(nvmax)
      double precision xrg(nlmax), yrg(nlmax), xds(nlmax), yds(nlmax)
      double precision xdl(ndmax), ydl(ndmax)
      integer ix(nmax), iy(nmax), ix2(nmax), iy2(nmax)
      integer is(nmax), icon(6,nvmax), ifl(nvmax), id(nvmax), ia(nhmax)
      integer iave(njmax), iaze(nkmax), ifun(nfmax), irun(nfmax)
      integer isun(nfmax), izun(nfmax), igun(nfmax), iwun(nfmax)
      integer ip(npmax)
c
      integer n, nv, nr, ivnxt, nvmz, nvmv, np, nl, lflag, lflug, nd
      integer icsfig, mhalf, mfull, isclp(2)
      double precision r215, deps, dscle, dfull, dfill, epf
c
      character*30 vninf, vsinf, voinf, lsinf, dlinf, rtinf,
     *             tninf, toinf, ccinf, npinf, trinf
      character*1 answ
      logical vnflg, vsflg, voflg, lsflg, dlflg, rtflg,
     *        toflg, ccflg, npflg, trflg
      logical prompt
c
      integer i, j, nvmn, nvmp, iase, ileft, iabe, idupl, inpro
      integer nvp, nls, nlp, ndl, nde, ndv, nds
      double precision xcor, ycor, xmax, xmin, ymax, ymin
      double precision xdes, ydes, xorg, yorg
c
      integer ng, ifir, isec, ithi, larg
      double precision derro, dist1, dist2, dist3
      double precision dista, disti, diffd, dnom
c
      double precision xne, yne, xnw, ynw, xsw, ysw, xse, yse
c
      write(*,*)' '
      write(*,*)'This program is for reading/computing a Delaunay ',
     *          'triangulation of'
      write(*,*)'a 2-dimensional set of points.'
      write(*,*)' '
      write(*,*)'This program is also for inserting points and/or ',
     *          'line segments'
      write(*,*)'into current Delaunay triangulation, and obtaining ',
     *          'a new Delaunay'
      write(*,*)'triangulation that contains the new points ',
     *          'and/or line segments.'
      write(*,*)' '
      write(*,*)'This program is also for deleting vertices from ',
     *          'current Delaunay'
      write(*,*)'triangulation, and obtaining a new Delaunay ',
     *          'triangulation that'
      write(*,*)'does not contain them.'
      write(*,*)' '
      write(*,*)'This program is also for computing circumcenters ',
     *          '(Voronoi vertices)'
      write(*,*)'of triangles in final Delaunay triangulation.'
      write(*,*)' '
      write(*,*)'If Voronoi vertices are computed and Delaunay ',
     *          'triangulation is'
      write(*,*)'unconstrained (no inserted line segments), its ',
     *          'dual (Voronoi'
      write(*,*)'diagram) is well defined, and then this program is ',
     *          'also for'
      write(*,*)'inserting edges of a rectangle into Voronoi diagram ',
     *          '(a rectangle'
      write(*,*)'that contains in its interior set of points for ',
     *          'which Voronoi'
      write(*,*)'diagram is defined) so that parts of Voronoi cells ',
     *          'outside'
      write(*,*)'rectangle (mainly unbounded cells) are ',
     *          'clipped by these edges.' 
c
c----------------------------------------------------------------------
c
c     inquire about existence of prompts file
c
      inquire(file = 'prompts_deltri', exist = prompt)
      write(*,*)' '
      if(prompt)then
         write(*,*)'User intervention not required: responses to ',
     *             'would-be-program-'
         write(*,*)'generated prompts will be read from file ',
     *             'prompts_deltri, thus'
         write(*,*)'prompts will be suppressed.'
         open (9, file = 'prompts_deltri')
      else
         write(*,*)'User intervention required: file named ',
     *             'prompts_deltri does not'
         write(*,*)'exist, thus program will generate prompts for ',
     *             'which user must'
         write(*,*)'provide responses. In order to avoid user ',
     *             'intervention a file'
         write(*,*)'named prompts_deltri must exist that contains ',
     *             'responses to the'
         write(*,*)'would-be prompts; prompts are then suppressed; ',
     *             'for this purpose'
         write(*,*)'the program will generate a file named ',
     *             'tentative_prompts_deltri'
         write(*,*)'that will contain the responses provided ',
     *             'by the user during the'
         write(*,*)'currrent run and which can then be renamed ',
     *             'prompts_deltri by the'
         write(*,*)'user for future runs.'
         open (10, file = 'tentative_prompts_deltri')
      endif
c
c----------------------------------------------------------------------
c
c     initialize
c
      vnflg=.false.
      vsflg=.false.
      voflg=.false.
      lsflg=.false.
      dlflg=.false.
      rtflg=.false.
      toflg=.false.
      ccflg=.false.
      npflg=.false.
      trflg=.false.

      do 20 i = 1, 30
          vninf(i:i) = ' '
          vsinf(i:i) = ' '
          voinf(i:i) = ' '
          lsinf(i:i) = ' '
          dlinf(i:i) = ' '
          rtinf(i:i) = ' '
          tninf(i:i) = ' '
          toinf(i:i) = ' '
          ccinf(i:i) = ' '
          npinf(i:i) = ' '
          trinf(i:i) = ' '
   20 continue

      do 25 i = 1, nmax
         is(i) = 0
   25 continue

   30 format(a1)
   50 format(a30)

      answ = ' '
      if(prompt)then
         write(*,*)' '
         read(9,30) answ
      else
         write(*,*)' '
         write(*,*)'Do you want Delaunay triangulation to be computed ',
     *             'from scratch?(y/n)'
         read(5,30) answ
         write(10,30) answ
      endif
      if(answ.eq.'y'.or.answ.eq.'Y') then
         vnflg=.true.
         write(*,*)'Delaunay triangulation will be computed from ',
     *             'scratch.'
         if(prompt)then
            read(9,50)vninf
         else
            write(*,*)' '
            write(*,*)'Specify name of file with points ',
     *                '(x-, y-coordinates) for which'
            write(*,*)'triangulation is to be computed:'
            read(5,50)vninf
            write(10,50)vninf
         endif
      else
         write(*,*)'A Delaunay triangulation will be read from a file.'
         if(prompt)then
            read(9,50)tninf
            read(9,50)vninf
         else
            write(*,*)' '
            write(*,*)'Specify name of file with triangulation ',
     *                'information:'
            read(5,50)tninf
            write(10,50)tninf
            write(*,*)'Specify name of file with points ',
     *                '(x-, y-coordinates) that are'
            write(*,*)'the vertices of this triangulation:'
            read(5,50)vninf
            write(10,50)vninf
         endif
      endif

      if(.not. vnflg) go to 55
      if(prompt)then
         read(9,*) icsfig
      else
         write(*,*)' '
         write(*,*)'Enter icsfig, i. e., # of significant figures ',
     *             'of decimal part'
         write(*,*)'of point coordinates, -1 < icsfig < 10 (total ',
     *             'number of significant'
         write(*,*)'figures should be at most 14 with at most 9 to ',
     *             'either the left or'
         write(*,*)'the right of the decimal point):'
         read(5,*) icsfig
         write(10,*) icsfig
      endif
   55 continue

      if(prompt) then
         read(9,30) answ
      else
         write(*,*)' '
         write(*,*)'Do you want an additional set of points to be ',
     *             'inserted as vertices'
         write(*,*)'into current triangulation?(y/n)'
         read(5,30) answ
         write(10,30) answ
      endif
      if(answ.eq.'y'.or.answ.eq.'Y') then
         vsflg=.true.
         if(prompt) then
            read(9,50) vsinf
         else
            write(*,*)'Specify name of file with additional points ',
     *                '(x-, y-coordinates):'
            read(5,50) vsinf
            write(10,50) vsinf
         endif
      endif

      if(.not. vnflg .and. .not. vsflg) go to 60
      if(prompt) then
         read(9,30) answ
      else
         write(*,*)' '
         write(*,*)'Do you want certain vertices in xy input files ',
     *             'not to be processed'
         write(*,*)'during the insertion of vertices into current ',
     *             'triangulation?(y/n)'
         read(5,30) answ
         write(10,30) answ
      endif
      if(answ.eq.'y'.or.answ.eq.'Y') then
         npflg=.true.
         if(prompt) then
            read(9,50) npinf
         else
            write(*,*)'Specify name of file with indices of points ',
     *                'not to be processed:'
            read(5,50) npinf
            write(10,50) npinf
         endif
      endif
   60 continue

      if(prompt) then
         read(9,30) answ
      else
         write(*,*)' '
         write(*,*)'Do you want a set of line segments to be ',
     *             'inserted as edges into'
         write(*,*)'current triangulation?(y/n)'
         read(5,30) answ
         write(10,30) answ
      endif
      if(answ.eq.'y'.or.answ.eq.'Y') then
         lsflg=.true.
         if(prompt) then
            read(9,50) lsinf
         else
            write(*,*)'Specify name of file with endpoints ',
     *                '(x-, y-coordinates)'
            write(*,*)'of line segments:'
            read(5,50) lsinf
            write(10,50) lsinf
         endif
      endif

      if(prompt) then
         read(9,30) answ
      else
         write(*,*)' '
         write(*,*)'Do you want a set of points to be deleted as ',
     *             'vertices from the'
         write(*,*)'current triangulation?(y/n)'
         read(5,30) answ
         write(10,30) answ
      endif
      if(answ.eq.'y'.or.answ.eq.'Y') then
         dlflg=.true.
         if(prompt) then
            read(9,50) dlinf
         else
            write(*,*)'Specify name of file with points ',
     *                '(x- ,y-coordinates) to be '
            write(*,*)'deleted as vertices from current triangulation:'
            read(5,50) dlinf
            write(10,50) dlinf
         endif
      endif

      if(prompt) then
         read(9,30) answ
      else
         write(*,*)' '
         write(*,*)'Do you want circumcenters (Voronoi vertices) of ',
     *             'triangles in'
         write(*,*)'final triangulation to be computed and saved ',
     *             'in a file?(y/n)'
         read(5,30) answ
         write(10,30) answ
      endif
      if(answ.eq.'y'.or.answ.eq.'Y') then
         ccflg=.true.
         if(prompt) then
            read(9,50) ccinf
         else
            write(*,*)'Specify name of ouput circumcenters file:'
            read(5,50) ccinf
            write(10,50) ccinf
         endif
      endif

      if(prompt) then
         read(9,30) answ
      else
         write(*,*)' '
         write(*,*)'Do you want edges of a rectangle to be inserted ',
     *             'into Voronoi'
         write(*,*)'diagram so that unbounded Voronoi cells are ',
     *             'clipped by these'
         write(*,*)'edges?'
         read(5,30) answ
         write(10,30) answ
      endif
      if(answ.eq.'y'.or.answ.eq.'Y') then
         rtflg=.true.
         if(prompt) then
            read(9,50) rtinf
         else
            write(*,*)'Specify name of file with 4 corners ',
     *                '(x- ,y-coordinates) of'
            write(*,*)'rectangle to be inserted into Voronoi diagram ',
     *                '(rectangle must'
            write(*,*)'contain in its interior set of points for ',
     *                'which diagram is'
            write(*,*)'being computed, and corners must be ',
     *                'in order NE, NW, SW, SE):'
            read(5,50) rtinf
            write(10,50) rtinf
         endif
      endif

      if(prompt) then
         read(9,30) answ
      else
         write(*,*)' '
         write(*,*)'Do you want information about final Delaunay ',
     *             'triangulation to be'
         write(*,*)'saved in a file?(y/n)'
         read(5,30) answ
         write(10,30) answ
      endif
      if(answ.eq.'y'.or.answ.eq.'Y') then
         toflg=.true.
         if(prompt) then
            read(9,50) toinf
         else
            write(*,*)'Specify name of output triangulation file:'
            read(5,50) toinf
            write(10,50) toinf
         endif
      endif

      if(prompt) then
         read(9,30) answ
      else
         write(*,*)' '
         write(*,*)'Do you want triangles in final triangulation to ',
     *             'be saved in a'
         write(*,*)'file by themselves?(y/n)'
         write(*,*)'(Each line in this file would simply contain 3 ',
     *             'integers: the indices'
         write(*,*)'of the vertices of a triangle in a ',
     *             'counterclowise direction around'
         write(*,*)'the triangle, one triangle per line, each ',
     *             'triangle in final'
         write(*,*)'triangulation in one and only one line)'
         read(5,30) answ
         write(10,30) answ
      endif
      if(answ.eq.'y'.or.answ.eq.'Y') then
         trflg=.true.
         write(*,*)'Triangles will be saved by themselves in a file.'
         if(prompt) then
            read(9,50) trinf
         else
            write(*,*)'Specify name of triangles file:'
            read(5,50) trinf
            write(10,50) trinf
         endif
      endif

      if(prompt) then
         read(9,30) answ
      else
         write(*,*)' '
         write(*,*)'Do you want set of points (x-, y-coordinates) ',
     *             'that are possibly'
         write(*,*)'vertices in final triangulation to be saved in ',
     *             'a file?(y/n)'
         write(*,*)'(Points that are vertices in final ',
     *             'triangulation that may or may'
         write(*,*)'not come from different sources will be put ',
     *             'together in this'
         write(*,*)'file in the order they were read for future use)'
         read(5,30) answ
         write(10,30) answ
      endif
      if(answ.eq.'y'.or.answ.eq.'Y') then
         voflg=.true.
         if(prompt) then
            read(9,50) voinf
         else
            write(*,*)'Specify name of output point file:'
            read(5,50) voinf
            write(10,50) voinf
         endif
      endif

c
      nvmz = -nvmax - 1
      nvmn = nvmz - 1
      nvmp = -nvmn
c     WRITE(*,*)'NVMAX=',NVMAX,' NVMZ=',NVMZ
c     WRITE(*,*)'NVMN=',NVMN,' NVMP=',NVMP
c
c     set tolerance, inserted line segments flag, clipping flag
c
      epf = 0.01d0
      lflag = 0
      lflug = 0
c
c*****************************************************************
c
c     read points file
c
      write(*,*)' '
      write(*,*)'Reading points file (x-,y-coordinates) ... '

      open(unit=11,file=vninf,status='old')
      nv = 0
  100 continue
      read(11,*,end=110) xcor, ycor
      nv = nv + 1
      if(nv .gt. nmax) stop 10
      x(nv) = xcor
      y(nv) = ycor
      if(nv .ne. 1) then
          if(xcor .gt. xmax) xmax = xcor
          if(xcor .lt. xmin) xmin = xcor
          if(ycor .gt. ymax) ymax = ycor
          if(ycor .lt. ymin) ymin = ycor
      else
          xmax = xcor
          xmin = xcor
          ymax = ycor
          ymin = ycor
      endif
      go to 100
  110 continue
      close(unit=11)
      if(nv.lt.3) stop 20
      write(*,*)' '
      write(*,*)'Number of points read =',nv
      n=nv
      iase=1
      iaze(1)=nv+1

      write(*,*)' '
      write(*,*)' xmin = ',xmin,' xmax = ',xmax
      write(*,*)' ymin = ',ymin,' ymax = ',ymax
c

      np = 0
      if(npflg)then
         write(*,*)' '
         write(*,*)'Reading indices of points not to be processed ... '
         open(unit=19,file=npinf,status='old')
         read(19,*)np
         if(np.gt.npmax) stop 30
         if(np.gt.0) then
            read(19,*)(ip(i),i=1,np)
            if(.not.vnflg) go to 130
            do 120 i=1,np
               if(ip(i).lt.1.or.ip(i).gt.nv) go to 120
               is(ip(i))=nvmz
  120       continue
  130       continue
         endif
         close(unit=19)
      endif

      if(.not.vnflg) go to 140
c
c*****************************************************************
c
c     transform coordinates into integer decomposition
c
      write(*,*)' '
      write(*,*)'Transforming input point coordinates into ',
     *          'integer decomposition ...'

      call transf(x, y, ix, iy, ix2, iy2, xmin, ymin, xmax, ymax,
     *            nv, icsfig, mhalf, mfull, isclp, r215, deps,
     *            dscle, dfull, dfill, epf)

      go to 200
c
c*****************************************************************
c
c     read triangulation
c
  140 continue
      write(*,*)' '
      write(*,*)'Reading triangulation information ... '

      open(unit=12,file=tninf,status='old')
      read(12,*) nv, ivnxt, nvmv, icsfig, lflag, lflug
      if(lflug.ne.0) then
         write(*,*)'Triangulation that was read corresponds to ',
     *             'clipped Voronoi'
         write(*,*)'diagram and this situation is not allowed.'
         write(*,*)'Program terminated.'
         stop
      endif
      if(ivnxt.gt.nvmax) stop 40
      if(nv.ne.n.or.ivnxt.gt.-nvmv-1) stop 50
      read(12,170)((icon(i,j),i=1,6),j=1,ivnxt)
      read(12,180)(is(i),i=1,nv)
      close(unit=12)
      if(nvmv.eq.nvmz)go to 158
      do 155 j=1,ivnxt
         do 153 i=1,3
            if(icon(i,j).eq.nvmv)icon(i,j)=nvmz
  153    continue
  155 continue
      do 157 i=1, nv
         if(is(i).eq.nvmv) is(i)=nvmz
  157 continue
  158 continue
      do 160 i=1, nv
         if(is(i).le.0) go to 160
         ileft=i
         go to 165
  160 continue
      stop 60
  165 continue
      iabe=1
      iave(1)=ivnxt+1
  170 format(6i10)
  180 format(5i10)
c
c*****************************************************************
c
c     transform coordinates into integer decomposition
c
      write(*,*)' '
      write(*,*)'Transforming input point coordinates into ',
     *          'integer decomposition ...'

      call transf(x, y, ix, iy, ix2, iy2, xmin, ymin, xmax, ymax,
     *            nv, icsfig, mhalf, mfull, isclp, r215, deps,
     *            dscle, dfull, dfill, epf)

      go to 480
c
c*****************************************************************
c
c     compute triangulation
c
  200 continue
c
      write(*,*)' '
      write(*,*)'Computing triangulation ... '
c
      call tricmp(x, y, ix, iy, ix2, iy2, is, icon, iave, iaze,
     *            epf, mhalf, mfull, isclp, nvmz, nv, iabe, iase,
     *            nmax, nvmax, ivnxt, idupl, inpro, r215, deps, ileft)
c
      write(*,*)' '
      write(*,*)'Number of vertices (points) in triangulation: ',nv
      write(*,*)'Number of triangles: ',ivnxt
c
      write(*,*)' '
      write(*,*)'Number of points read that were duplications =',idupl
      write(*,*)'Number of points read that were not processed =',inpro

  480 continue
      if(.not.vsflg) go to 600
c
c**********************************************************************
c
c     insert additional set of points
c
      write(*,*)' '
      write(*,*)'Reading and inserting additional set of points ',
     *          '(x-,y-coordinates) ...'

      open(unit=13,file=vsinf,status='old')
  500 continue
      read(13,*,end=550) xcor, ycor
      nv = nv + 1
      if(nv .gt. nmax) stop 70
      x(nv) = xcor
      y(nv) = ycor
      if(xcor .gt. xmax) xmax = xcor
      if(xcor .lt. xmin) xmin = xcor
      if(ycor .gt. ymax) ymax = ycor
      if(ycor .lt. ymin) ymin = ycor
      go to 500
  550 continue
      close(unit=13)
      nvp=nv-n
      if(nvp.le.0)then
         write(*,*)' '
         write(*,*)'main: (warning) empty additional points file'
         go to 600
      endif
      write(*,*)' '
      write(*,*)'Number of additional points read =',nvp
      iase=1
      iaze(1)=nv+1

      write(*,*)' '
      write(*,*)' xmin = ',xmin,' xmax = ',xmax
      write(*,*)' ymin = ',ymin,' ymax = ',ymax

      call triadd(x, y, ix, iy, ix2, iy2, is, icon, ip, iave, iaze,
     *            epf, icsfig, mhalf, mfull, isclp, nvmz, n, nv, np,
     *            iabe, iase, nmax, nvmax, ivnxt, idupl, inpro, ileft,
     *            dscle, dfull, dfill)
c
      write(*,*)' '
      write(*,*)'Number of vertices (points) in triangulation: ',nv
      write(*,*)'Number of triangles: ',ivnxt
c
      write(*,*)' '
      write(*,*)'Number of additional points read that were ',
     *          'duplications =',idupl
      write(*,*)'Number of additional points read that were ',
     *          'not processed =',inpro
      n=nv

  600 continue

      if(.not.lsflg) go to 700
c
c**********************************************************************
c
c     insert line segments
c
      write(*,*)' '
      write(*,*)'Reading endpoints (x-,y-coordinates) of line ',
     *          'segments ... '

      open(unit=14,file=lsinf,status='old')
      nl = 0
  620 continue
      read(14,*,end=680) xorg, yorg, xdes, ydes
      nl = nl + 1
      if(nl .gt. nlmax) stop 80
      xrg(nl) = xorg
      yrg(nl) = yorg
      xds(nl) = xdes
      yds(nl) = ydes
      go to 620
  680 continue
      close(unit=14)
      if(nl.eq.0) go to 700
c
c     set inserted line segments flag
c
      lflag = 1
      write(*,*)' '
      write(*,*)'Inserting line segments ... '

      call linins(x, y, xrg, yrg, xds, yds, ix, iy, ix2, iy2, is, icon,
     *            iave, iaze, ia, ifun, irun, isun, izun, igun, iwun,
     *            epf, icsfig, mhalf, mfull, isclp, nvmz, nv, nl,
     *            iabe, iase, nmax, nvmax, nfmax, nhmax, ivnxt, nvmp,
     *            nvmn, nls, nlp, ileft, dscle, dfull, dfill)
      if(nls.eq.0.and.nlp.eq.0)then
         write(*,*)' '
         write(*,*)'main: (warning) empty line segment file ...'
         go to 700
      endif
      write(*,*)' '
      write(*,*)'Number of inserted line segments =',nls
      if(nlp.ne.0) then
         write(*,*)' '
         write(*,*)'Number of degenerate line segments (endpoints ',
     *             'identical) =',nlp
      endif
      nvp=nv-n
      write(*,*)' '
      write(*,*)'Number of inserted points due to line segments =',nvp
      n=nv

  700 continue
c     WRITE(*,*)'700 IABE=',IABE,' IASE=',IASE
      if(.not.dlflg) go to 800
c
c**********************************************************************
c
c     delete set of points from triangulation
c
      write(*,*)' '
      write(*,*)'Reading points or vertices to be deleted ... '

      nd = 0
      open(unit=15,file=dlinf,status='old')
  720 continue
      read(15,*,end=780) xdes, ydes
      nd = nd + 1
      if(nd .gt. ndmax) stop 90
      xdl(nd) = xdes
      ydl(nd) = ydes
      go to 720
  780 continue
      close(unit=15)
      if(nd.eq.0) go to 800

      write(*,*)' '
      write(*,*)'Deleting points or vertices from triangulation ... '

      call pntdel(x, y, xdl, ydl, ix, iy, ix2, iy2, is, icon, iave,
     *            iaze, ia, ifun, irun, isun, izun, igun, iwun, epf,
     *            icsfig, mhalf, mfull, isclp, nvmz, nv, nd, iabe, iase,
     *            nmax, nvmax, nfmax, nhmax, njmax, nkmax, ivnxt, nvmp,
     *            nvmn, ndl, nde, ndv, nds, ileft, dscle, dfull, dfill)

      if(ndl.eq.0.and.nde.eq.0.and.ndv.eq.0.and.nds.eq.0)then
        write(*,*)' '
        write(*,*)'main: (warning) empty file with points to be deleted'
        go to 800
      endif
      write(*,*)' '
      write(*,*)'Number of deleted points =',ndl
      if(nde.ne.0) then
         write(*,*)' '
         write(*,*)'Number of points that were not deleted ',
     *             'because they were'
         write(*,*)'in inserted line segments =',nde
      endif
      if(ndv.ne.0) then
         write(*,*)' '
         write(*,*)'Number of points that were not deleted ',
     *             'because they were'
         write(*,*)'not vertices in triangulation =',ndv
      endif
      if(nds.ne.0) then
         write(*,*)' '
         write(*,*)'Number of points that were not deleted ',
     *             'because triangulation'
         write(*,*)'consists of one single triangle =',nds
      endif
c
c     compress array icon
c
      call comprs(is, icon, id, n, ivnxt, nmax, nvmax, nvmz)

  800 continue
c     WRITE(*,*)'800 IABE=',IABE,' IASE=',IASE
      if(.not.ccflg) go to 1040
c
c********************************************************************
c
c     compute circumcenters (x-,y-coordinates) of triangles,
c     i. e., Voronoi vertices, in final triangulation
c
      write(*,*)' '
      write(*,*)'Computing circumcenters (x-,y-coordinates) ',
     *          'of triangles,'
      write(*,*)'i. e., Voronoi vertices, in final triangulation ... '

      call vorcmp(ix, iy, ix2, iy2, vx, vy, icon, ifl, ivnxt, dscle,
     *            mhalf, mfull, isclp)
c
c     test Voronoi vertices
c
      write(*,*)' '
      write(*,*)'Testing circumcenters/Voronoi vertices ...'

      do 1010 i = 1, ivnxt
         if(ifl(i).eq.0) go to 1010
         go to 1015
 1010 continue
      write(*,*)' '
      write(*,*)'All Voronoi vertices marked as too large.'
      go to 1030
 1015 continue
      ng = 0
      derro = 0.0d0
      larg = 0
      do 1020 i = 1, ivnxt
         if(ifl(i).eq.0) then
            larg = larg+1
            go to 1020
         endif
         ifir = icon(4,i)
         isec = icon(5,i)
         ithi = icon(6,i)
         if(ifir.le.ng .or. isec.le.ng .or. ithi.le.ng) stop 100
c
         xcor = x(ifir)-vx(i)
         ycor = y(ifir)-vy(i)
         dnom=dmax1(dabs(xcor),dabs(ycor))
         if(dnom.gt.1.0d0) then
            xcor = xcor/dnom
            ycor = ycor/dnom
            dist1 = dnom*dsqrt(dabs(xcor**2 + ycor**2))
         else
            dist1 = dsqrt(dabs(xcor**2 + ycor**2))
         endif
c
         xcor = x(isec)-vx(i)
         ycor = y(isec)-vy(i)
         dnom=dmax1(dabs(xcor),dabs(ycor))
         if(dnom.gt.1.0d0) then
            xcor = xcor/dnom
            ycor = ycor/dnom
            dist2 = dnom*dsqrt(dabs(xcor**2 + ycor**2))
         else
            dist2 = dsqrt(dabs(xcor**2 + ycor**2))
         endif
c
         xcor = x(ithi)-vx(i)
         ycor = y(ithi)-vy(i)
         dnom=dmax1(dabs(xcor),dabs(ycor))
         if(dnom.gt.1.0d0) then
            xcor = xcor/dnom
            ycor = ycor/dnom
            dist3 = dnom*dsqrt(dabs(xcor**2 + ycor**2))
         else
            dist3 = dsqrt(dabs(xcor**2 + ycor**2))
         endif
c
         dista = dmax1(dist1,dist2,dist3)
         disti = dmin1(dist1,dist2,dist3)
         diffd = dista-disti
         if(diffd.gt.derro) derro = diffd
 1020 continue
      write(*,*)' '
      write(*,*)'Number of Voronoi vertices marked as too ',
     *          'large = ',larg
      write(*,*)' '
      write(*,*)'Voronoi distance error = ',derro
      write(*,*)'(This number should be close to zero).'
 1030 continue

 1040 continue
      if(.not.rtflg) go to 1500
c
c*********************************************************************
c
c     insert edges of a rectangle into Voronoi diagram so that
c     parts of Voronoi cells outside rectangle are clipped by
c     these edges
c
      if(.not.ccflg) then
         write(*,*)' '
         write(*,*)'Voronoi vertices were not computed'
         write(*,*)'Edges of rectagle will not be inserted.'
         rtflg = .false.
         go to 1500
      endif
      write(*,*)' '
      write(*,*)'Inserting edges of a rectangle into Voronoi diagram ',
     *          'so that'
      write(*,*)'Voronoi cell parts outside rectangle are clipped ',
     *          'by edges ...'
      open(unit=21,file=rtinf,status='old')
      read(21,*) xne, yne
      read(21,*) xnw, ynw
      read(21,*) xsw, ysw
      read(21,*) xse, yse
      if(xne.ne.xse .or. xnw.ne.xsw .or. yne.ne.ynw .or.
     *   ysw.ne.yse) stop 110
c     write(*,*)' xmin = ',xmin,' xmax = ',xmax
c     write(*,*)' ymin = ',ymin,' ymax = ',ymax
c     write(*,*)' xne = ',xne,' xnw = ',xnw
c     write(*,*)' yne = ',yne,' yse = ',yse
      if(xne.le.xnw .or. yne.le.yse) stop 120
      if(xne.le.xmax .or. xnw.ge.xmin .or. yne.le.ymax .or.
     *   yse.ge.ymin) stop 130
c
c     test that Delaunay triangulation is not constrained
c
      if(lflag.eq.1) then
         write(*,*)' '
         write(*,*)'Delaunay triangulation appears constrained'
         write(*,*)'Edges of rectangle will not be inserted.'
         rtflg = .false.
         go to 1500
      endif
      do 1050 j = 1, ivnxt
         if(icon(1,j).ge.0 .and. icon(2,j).ge.0 .and.
     *      icon(3,j).ge.0) go to 1050
         stop 135
 1050 continue
c
c     set clipping flag
c
      lflug = 1
c
c     identify Voronoi cells parts of which intersect the exterior
c     of rectangle, parts that will be clipped during insertion
c     of edges of rectangle into Voronoi diagram
c
      write(*,*)' '
      write(*,*)'Identifying Voronoi cells to be clipped ...'
      call vclidn(vx, vy, x, y, icon, is, nr, nv, ivnxt, nmax,
     *            xne, yne, xnw, ynw, xsw, ysw, xse, yse)
c
c     insert into Delaunay triangulation additional set of points
c     that are reflections (computed above in routine vdlidn) of
c     points whose Voronoi cells are to be clipped in the proper
c     order along edges of rectangle in such a way that the final
c     Voronoi diagram will contain the edges of the rectangle as
c     the union of edges of Voronoi cells
c
      nvp = nr-nv
      if(nvp.le.0) stop 140
      write(*,*)' '
      write(*,*)'Number of computed reflection points =',nvp
      iase=1
      iaze(1)=nr+1
      write(*,*)' '
      write(*,*)'Inserting reflection points ...'

      call triadd(x, y, ix, iy, ix2, iy2, is, icon, ip, iave, iaze,
     *            epf, icsfig, mhalf, mfull, isclp, nvmz, nv, nr, np,
     *            iabe, iase, nmax, nvmax, ivnxt, idupl, inpro, ileft,
     *            dscle, dfull, dfill)
      n = nv
      nv = nr
c
      write(*,*)' '
      write(*,*)'Number of vertices (points) in triangulation: ',nv
      write(*,*)'Number of triangles: ',ivnxt
c
      write(*,*)' '
      write(*,*)'Number of reflection points that were ',
     *          'duplications =',idupl
      write(*,*)'Number of reflection points that were ',
     *          'not processed =',inpro
c
c     compute circumcenters (x-,y-coordinates) of triangles,
c     i. e., Voronoi vertices, in final triangulation that
c     includes reflection points
c
      write(*,*)' '
      write(*,*)'Computing circumcenters (x-,y-coordinates) ',
     *          'of triangles,'
      write(*,*)'i. e., Voronoi vertices, in final triangulation that'
      write(*,*)'includes reflection points ... '

      call vorcmp(ix, iy, ix2, iy2, vx, vy, icon, ifl, ivnxt, dscle,
     *            mhalf, mfull, isclp)
c
c     test Voronoi vertices
c
      write(*,*)' '
      write(*,*)'Testing circumcenters/Voronoi vertices ...'

      do 1060 i = 1, ivnxt
         if(ifl(i).eq.0) go to 1060
         go to 1070
 1060 continue
      write(*,*)' '
      write(*,*)'All Voronoi vertices marked as too large.'
      go to 1090
 1070 continue
      ng = 0
      derro = 0.0d0
      larg = 0
      do 1080 i = 1, ivnxt
         if(ifl(i).eq.0) then
            larg = larg+1
            go to 1080
         endif
         ifir = icon(4,i)
         isec = icon(5,i)
         ithi = icon(6,i)
         if(ifir.le.ng .or. isec.le.ng .or. ithi.le.ng) stop 150
c
         xcor = x(ifir)-vx(i)
         ycor = y(ifir)-vy(i)
         dnom=dmax1(dabs(xcor),dabs(ycor))
         if(dnom.gt.1.0d0) then
            xcor = xcor/dnom
            ycor = ycor/dnom
            dist1 = dnom*dsqrt(dabs(xcor**2 + ycor**2))
         else
            dist1 = dsqrt(dabs(xcor**2 + ycor**2))
         endif
c
         xcor = x(isec)-vx(i)
         ycor = y(isec)-vy(i)
         dnom=dmax1(dabs(xcor),dabs(ycor))
         if(dnom.gt.1.0d0) then
            xcor = xcor/dnom
            ycor = ycor/dnom
            dist2 = dnom*dsqrt(dabs(xcor**2 + ycor**2))
         else
            dist2 = dsqrt(dabs(xcor**2 + ycor**2))
         endif
c
         xcor = x(ithi)-vx(i)
         ycor = y(ithi)-vy(i)
         dnom=dmax1(dabs(xcor),dabs(ycor))
         if(dnom.gt.1.0d0) then
            xcor = xcor/dnom
            ycor = ycor/dnom
            dist3 = dnom*dsqrt(dabs(xcor**2 + ycor**2))
         else
            dist3 = dsqrt(dabs(xcor**2 + ycor**2))
         endif
c
         dista = dmax1(dist1,dist2,dist3)
         disti = dmin1(dist1,dist2,dist3)
         diffd = dista-disti
         if(diffd.gt.derro) derro = diffd
 1080 continue
      write(*,*)' '
      write(*,*)'Number of Voronoi vertices marked as too ',
     *          'large = ',larg
      write(*,*)' '
      write(*,*)'Voronoi distance error = ',derro
      write(*,*)'(This number should be close to zero).'
 1090 continue

 1500 continue
c
c*********************************************************************
c
c     test consistency of triangulation
c
      write(*,*)' '
      write(*,*)'Testing consistency of final triangulation ... '
      call contst(icon, is, id, nv, ivnxt, nvmz)
c
c*********************************************************************
c
c     test circle criterion
c
      write(*,*)' '
      write(*,*)'Testing circle criterion ... '
      call circri(x, y, ix, iy, ix2, iy2, icon, mhalf, mfull,
     *            isclp, ivnxt, epf)
c
c*********************************************************************
c
c     test orientation of triangles
c
      write(*,*)' '
      write(*,*)'Testing orientation of triangles ... '
      call orient(x, y, ix, iy, ix2, iy2, icon, ivnxt, mhalf, mfull,
     *            isclp, epf)
c
c*********************************************************************
c
c     test convexity of union of triangles
c
      write(*,*)' '
      write(*,*)'Testing convexity of union of triangles ... '
      call convex(x, y, ix, iy, ix2, iy2, icon, is, nv, ivnxt,
     *             nvmz, mhalf, mfull, isclp, epf)

c
c*********************************************************************
c
c     compress arrays icon, vx, vy, ifl, so that only Voronoi vertices
c     inside or on boundary of rectangle are kept
c
      if(rtflg) then
         nv = n
         call cmprs2(is, icon, id, vx, vy, ifl, nr, nv, ivnxt,
     *               nmax, nvmax)
      endif

      if(.not.toflg) go to 2150
c
c**********************************************************************
c
c     save information about final Delaunay triangulation
c
      write(*,*)' '
      write(*,*)'Saving information about final Delaunay ',
     *          'triangulation ... '

c     open(unit=16,file=toinf,status='new')
      open(unit=16,file=toinf)
      write(16,*) nv, ivnxt, nvmz, icsfig, lflag, lflug
      write(16,170)((icon(i,j),i=1,6),j=1,ivnxt)
      write(16,180)(is(i),i=1,nv)
      close(unit=16)

 2150 continue
      if(.not.trflg) go to 2200
c
c*********************************************************************
c
c     save vertices of triangles by themselves
c
      write(*,*)' '
      write(*,*)'Saving vertices of triangles by themselves ... '

c     open(unit=20,file=trinf,status='new')
      open(unit=20,file=trinf)
      write(20,2190)((icon(i,j),i=4,6),j=1,ivnxt)
      close(unit=20)
 2190 format(3i10)

 2200 continue
      if(.not.voflg) go to 2300
c
c**********************************************************************
c
c     save set of points (x-,y-coordinates) that are possibly
c     vertices in final triangulation
c
      write(*,*)' '
      write(*,*)'Saving points (x-,y-coordinates) that are possibly ',
     *          'vertices in'
      write(*,*)'final triangulation ... '
c     open(unit=17,file=voinf,status='new')
      open(unit=17,file=voinf)
      do 2250 i=1,nv
         write(17,*)x(i),y(i)
 2250 continue
      close(17)

 2300 continue
      if(.not.ccflg) go to 2380
c
c********************************************************************
c
c     save circumcenters (x-,y-coordinates) of triangles,
c     i. e., Voronoi vertices, in final triangulation
c
      write(*,*)' '
      write(*,*)'Saving circumcenters/Voronoi vertices ... '

c     open(unit=18,file=ccinf,status='new')
      open(unit=18,file=ccinf)
      write(18,*) nv, ivnxt, lflag, lflug
      write(18,2350) (ifl(i),i=1,ivnxt)
      do 2340 i=1,ivnxt
         write(18,*) vx(i),vy(i)
 2340 continue
 2350 format(40(1x,i1))
      close(unit=18)

 2380 continue
c********************************************************************
c
c     computations are completed
c
      write(*,*)' '
      write(*,*)'Computations completed.'
      write(*,*)' '
      write(*,*)'Number of vertices (points) in final ',
     *          'triangulation: ',nv
      write(*,*)'Number of triangles: ',ivnxt
      write(*,*)' '
c
c**********************************************************************
c
c     the following write's are for testing purposes only
c     they should not be used for big runs
c
c     write(*,*)' '
c     write(*,*)'Number of points: '
c     write(*,*) nv
c     write(*,*)'Triangle pointers, x and y coordinates for points: '
c     do 3100 i = 1, nv
c         write(*,*)' pnt: ',i,' pnter: ',is(i),
c    *               ' x: ',x(i),' y: ',y(i)
c3100 continue
c     write(*,*)'Number of triangles: '
c     write(*,*) ivnxt
c     write(*,*)'Triangles: '
c     do 3200 i = 1, ivnxt
c         write(*,3400) i, (icon(j,i),j=1,6)
c3200 continue
c3400 format('tri: ',i10,' ic: ',6i10)
c
      stop
      end
*TRANSF
c
c     subroutine transf to -
c
c     transform double precision coordinates into integer decomposition
c 
      subroutine transf(x, y, ix, iy, ix2, iy2, xmin, ymin, xmax, ymax,
     *                  nv, icsfig, mhalf, mfull, isclp, r215, deps,
     *                  dscle, dfull, dfill, epf)
c
      double precision x(*), y(*)
      integer ix(*), iy(*), ix2(*), iy2(*)
      integer nv, icsfig
      integer mhalf, mfull, ibsfig, itsfig
      double precision epf, wbig
      integer isclp(*)
      double precision xmin, xmax, ymin, ymax
      double precision dscle, dfull, dfill, decml, deps, r215
      integer i, isclu, isgcl
c
c     initialize word lengths
c
      mhalf=32768
      mfull=1073741824
      if(mfull.ne.mhalf**2) stop 200
      deps = dble(0.9)
      r215 = dble(mhalf)
c
c     test # of significant figures of nondecimal part of coordinates
c
      wbig = 0.0d0
      if(wbig .lt. dabs(xmax)) wbig = dabs(xmax)
      if(wbig .lt. dabs(xmin)) wbig = dabs(xmin)
      if(wbig .lt. dabs(ymax)) wbig = dabs(ymax)
      if(wbig .lt. dabs(ymin)) wbig = dabs(ymin)
      wbig = wbig + epf
c     WRITE(*,*)'COORDINATES WBIG=',WBIG
      ibsfig = 0
  180 continue
      ibsfig = ibsfig+1
      wbig = wbig/10.0d0
      if(wbig .ge. 1.0d0) go to 180
      if(ibsfig.gt.9) then
         write(*,*)' '
         write(*,*)'Number of significant figures of largest ',
     *             'nondecimal part of'
         write(*,*)'a point coordinate appears to be greater than 9.'
         write(*,*)'Program is terminated.'
         write(*,*)' '
         stop 210
      endif
      itsfig = ibsfig + icsfig
c     WRITE(*,*)'ITSFIG=',ITSFIG,' IBSFIG=',IBSFIG,' ICSFIG=',ICSFIG
      if(itsfig.gt.14) then
         write(*,*)' '
         write(*,*)'For this execution of the program the largest ',
     *             'total number of'
         write(*,*)'significant figures ',
     *             'that a coordinate requires appears to be ',itsfig
         write(*,*)'Since the maximum allowed is 14, the number of ',
     *             'significant figures'
         write(*,*)'of the decimal part of the point coordinates for ',
     *             'this run is '
         write(*,*)'decreased accordingly.'
         icsfig = 14 - ibsfig
         write(*,*)' '
      endif
c
c     test number of significant figures of decimal part of coordinates
c
      if(icsfig.lt.0 .or. icsfig.gt.9) stop 220
      isclu = 1
      dscle = 1.0d0
      if(icsfig.eq.0) go to 220
      do 210 i = 1, icsfig
         isclu = 10*isclu
         dscle = 10.0d0*dscle
  210 continue
  220 continue
      if(iabs(isclu).ge.mfull) stop 230
      call decomp(isclp, isgcl, isclu, mhalf)
      if(isgcl.ne.1) stop 240
c
c     test lengths of x, y-coordinates, shift and make them integers
c
      dfull = dble(mfull)
      if(dscle.lt.deps) stop 250
      dfill=dfull/dscle
      do 235 i = 1, nv
         ix2(i) = 0
         iy2(i) = 0
         if(dabs(x(i)).lt.dfill) then
            ix(i) = idnint(dscle*x(i))
            if(iabs(ix(i)).lt.mfull) then
               x(i) = dble(ix(i))/dscle
               go to 225
            endif
         endif
         if(dabs(x(i)).ge.dfull) stop 260
         ix(i) = idint(x(i))
         if(iabs(ix(i)).ge.mfull) stop 270
         decml = (x(i) - dint(x(i)))*dscle
         if(dabs(decml).ge.dfull) stop 280
         ix2(i) = idnint(decml)
         if(iabs(ix2(i)).ge.mfull) stop 290
         if(iabs(ix2(i)).eq.0) then
            x(i) = dble(ix(i))
            ix2(i) = mfull
         else
            x(i) = dble(ix(i)) + (dble(ix2(i))/dscle)
         endif
  225    continue
         if(dabs(y(i)).lt.dfill) then
            iy(i) = idnint(dscle*y(i))
            if(iabs(iy(i)).lt.mfull) then
               y(i) = dble(iy(i))/dscle
               go to 235
            endif
         endif
         if(dabs(y(i)).ge.dfull) stop 300
         iy(i) = idint(y(i))
         if(iabs(iy(i)).ge.mfull) stop 310
         decml = (y(i) - dint(y(i)))*dscle
         if(dabs(decml).ge.dfull) stop 320
         iy2(i) = idnint(decml)
         if(iabs(iy2(i)).ge.mfull) stop 330
         if(iabs(iy2(i)).eq.0) then
            y(i) = dble(iy(i))
            iy2(i) = mfull
         else
            y(i) = dble(iy(i)) + (dble(iy2(i))/dscle)
         endif
  235 continue
c
      return
      end
*PNTDEC
c
c     subroutine pntdec -
c
c     to decompose double precision coordinates of one point
c     into integers
c
c
      subroutine pntdec(xi, yi, x, y, x2, y2, n, mfull, icsfig,
     *                  dscle, dfull, dfill, epf)
c
      double precision xi(*), yi(*)
      integer x(*), y(*), x2(*), y2(*)
      integer n, mfull, icsfig, ibsfig, itsfig
      double precision xdes, ydes, wbig, epf
      double precision dscle, dfull, dfill, decml
c
c     test # of significant figures of nondecimal part of coordinates
c
      xdes = xi(n)
      ydes = yi(n)
      wbig = 0.0d0
      if(wbig .lt. dabs(xdes)) wbig = dabs(xdes)
      if(wbig .lt. dabs(xdes)) wbig = dabs(xdes)
      if(wbig .lt. dabs(ydes)) wbig = dabs(ydes)
      if(wbig .lt. dabs(ydes)) wbig = dabs(ydes)
      wbig = wbig + epf
c     WRITE(*,*)'COORDINATES WBIG=',WBIG
      ibsfig = 0
  180 continue
      ibsfig = ibsfig+1
      wbig = wbig/10.0d0
      if(wbig .ge. 1.0d0) go to 180
      if(ibsfig.gt.9) then
         write(*,*)' '
         write(*,*)'Number of significant figures of nondecimal',
     *             'part of a'
         write(*,*)'new point coordinate appears to be greater than 9.'
         write(*,*)'Program is terminated.'
         write(*,*)' '
         stop 340
      endif
      itsfig = ibsfig + icsfig
c     WRITE(*,*)'ITSFIG=',ITSFIG,' IBSFIG=',IBSFIG,' ICSFIG=',ICSFIG
      if(itsfig.gt.14) then
         write(*,*)' '
         write(*,*)'Total number of significant figures of a new ',
     *             'point coordinate'
         write(*,*)'appears to be greater than 14.'
         write(*,*)'Program is terminated.'
         write(*,*)' '
         stop 350
      endif
c
      x2(n) = 0
      y2(n) = 0
      if(dabs(xi(n)).lt.dfill) then
         x(n) = idnint(dscle*xi(n))
         if(iabs(x(n)).lt.mfull) then
            xi(n) = dble(x(n))/dscle
            go to 1110
         endif
      endif
      if(dabs(xi(n)).ge.dfull) stop 360
      x(n) = idint(xi(n))
      if(iabs(x(n)).ge.mfull) stop 370
      decml = (xi(n) - dint(xi(n)))*dscle
      if(dabs(decml).ge.dfull) stop 380
      x2(n) = idnint(decml)
      if(iabs(x2(n)).ge.mfull) stop 390
      if(iabs(x2(n)).eq.0) then
         xi(n) = dble(x(n))
         x2(n) = mfull
      else
         xi(n) = dble(x(n)) + (dble(x2(n))/dscle)
      endif
 1110 continue
      if(dabs(yi(n)).lt.dfill) then
         y(n) = idnint(dscle*yi(n))
         if(iabs(y(n)).lt.mfull) then
            yi(n) = dble(y(n))/dscle
            go to 1120
         endif
      endif
      if(dabs(yi(n)).ge.dfull) stop 400
      y(n) = idint(yi(n))
      if(iabs(y(n)).ge.mfull) stop 410
      decml = (yi(n) - dint(yi(n)))*dscle
      if(dabs(decml).ge.dfull) stop 420
      y2(n) = idnint(decml)
      if(iabs(y2(n)).ge.mfull) stop 430
      if(iabs(y2(n)).eq.0) then
         yi(n) = dble(y(n))
         y2(n) = mfull
      else
         yi(n) = dble(y(n)) + (dble(y2(n))/dscle)
      endif
 1120 continue
c
      return
      end
*TRICMP
c
c     driver subroutine tricmp to -
c
c     compute a Delaunay triangulation from scratch
c
      subroutine tricmp(x, y, ix, iy, ix2, iy2, is, icon, iave, iaze,
     *                  epf, mhalf, mfull, isclp, nvmz, nv, iabe, iase,
     *                  nmax, nvmax, ivnxt, idupl, inpro, r215, deps,
     *                  ileft)
c
      double precision x(*), y(*)
      integer ix(*), iy(*), ix2(*), iy2(*)
      integer is(*), icon(6,*), iave(*), iaze(*)
      double precision epf, r215, deps
      integer mhalf, mfull, nvmz, nv, iabe, iase, nmax, nvmax, ivnxt
      integer isclp(*)
c
      integer nkmax
      parameter (nkmax = 30)
      integer iax(nkmax), iay(nkmax), isgax, isgay, ikax, ikay
      integer iox(nkmax), ioy(nkmax), isgox, isgoy, ikox, ikoy
      integer isga, ifreq
      double precision xtemp, ytemp, delx, dely
      double precision xmax, xmin, dotx, dnux, dot
      integer i, i1, i2, i3, it, ileft, iright, insrt, isst
      integer idupl, inpro, isite, itype, ipoint
c
c     obtain initial triangle
c
      do 95 i = 1, nv
         if(is(i).eq.0) go to 98
   95 continue
      stop 440
c
   98 continue
      i1 = i
      i2 = i
      xmin = x(i1)
      xmax = x(i2)
      do 100 i = 1, nv
         if(is(i).ne.0) go to 100
         if(x(i) .lt. xmin) then
            i1 = i
            xmin = x(i1)
         endif
         if(x(i) .gt. xmax) then
            i2 = i
            xmax = x(i2)
         endif
  100 continue
      if(i1.eq.i2) stop 450
c
      call prpvec(ix, iy, ix2, iy2, i1, i2, mhalf, mfull, isclp,
     *            iax, isgax, ikax, iay, isgay, ikay)
      if(isgay.le.0) stop 460
      call doubnm(iax, isgax, ikax, r215, delx)
      call doubnm(iay, isgay, ikay, r215, dely)
      dnux = dmax1(dabs(delx),dabs(dely))
      if(dnux.lt.deps) stop 470
      delx = delx/dnux
      dely = dely/dnux
c
      i3 = 0
      dotx = 0.0d0
      do 350 i = 1, nv
         if(i.eq.i1 .or. i.eq.i2 .or. is(i).ne.0) go to 350
         call vectrd(ix, iy, ix2, iy2, i1, i, mhalf, mfull, isclp,
     *               iox, isgox, ikox, ioy, isgoy, ikoy)
         call doubnm(iox, isgox, ikox, r215, xtemp)
         call doubnm(ioy, isgoy, ikoy, r215, ytemp)
         dot = delx*xtemp + dely*ytemp
c        WRITE(*,*)'I=',I,' DOT =',DOT
         if(dabs(dot).gt.dotx) then
            dotx = dabs(dot)
            i3 = i
         endif
  350 continue
c
      if(i3 .eq. 0) then
         write(*,*)' '
         write(*,*)'Points appear colinear or close to colinear'
         write(*,*)'Program terminated'
         write(*,*)' '
         stop 480
      endif
c
      call distsn(ix, iy, ix2, iy2, i1, i3, mhalf, mfull, isclp,
     *            iax, isgax, ikax, iay, isgay, ikay, isga)
c
      if(isga .eq. 0) stop 490
      if(isga .gt. 0) go to 400
      it = i2
      i2 = i3
      i3 = it
c
  400 continue
c     write(*,*)' '
c     write(*,*)'Initial triangle has vertices:',i1,i2,i3
      ivnxt = 1
      if(ivnxt.gt.nvmax) stop 500
c
      icon(1,1) = 0
      icon(2,1) = 0
      icon(3,1) = 0
      icon(4,1) = i1
      icon(5,1) = i2
      icon(6,1) = i3
c
      is(i1) = 1
      is(i2) = 1
      is(i3) = 1
c
c     insert other points
c
      iabe=1
      iave(1)=ivnxt+1
      ileft=i1
      insrt=1
      idupl=0
      inpro=0
      if(nv.eq.3) go to 480
c
      ifreq = 10000
      if(ifreq.le.0) stop 503
      if(ifreq.le.nv) write(*,*)' '
      do 450 isst = 1, nv
         isite = isst
         if(is(isite).eq.nvmz) inpro=inpro+1
         if(is(isite).ne.0) go to 450
         ipoint = isite
c        WRITE(*,*)'IPOINT=',IPOINT
         call pntins(x, y, ix, iy, ix2, iy2, is, icon, iave, iaze,
     *               epf, mhalf, mfull, isclp, nvmz, ileft, iright,
     *               nv, iabe, iase, isite, nmax, nvmax, ivnxt,
     *               insrt, itype, ipoint)
         if(itype.eq.1.or.itype.eq.-1) then
            is(isite) = -iright
            idupl=idupl+1
         endif
c
         ileft = iright
         if(isite.le.(isite/ifreq)*ifreq)
     *      write(*,*)'Number of points inserted =',isite
  450 continue
c
  480 continue
      return
      end
*TRIADD
c
c     driver subroutine triadd to -
c
c     insert additional points into current Delaunay triangulation
c
      subroutine triadd(x, y, ix, iy, ix2, iy2, is, icon, ip, iave,
     *                 iaze, epf, icsfig, mhalf, mfull, isclp, nvmz,
     *                 n, nv, np, iabe, iase, nmax, nvmax, ivnxt,
     *                 idupl, inpro, ileft, dscle, dfull, dfill)
c
      double precision x(*), y(*)
      integer ix(*), iy(*), ix2(*), iy2(*)
      integer is(*), icon(6,*), ip(*), iave(*), iaze(*)
      double precision epf, dscle, dfull, dfill
      integer icsfig, mhalf, mfull, isclp(*)
      integer nvmz, n, nv, np, iabe, iase, nmax, nvmax, ivnxt
c
      integer ileft, iright, insrt, isst, ifreq
      integer idupl, inpro, isite, itype, i, ipoint
c
      if(np.gt.0)then
         do 575 i=1,np
            if(ip(i).lt.n+1.or.ip(i).gt.nv) go to 575
            is(ip(i))=nvmz
  575    continue
      endif

      insrt=1
      idupl=0
      inpro=0
c
      ifreq = 10000
      if(ifreq.le.0) stop 505
      if(n.lt.(nv/ifreq)*ifreq) write(*,*)' '
      do 580 isst = n+1, nv
         isite = isst
         if(is(isite).eq.nvmz) inpro=inpro+1
         if(is(isite).ne.0) go to 580
         ipoint = isite
         call pntdec(x, y, ix, iy, ix2, iy2, ipoint, mfull, icsfig,
     *               dscle, dfull, dfill, epf)
         call pntins(x, y, ix, iy, ix2, iy2, is, icon, iave, iaze,
     *               epf, mhalf, mfull, isclp, nvmz, ileft, iright,
     *               nv, iabe, iase, isite, nmax, nvmax, ivnxt,
     *               insrt, itype, ipoint)
         if(itype.eq.1.or.itype.eq.-1) then
            is(isite) = -iright
            idupl=idupl+1
         endif
c
         ileft = iright
         if(isite.le.(isite/ifreq)*ifreq)
     *      write(*,*)'Number of points inserted =',isite
  580 continue
      return
      end
*LININS
c
c     driver subroutine linins to -
c
c     insert line segments into Delaunay triangulation and maintain
c     Delaunay property
c
      subroutine linins(x, y, xrg, yrg, xds, yds, ix, iy, ix2, iy2, is,
     *                  icon, iave, iaze, ia, ifun, irun, isun, izun,
     *                  igun, iwun, epf, icsfig, mhalf, mfull, isclp,
     *                  nvmz, nv, nl, iabe, iase, nmax, nvmax, nfmax,
     *                  nhmax, ivnxt, nvmp, nvmn, nls, nlp, ileft,
     *                  dscle, dfull, dfill)
c
      double precision x(*), y(*), xrg(*), yrg(*), xds(*), yds(*)
      integer ix(*), iy(*), ix2(*), iy2(*)
      integer is(*), icon(6,*), iave(*), iaze(*), ia(*)
      integer ifun(*), irun(*), isun(*), izun(*)
      integer igun(*), iwun(*)
      double precision epf, dscle, dfull, dfill
      integer icsfig, mhalf, mfull, isclp(*)
      integer nvmz, nv, nl, iabe, iase, nmax, nvmax, ivnxt
      integer nfmax, nhmax, nvmp, nvmn, ileft, iright
      integer nls, nlp, ifreq
c
      integer i, insrt, isite, itype, ipoint
c
      insrt=1
      nls=0
      nlp=0
      ifreq = 1000
      if(ifreq.le.0) stop 507
      if(ifreq.le.nl) write(*,*)' '
      do 600 i = 1, nl
         isite=0
         ipoint = nv + 1
         if(ipoint.gt.nmax) stop 510
         x(ipoint) = xrg(i)
         y(ipoint) = yrg(i)
         call pntdec(x, y, ix, iy, ix2, iy2, ipoint, mfull, icsfig,
     *               dscle, dfull, dfill, epf)
         call pntins(x, y, ix, iy, ix2, iy2, is, icon, iave, iaze,
     *               epf, mhalf, mfull, isclp, nvmz, ileft, iright,
     *               nv, iabe, iase, isite, nmax, nvmax, ivnxt,
     *               insrt, itype, ipoint)
c
         ileft = iright
         isite=0
         ipoint = nv + 1
         if(ipoint.gt.nmax) stop 520
         x(ipoint) = xds(i)
         y(ipoint) = yds(i)
         call pntdec(x, y, ix, iy, ix2, iy2, ipoint, mfull, icsfig,
     *               dscle, dfull, dfill, epf)
         call pntins(x, y, ix, iy, ix2, iy2, is, icon, iave, iaze,
     *               epf, mhalf, mfull, isclp, nvmz, ileft, iright,
     *               nv, iabe, iase, isite, nmax, nvmax, ivnxt,
     *               insrt, itype, ipoint)
c
         call reconc(x, y, ix, iy, ix2, iy2, is, icon, ifun, irun,
     *               isun, izun, igun, iwun, ia, epf, mhalf, mfull,
     *               isclp, nfmax, nhmax, nvmz, nvmp, nvmn,
     *               ileft, iright, nv)
c
         if(ileft .ne. iright) then
            ileft = iright
            nls=nls+1
         else
            nlp=nlp+1
         endif
         if(i.le.(i/ifreq)*ifreq)
     *      write(*,*)'Number of line segments inserted =',i
  600 continue
c
      return
      end
*PNTDEL
c
c     driver subroutine pntdel to -
c
c     delete points or vertices in current Delaunay triangulation
c
      subroutine pntdel(x, y, xdl, ydl, ix, iy, ix2, iy2, is, icon,
     *                  iave, iaze, ia, ifun, irun, isun, izun, igun,
     *                  iwun, epf, icsfig, mhalf, mfull, isclp, nvmz,
     *                  nv, nd, iabe, iase, nmax, nvmax, nfmax, nhmax,
     *                  njmax, nkmax, ivnxt, nvmp, nvmn, ndl, nde, ndv,
     *                  nds, ileft, dscle, dfull, dfill)
c
      double precision x(*), y(*), xdl(*), ydl(*)
      integer ix(*), iy(*), ix2(*), iy2(*)
      integer is(*), icon(6,*), iave(*), iaze(*), ia(*)
      integer ifun(*), irun(*), isun(*), izun(*)
      integer igun(*), iwun(*)
      double precision epf, dscle, dfull, dfill
      integer icsfig, mhalf, mfull, isclp(*)
      integer nvmz, nv, nd, iabe, iase
      integer nmax, nvmax, nfmax, nhmax, njmax, nkmax
      integer ivnxt, nvmp, nvmn, ndl, nde, ndv, nds
c
      integer ileft, iright, insrt, ifreq
      integer isite, itype, i, isucs, isnew, ipoint
c
      insrt = 0
      isite = 0
      ndl=0
      nde=0
      ndv=0
      nds=0
      ipoint = nv + 1
      if(ipoint.gt.nmax) stop 530
c
      write(*,*)' '
      ifreq = 1000
      if(ifreq.le.0) stop 535
      if(ifreq.le.nd) write(*,*)' '
      do 720 i = 1, nd
         x(ipoint) = xdl(i)
         y(ipoint) = ydl(i)
         call pntdec(x, y, ix, iy, ix2, iy2, ipoint, mfull, icsfig,
     *               dscle, dfull, dfill, epf)
         isucs = -1
         call pntins(x, y, ix, iy, ix2, iy2, is, icon, iave, iaze,
     *               epf, mhalf, mfull, isclp, nvmz, ileft, iright,
     *               nv, iabe, iase, isite, nmax, nvmax, ivnxt,
     *               insrt, itype, ipoint)
c
         if(itype .eq. 1 .or. itype .eq. -1) then
            call deleli(x, y, ix, iy, ix2, iy2, is, icon, ifun, irun,
     *                  isun, izun, igun, iwun, ia, iave, iaze, epf,
     *                  mhalf, mfull, isclp, nfmax, nhmax, njmax,
     *                  nkmax, nvmz, nvmp, nvmn, iright, nv, isucs,
     *                  isnew, iabe, iase)
         endif
c
         if(isucs .eq. 1) then
            ndl=ndl+1
            ileft = isnew
         elseif(isucs .eq. 0) then
            nde=nde+1
         elseif(isucs .eq. -1) then
            ndv=ndv+1
         else
            nds=nds+1
         endif
         if(i.le.(i/ifreq)*ifreq)
     *      write(*,*)'Number of points deleted =',i
  720 continue
c
      return
      end
*VORCMP
c
c     This subroutine computes Voronoi vertices.
c 
c
      subroutine vorcmp(ix, iy, ix2, iy2, xp, yp, icon, ifl, nt,
     *                  dscle, mhalf, mfull, isclp)
c
      integer ix(*), iy(*), ix2(*), iy2(*)
      double precision xp(*), yp(*)
      integer icon(6,*), ifl(*)
      integer nt, mhalf, mfull, nkmax, isclp(*)
      parameter(nkmax=30)
      integer io(nkmax), inx(nkmax), iny(nkmax)
      double precision r215, dnom, xnum, ynum
      double precision deps, dscle, dscl2, dfull
c
      integer ng, i, ifir, isec, ithi, isgo, iko
      integer isgnx, isgny, iknx, ikny, ifreq
c
c     compute vertices
c
      ng = 0
      dfull = dble(mfull)
      r215 = dble(mhalf)
      deps = dble(0.9)
      if(dscle.lt.deps) stop 540
      dscl2 = 2.0d0*dscle
      ifreq = 10000
      if(ifreq.le.0) stop 545
      if(ifreq.le.nt) write(*,*)' '
      do 1000 i = 1, nt
         if(i.le.(i/ifreq)*ifreq)
     *   write(*,*)'Number of computed vertices =',i
         ifir = icon(4,i)
         isec = icon(5,i)
         ithi = icon(6,i)
         if(ifir.le.ng .or. isec.le.ng .or. ithi.le.ng) stop 550
         call numden(ix, iy, ix2, iy2, ifir, isec, ithi, mhalf, mfull,
     *               isclp, io, isgo, iko, inx, iny, isgnx, isgny,
     *               iknx, ikny)
         if(isgo.le.0) stop 560
         call doubnm(io, isgo, iko, r215, dnom)
         if(dnom.lt.deps) stop 570
         call doubnm(inx, isgnx, iknx, r215, xnum)
         call doubnm(iny, isgny, ikny, r215, ynum)
         xp(i) = (xnum/dnom)/dscl2
         yp(i) = (ynum/dnom)/dscl2
         ifl(i) = 1
 1000 continue
c
c     identify vertices whose sizes are too big
c
      do 1600 i = 1, nt
         dnom = dmax1(dabs(xp(i)),dabs(yp(i)))
         if(dnom.gt.dfull) ifl(i) = 0
 1600 continue
c
      return
      end
*VCLIDN
c
c     This subroutine identifies Voronoi cells, parts of which
c     intersect the exterior of a rectangle, parts that will
c     be clipped during insertion of edges of rectangle into
c     Voronoi diagram.
c
      subroutine vclidn(vx, vy, x, y, icon, is, nr, nv, ivnxt, nmax,
     *                  xne, yne, xnw, ynw, xsw, ysw, xse, yse)
c
      double precision vx(*), vy(*), x(*), y(*)
      integer icon(6,*), is(*)
      integer nr, nv, ivnxt, nmax 
      double precision xne, yne, xnw, ynw, xsw, ysw, xse, yse
c
      double precision dne, dnw, dsw, dse, dst
      integer ine, inw, isw, ise, icr, ikr, isr, ntot
      integer islst, iscur, isnxt, i, ileft
c
c     identify points closest to corners of rectangle, i. e., Voronoi
c     cells that contain corners
c
      do 1060 i = 1, nv
         if(is(i).le.0) go to 1060
         ine = i
         inw = i
         isw = i
         ise = i
         dne = (x(i)-xne)**2 + (y(i)-yne)**2
         dnw = (x(i)-xnw)**2 + (y(i)-ynw)**2
         dsw = (x(i)-xsw)**2 + (y(i)-ysw)**2
         dse = (x(i)-xse)**2 + (y(i)-yse)**2
         ileft = i+1
         go to 1070
 1060 continue
      stop 580
 1070 continue
c
      if(ileft.gt.nv) stop 590
      do 1080 i = ileft, nv
         if(is(i).le.0) go to 1080
         dst = (x(i)-xne)**2 + (y(i)-yne)**2
         if(dst.lt.dne) then
            dne = dst
            ine = i
         endif
         dst = (x(i)-xnw)**2 + (y(i)-ynw)**2
         if(dst.lt.dnw) then
            dnw = dst
            inw = i
         endif
         dst = (x(i)-xsw)**2 + (y(i)-ysw)**2
         if(dst.lt.dsw) then
            dsw = dst
            isw = i
         endif
         dst = (x(i)-xse)**2 + (y(i)-yse)**2
         if(dst.lt.dse) then
            dse = dst
            ise = i
         endif
 1080 continue
c
c     WRITE(*,*)'INE =',INE,' DNE =',DNE
c     WRITE(*,*)'INW =',INW,' DNW =',DNW
c     WRITE(*,*)'ISW =',ISW,' DSW =',DSW
c     WRITE(*,*)'ISE =',ISE,' DSE =',DSE
c     WRITE(*,*)' '
c     WRITE(*,*)'INE=',INE,' ISE=',ISE,' ISW=',ISW,' INW=',INW
c
c     go around edges of rectangle in clockwise direction
c     starting at NE corner and ending at the same corner
c
c     process NE corner and then NE-SE edge
c
      nr = nv + 1
      if(nr .gt. nmax) stop 600
      x(nr) = 2.0d0*xne - x(ine)
      y(nr) = y(ine)
      ntot = 1
c     WRITE(*,*)' '
c     WRITE(*,*)'INE=',INE
      if(ine.eq.ise) go to 1900
      isr  = ine
      iscur= is(ine)
      if(iscur.le.0 .or. iscur.gt.ivnxt) stop 610
      if(icon(4,iscur).eq.ine) then
         islst = icon(3,iscur)
         isnxt = icon(2,iscur)
         ikr = icon(5,iscur)
         icr = icon(6,iscur)
      elseif(icon(5,iscur).eq.ine) then
         islst = icon(1,iscur)
         isnxt = icon(3,iscur)
         ikr = icon(6,iscur)
         icr = icon(4,iscur)
      elseif(icon(6,iscur).eq.ine) then
         islst = icon(2,iscur)
         isnxt = icon(1,iscur)
         ikr = icon(4,iscur)
         icr = icon(5,iscur)
      else
         stop 620
      endif
      if(islst.ne.0) islst = iscur
c
c     test for vertex in initial Voronoi cell left of NE-SE line
c
      if(vx(iscur).le.xne) go to 1300
c
c     test for bounded initial Voronoi cell
c
      if(islst.ne.0) go to 1100
c
c     test for unbounded edge moving into rectangle
c
      if(y(ikr).ge.y(ine)) go to 1100
      nr = nr + 1
      if(nr .gt. nmax) stop 630
      x(nr) = 2.0d0*xne - x(ikr)
      y(nr) = y(ikr)
      ntot = 2
c     WRITE(*,*)'IKR 1 =',IKR
      if(ikr.eq.ise) go to 1900
      isr  = ikr
      islst = iscur
      go to 1500
c
c     find vertex of Voronoi cell left of NE-SE line
c     (in counterclockwise direction around isr)
c
 1100 continue
      if(isnxt.eq.0 .or. isnxt.eq.islst) stop 640
      iscur = isnxt
      if(vx(iscur).le.xne) then
         if(islst.ne.0) islst = iscur
         go to 1200
      endif
      if(icon(4,iscur).eq.isr) then
         isnxt = icon(2,iscur)
      elseif(icon(5,iscur).eq.isr) then
         isnxt = icon(3,iscur)
      elseif(icon(6,iscur).eq.isr) then
         isnxt = icon(1,iscur)
      else
         stop 650
      endif
      go to 1100
c
c     find vertex of Voronoi cell right of NE-SE line
c     (in counterclockwise direction around isr)
c
 1200 continue
      if(icon(4,iscur).eq.isr) then
         isnxt = icon(2,iscur)
         icr = icon(6,iscur)
      elseif(icon(5,iscur).eq.isr) then
         isnxt = icon(3,iscur)
         icr = icon(4,iscur)
      elseif(icon(6,iscur).eq.isr) then
         isnxt = icon(1,iscur)
         icr = icon(5,iscur)
      else
         stop 660
      endif
c
 1300 continue
      if(isnxt.eq.0) then
         if(islst.ne.0) stop 670
         go to 1400
      endif
      if(isnxt.eq.islst) stop 680
      iscur = isnxt
      if(vx(iscur).le.xne) go to 1200
 1400 continue
      nr = nr + 1
      if(nr .gt. nmax) stop 690
      x(nr) = 2.0d0*xne - x(icr)
      y(nr) = y(icr)
      ntot = ntot + 1
      if(ntot.gt.nv) stop 700
c     WRITE(*,*)'ICR=',ICR
      if(icr.eq.ise) go to 1900
      isr  = icr
      if(isnxt.eq.0) go to 1200
      islst = iscur
c
c     test vertices around isr in clockwise direction
c
 1500 continue
      if(icon(4,iscur).eq.isr) then
         isnxt = icon(3,iscur)
         ikr = icon(5,iscur)
      elseif(icon(5,iscur).eq.isr) then
         isnxt = icon(1,iscur)
         ikr = icon(6,iscur)
      elseif(icon(6,iscur).eq.isr) then
         isnxt = icon(2,iscur)
         ikr = icon(4,iscur)
      else
         stop 720
      endif
c
      if(isnxt.eq.0 .and. y(ikr).ge.y(isr)) go to 1800
      if(isnxt.eq.0) go to 1600
      if(isnxt.eq.islst) stop 740
      if(vx(isnxt).le.xne) go to 1600
      iscur = isnxt
      go to 1500
 1600 continue
      nr = nr + 1
      if(nr .gt. nmax) stop 750
      x(nr) = 2.0d0*xne - x(ikr)
      y(nr) = y(ikr)
      ntot = ntot + 1
      if(ntot.gt.nv) stop 760
c     WRITE(*,*)'IKR 2 =',IKR
      if(ikr.eq.ise) go to 1900
      isr  = ikr
      islst = iscur
      go to 1500
c
c     switch directions (to counterclockwise) around isr
c
 1800 continue
      if(icon(4,islst).eq.isr) then
         isnxt = icon(2,islst)
      elseif(icon(5,islst).eq.isr) then
         isnxt = icon(3,islst)
      elseif(icon(6,islst).eq.isr) then
         isnxt = icon(1,islst)
      else
         stop 770
      endif
      islst = 0
      go to 1100
c
c     process SE corner and then SE-SW edge
c
 1900 continue
      nr = nr + 1
      if(nr .gt. nmax) stop 800
      y(nr) = 2.0d0*yse - y(ise)
      x(nr) = x(ise)
      ntot = 1
c     WRITE(*,*)' '
c     WRITE(*,*)'ISE=',ISE
      if(ise.eq.isw) go to 2900
      isr  = ise
      iscur= is(ise)
      if(iscur.le.0 .or. iscur.gt.ivnxt) stop 810
      if(icon(4,iscur).eq.ise) then
         islst = icon(3,iscur)
         isnxt = icon(2,iscur)
         ikr = icon(5,iscur)
         icr = icon(6,iscur)
      elseif(icon(5,iscur).eq.ise) then
         islst = icon(1,iscur)
         isnxt = icon(3,iscur)
         ikr = icon(6,iscur)
         icr = icon(4,iscur)
      elseif(icon(6,iscur).eq.ise) then
         islst = icon(2,iscur)
         isnxt = icon(1,iscur)
         ikr = icon(4,iscur)
         icr = icon(5,iscur)
      else
         stop 820
      endif
      if(islst.ne.0) islst = iscur
c
c     test for vertex in initial Voronoi cell above SE-SW line
c
      if(vy(iscur).ge.yse) go to 2300
c
c     test for bounded initial Voronoi cell
c
      if(islst.ne.0) go to 2100
c
c     test for unbounded edge moving into rectangle
c
      if(x(ikr).ge.x(ise)) go to 2100
      nr = nr + 1
      if(nr .gt. nmax) stop 830
      y(nr) = 2.0d0*yse - y(ikr)
      x(nr) = x(ikr)
      ntot = 2
c     WRITE(*,*)'IKR 1 =',IKR
      if(ikr.eq.isw) go to 2900
      isr  = ikr
      islst = iscur
      go to 2500
c
c     find vertex of Voronoi cell above SE-SW line
c     (in counterclockwise direction around isr)
c
 2100 continue
      if(isnxt.eq.0 .or. isnxt.eq.islst) stop 840
      iscur = isnxt
      if(vy(iscur).ge.yse) then
         if(islst.ne.0) islst = iscur
         go to 2200
      endif
      if(icon(4,iscur).eq.isr) then
         isnxt = icon(2,iscur)
      elseif(icon(5,iscur).eq.isr) then
         isnxt = icon(3,iscur)
      elseif(icon(6,iscur).eq.isr) then
         isnxt = icon(1,iscur)
      else
         stop 850
      endif
      go to 2100
c
c     find vertex of Voronoi cell below SE-SW line
c     (in counterclockwise direction around isr)
c
 2200 continue
      if(icon(4,iscur).eq.isr) then
         isnxt = icon(2,iscur)
         icr = icon(6,iscur)
      elseif(icon(5,iscur).eq.isr) then
         isnxt = icon(3,iscur)
         icr = icon(4,iscur)
      elseif(icon(6,iscur).eq.isr) then
         isnxt = icon(1,iscur)
         icr = icon(5,iscur)
      else
         stop 860
      endif
c
 2300 continue
      if(isnxt.eq.0) then
         if(islst.ne.0) stop 870
         go to 2400
      endif
      if(isnxt.eq.islst) stop 880
      iscur = isnxt
      if(vy(iscur).ge.yse) go to 2200
 2400 continue
      nr = nr + 1
      if(nr .gt. nmax) stop 890
      y(nr) = 2.0d0*yse - y(icr)
      x(nr) = x(icr)
      ntot = ntot + 1
      if(ntot.gt.nv) stop 900
c     WRITE(*,*)'ICR=',ICR
      if(icr.eq.isw) go to 2900
      isr  = icr
      if(isnxt.eq.0) go to 2200
      islst = iscur
c
c     test vertices around isr in clockwise direction
c
 2500 continue
      if(icon(4,iscur).eq.isr) then
         isnxt = icon(3,iscur)
         ikr = icon(5,iscur)
      elseif(icon(5,iscur).eq.isr) then
         isnxt = icon(1,iscur)
         ikr = icon(6,iscur)
      elseif(icon(6,iscur).eq.isr) then
         isnxt = icon(2,iscur)
         ikr = icon(4,iscur)
      else
         stop 920
      endif
c
      if(isnxt.eq.0 .and. x(ikr).ge.x(isr)) go to 2800
      if(isnxt.eq.0) go to 2600
      if(isnxt.eq.islst) stop 940
      if(vy(isnxt).ge.yse) go to 2600
      iscur = isnxt
      go to 2500
 2600 continue
      nr = nr + 1
      if(nr .gt. nmax) stop 950
      y(nr) = 2.0d0*yse - y(ikr)
      x(nr) = x(ikr)
      ntot = ntot + 1
      if(ntot.gt.nv) stop 960
c     WRITE(*,*)'IKR 2 =',IKR
      if(ikr.eq.isw) go to 2900
      isr  = ikr
      islst = iscur
      go to 2500
c
c     switch directions (to counterclockwise) around isr
c
 2800 continue
      if(icon(4,islst).eq.isr) then
         isnxt = icon(2,islst)
      elseif(icon(5,islst).eq.isr) then
         isnxt = icon(3,islst)
      elseif(icon(6,islst).eq.isr) then
         isnxt = icon(1,islst)
      else
         stop 970
      endif
      islst = 0
      go to 2100
c
c     process SW corner and then SW-NW edge
c
 2900 continue
      nr = nr + 1
      if(nr .gt. nmax) stop 1000
      x(nr) = 2.0d0*xsw - x(isw)
      y(nr) = y(isw)
      ntot = 1
c     WRITE(*,*)' '
c     WRITE(*,*)'ISW=',ISW
      if(isw.eq.inw) go to 3900
      isr  = isw
      iscur= is(isw)
      if(iscur.le.0 .or. iscur.gt.ivnxt) stop 1010
      if(icon(4,iscur).eq.isw) then
         islst = icon(3,iscur)
         isnxt = icon(2,iscur)
         ikr = icon(5,iscur)
         icr = icon(6,iscur)
      elseif(icon(5,iscur).eq.isw) then
         islst = icon(1,iscur)
         isnxt = icon(3,iscur)
         ikr = icon(6,iscur)
         icr = icon(4,iscur)
      elseif(icon(6,iscur).eq.isw) then
         islst = icon(2,iscur)
         isnxt = icon(1,iscur)
         ikr = icon(4,iscur)
         icr = icon(5,iscur)
      else
         stop 1020
      endif
      if(islst.ne.0) islst = iscur
c
c     test for vertex in initial Voronoi cell right of SW-NW line
c
      if(vx(iscur).ge.xsw) go to 3300
c
c     test for bounded initial Voronoi cell
c
      if(islst.ne.0) go to 3100
c
c     test for unbounded edge moving into rectangle
c
      if(y(ikr).le.y(isw)) go to 3100
      nr = nr + 1
      if(nr .gt. nmax) stop 1030
      x(nr) = 2.0d0*xsw - x(ikr)
      y(nr) = y(ikr)
      ntot = 2
c     WRITE(*,*)'IKR 1 =',IKR
      if(ikr.eq.inw) go to 3900
      isr  = ikr
      islst = iscur
      go to 3500
c
c     find vertex of Voronoi cell right of SW-NW line
c     (in counterclockwise direction around isr)
c
 3100 continue
      if(isnxt.eq.0 .or. isnxt.eq.islst) stop 1040
      iscur = isnxt
      if(vx(iscur).ge.xsw) then
         if(islst.ne.0) islst = iscur
         go to 3200
      endif
      if(icon(4,iscur).eq.isr) then
         isnxt = icon(2,iscur)
      elseif(icon(5,iscur).eq.isr) then
         isnxt = icon(3,iscur)
      elseif(icon(6,iscur).eq.isr) then
         isnxt = icon(1,iscur)
      else
         stop 1050
      endif
      go to 3100
c
c     find vertex of Voronoi cell left of SW-NW line
c     (in counterclockwise direction around isr)
c
 3200 continue
      if(icon(4,iscur).eq.isr) then
         isnxt = icon(2,iscur)
         icr = icon(6,iscur)
      elseif(icon(5,iscur).eq.isr) then
         isnxt = icon(3,iscur)
         icr = icon(4,iscur)
      elseif(icon(6,iscur).eq.isr) then
         isnxt = icon(1,iscur)
         icr = icon(5,iscur)
      else
         stop 1060
      endif
c
 3300 continue
      if(isnxt.eq.0) then
         if(islst.ne.0) stop 1070
         go to 3400
      endif
      if(isnxt.eq.islst) stop 1080
      iscur = isnxt
      if(vx(iscur).ge.xsw) go to 3200
 3400 continue
      nr = nr + 1
      if(nr .gt. nmax) stop 1090
      x(nr) = 2.0d0*xsw - x(icr)
      y(nr) = y(icr)
      ntot = ntot + 1
      if(ntot.gt.nv) stop 1100
c     WRITE(*,*)'ICR=',ICR
      if(icr.eq.inw) go to 3900
      isr  = icr
      if(isnxt.eq.0) go to 3200
      islst = iscur
c
c     test vertices around isr in clockwise direction
c
 3500 continue
      if(icon(4,iscur).eq.isr) then
         isnxt = icon(3,iscur)
         ikr = icon(5,iscur)
      elseif(icon(5,iscur).eq.isr) then
         isnxt = icon(1,iscur)
         ikr = icon(6,iscur)
      elseif(icon(6,iscur).eq.isr) then
         isnxt = icon(2,iscur)
         ikr = icon(4,iscur)
      else
         stop 1120
      endif
c
      if(isnxt.eq.0 .and. y(ikr).le.y(isr)) go to 3800
      if(isnxt.eq.0) go to 3600
      if(isnxt.eq.islst) stop 1140
      if(vx(isnxt).ge.xsw) go to 3600
      iscur = isnxt
      go to 3500
 3600 continue
      nr = nr + 1
      if(nr .gt. nmax) stop 1150
      x(nr) = 2.0d0*xsw - x(ikr)
      y(nr) = y(ikr)
      ntot = ntot + 1
      if(ntot.gt.nv) stop 1160
c     WRITE(*,*)'IKR 2 =',IKR
      if(ikr.eq.inw) go to 3900
      isr  = ikr
      islst = iscur
      go to 3500
c
c     switch directions (to counterclockwise) around isr
c
 3800 continue
      if(icon(4,islst).eq.isr) then
         isnxt = icon(2,islst)
      elseif(icon(5,islst).eq.isr) then
         isnxt = icon(3,islst)
      elseif(icon(6,islst).eq.isr) then
         isnxt = icon(1,islst)
      else
         stop 1170
      endif
      islst = 0
      go to 3100
c
c     process NW corner and then NW-NE edge
c
 3900 continue
      nr = nr + 1
      if(nr .gt. nmax) stop 1200
      y(nr) = 2.0d0*ynw - y(inw)
      x(nr) = x(inw)
      ntot = 1
c     WRITE(*,*)' '
c     WRITE(*,*)'INW=',INW
      if(inw.eq.ine) go to 4900
      isr  = inw
      iscur= is(inw)
      if(iscur.le.0 .or. iscur.gt.ivnxt) stop 1210
      if(icon(4,iscur).eq.inw) then
         islst = icon(3,iscur)
         isnxt = icon(2,iscur)
         ikr = icon(5,iscur)
         icr = icon(6,iscur)
      elseif(icon(5,iscur).eq.inw) then
         islst = icon(1,iscur)
         isnxt = icon(3,iscur)
         ikr = icon(6,iscur)
         icr = icon(4,iscur)
      elseif(icon(6,iscur).eq.inw) then
         islst = icon(2,iscur)
         isnxt = icon(1,iscur)
         ikr = icon(4,iscur)
         icr = icon(5,iscur)
      else
         stop 1220
      endif
      if(islst.ne.0) islst = iscur
c
c     test for vertex in initial Voronoi cell below NW-NE line
c
      if(vy(iscur).le.ynw) go to 4300
c
c     test for bounded initial Voronoi cell
c
      if(islst.ne.0) go to 4100
c
c     test for unbounded edge moving into rectangle
c
      if(x(ikr).le.x(inw)) go to 4100
      nr = nr + 1
      if(nr .gt. nmax) stop 1230
      y(nr) = 2.0d0*ynw - y(ikr)
      x(nr) = x(ikr)
      ntot = 2
c     WRITE(*,*)'IKR 1 =',IKR
      if(ikr.eq.ine) go to 4900
      isr  = ikr
      islst = iscur
      go to 4500
c
c     find vertex of Voronoi cell below NW-NE line
c     (in counterclockwise direction around isr)
c
 4100 continue
      if(isnxt.eq.0 .or. isnxt.eq.islst) stop 1240
      iscur = isnxt
      if(vy(iscur).le.ynw) then
         if(islst.ne.0) islst = iscur
         go to 4200
      endif
      if(icon(4,iscur).eq.isr) then
         isnxt = icon(2,iscur)
      elseif(icon(5,iscur).eq.isr) then
         isnxt = icon(3,iscur)
      elseif(icon(6,iscur).eq.isr) then
         isnxt = icon(1,iscur)
      else
         stop 1250
      endif
      go to 4100
c
c     find vertex of Voronoi cell above NW-NE line
c     (in counterclockwise direction around isr)
c
 4200 continue
      if(icon(4,iscur).eq.isr) then
         isnxt = icon(2,iscur)
         icr = icon(6,iscur)
      elseif(icon(5,iscur).eq.isr) then
         isnxt = icon(3,iscur)
         icr = icon(4,iscur)
      elseif(icon(6,iscur).eq.isr) then
         isnxt = icon(1,iscur)
         icr = icon(5,iscur)
      else
         stop 1260
      endif
c
 4300 continue
      if(isnxt.eq.0) then
         if(islst.ne.0) stop 1270
         go to 4400
      endif
      if(isnxt.eq.islst) stop 1280
      iscur = isnxt
      if(vy(iscur).le.ynw) go to 4200
 4400 continue
      nr = nr + 1
      if(nr .gt. nmax) stop 1290
      y(nr) = 2.0d0*ynw - y(icr)
      x(nr) = x(icr)
      ntot = ntot + 1
      if(ntot.gt.nv) stop 1300
c     WRITE(*,*)'ICR=',ICR
      if(icr.eq.ine) go to 4900
      isr  = icr
      if(isnxt.eq.0) go to 4200
      islst = iscur
c
c     test vertices around isr in clockwise direction
c
 4500 continue
      if(icon(4,iscur).eq.isr) then
         isnxt = icon(3,iscur)
         ikr = icon(5,iscur)
      elseif(icon(5,iscur).eq.isr) then
         isnxt = icon(1,iscur)
         ikr = icon(6,iscur)
      elseif(icon(6,iscur).eq.isr) then
         isnxt = icon(2,iscur)
         ikr = icon(4,iscur)
      else
         stop 1320
      endif
c
      if(isnxt.eq.0 .and. x(ikr).le.x(isr)) go to 4800
      if(isnxt.eq.0) go to 4600
      if(isnxt.eq.islst) stop 1340
      if(vy(isnxt).le.ynw) go to 4600
      iscur = isnxt
      go to 4500
 4600 continue
      nr = nr + 1
      if(nr .gt. nmax) stop 1350
      y(nr) = 2.0d0*ynw - y(ikr)
      x(nr) = x(ikr)
      ntot = ntot + 1
      if(ntot.gt.nv) stop 1360
c     WRITE(*,*)'IKR 2 =',IKR
      if(ikr.eq.ine) go to 4900
      isr  = ikr
      islst = iscur
      go to 4500
c
c     switch directions (to counterclockwise) around isr
c
 4800 continue
      if(icon(4,islst).eq.isr) then
         isnxt = icon(2,islst)
      elseif(icon(5,islst).eq.isr) then
         isnxt = icon(3,islst)
      elseif(icon(6,islst).eq.isr) then
         isnxt = icon(1,islst)
      else
         stop 1370
      endif
      islst = 0
      go to 4100
c
 4900 continue
c
      return
      end
*PNTINS
c
c     subroutine pntins to -
c
c     insert a point into a Delaunay triangulation and obtain a
c     Delaunay triangulation that contains the inserted point
c
c     February 11, 1991
c
      subroutine pntins(x, y, ix, iy, ix2, iy2, is, icon, iave, iaze,
     *                  epf, mhalf, mfull, isclp, nvmz, ileft, iright,
     *                  n, iabe, iase, isite, nmax, nvmax, ivnxt,
     *                  insrt, itype, ipoint)
c
      integer mhalf, mfull, isclp(*)
      integer nvmz, ileft, iright, n, iabe, iase, isite
      integer nmax, nvmax, ivnxt, insrt, itype, ipoint
      double precision epf
      double precision x(*), y(*)
      integer ix(*), iy(*), ix2(*), iy2(*)
      integer is(*), icon(6,*), iave(*), iaze(*)
c
      integer ivfnd, ivcur, ivaft
c
c     find position of point in triangulation
c
      call gettri(x, y, ix, iy, ix2, iy2, is, icon, epf,
     *            mhalf, mfull, isclp, nvmz, ileft, iright,
     *            n, itype, ivfnd, ipoint)
c
c     insert point
c
      if(itype .eq. -1 .or. itype .eq. 1) go to 400
      if(insrt .eq. 0) go to 400
      if(isite.ne.0) go to 200
      isite = iaze(iase)
      if(iase .eq. 1) then
          n = isite
          if(n. gt. nmax) stop 1700
          iaze(1) = n + 1
      else
          iase = iase - 1
      endif
      x(isite) = x(ipoint)
      y(isite) = y(ipoint)
      ix(isite) = ix(ipoint)
      iy(isite) = iy(ipoint)
      ix2(isite) = ix2(ipoint)
      iy2(isite) = iy2(ipoint)
  200 continue
      if(itype .eq. 2) then
          call edgint(icon, is, iave, ivnxt, nvmax, nvmz, ivfnd, isite,
     *                iright, iabe)
      elseif(itype .eq. 3) then
          call intins(icon, is, iave, ivnxt, nvmax, nvmz, ivfnd, isite,
     *                iabe)
      else
          call outhul(x, y, ix, iy, ix2, iy2, icon, is, iave,
     *                ivnxt, nvmax, nvmz, mhalf, mfull, isclp,
     *                ivfnd, isite, iright, iabe, epf)
      endif
      iright = isite
c
c     optimize triangulation with respect to isite
c
      ivcur = ivfnd
  300 continue
      ivaft = icon(2,ivcur)
      call opttri(x, y, ix, iy, ix2, iy2, icon, is, epf,
     *            ivcur, isite, nvmz, mhalf, mfull, isclp)
      if(ivaft .eq. 0 .or. ivaft .eq. nvmz) go to 400
      ivcur = iabs(ivaft)
      if(ivcur .ne. ivfnd) go to 300
c
  400 continue
c
      return
      end
*GETTRI
c
c     subroutine gettri to -
c
c     obtain position in a triangulation for
c     a point (xdes,ydes) by moving through the triangulation from
c     a known vertex (ileft) in the triangulation to the point
c
c     the following holds on output: 
c
c     if itype equals -1 then (xdes,ydes) is exactly the vertex with
c     index ileft
c
c     if itype equals 0 then (xdes,ydes) is a point outside the convex
c     hull of the triangulation and iright is a vertex of the convex
c     hull visible from (xdes,ydes)
c
c     if itype equals 1 then (xdes,ydes) is exactly the vertex in the
c     triangulation with index iright and iright is different from
c     ileft
c
c     if itype equals 2 then (xdes,ydes) is a point in the interior of
c     a side of the triangle in the triangulation with index ivfnd and
c     iright is the first vertex of this triangle which is found when
c     moving from the point in a counterclockwise direction around the
c     boundary of the triangle
c
c     if itype equals 3 then (xdes,ydes) is a point in the interior of
c     the triangle in the triangulation with index ivfnd
c
c     February 6, 1991
c
      subroutine gettri(x, y, ix, iy, ix2, iy2, is, icon, epf,
     *                  mhalf, mfull, isclp, nvmz, ileft,
     *                  iright, n, itype, ivfnd, ipoint)
c
      integer mhalf, mfull, isclp(*), ipoint
      integer nvmz, ileft, iright, n, itype, ivfnd
      double precision epf, xdes, ydes
      double precision x(*), y(*)
      integer ix(*), iy(*), ix2(*), iy2(*)
      integer is(*), icon(6,*)
c
      integer nkmax
      parameter (nkmax=30)
      integer iax(nkmax), iay(nkmax), isgax, isgay, ikax, ikay
      integer iox(nkmax), ioy(nkmax), isgox, isgoy, ikox, ikoy
      integer isi1, isi2, isi3, iloft, iperp, nsgax, isga, i
      integer ivini, ivcur, ivadj, ivfol, ivpre, ivlst
      double precision delx, dely, xtemp, ytemp, dls1, dls2, dls3
      double precision xtamp, ytamp, dlg1, dlg2, dlg3, dlen, daf
      double precision dot1, dot2, dot3, dut2
      double precision xtump, ytump, dalx, daly
c
c     compute direction parameters
c
      if(ileft .lt. 1 .or. ileft .gt. n) stop 1800
      iperp = 0
      xdes = x(ipoint)
      ydes = y(ipoint)
      delx = y(ileft) - ydes
      dely = xdes - x(ileft)
      dlen = dsqrt(delx**2 + dely**2)
      if(dlen .lt. epf) iperp = 1
      call prpvec(ix, iy, ix2, iy2, ileft, ipoint, mhalf, mfull, isclp,
     *            iax, isgax, ikax, iay, isgay, ikay)
      if(isgax.eq.0 .and. isgay.eq.0) then
         itype = -1
         iright = ileft
         go to 2000
      endif
      nsgax = -isgax
      isi2 = ileft
      itype = 0
c
  100 continue
      iloft = isi2
      ivini = is(iloft)
      ivcur = ivini
      if(icon(4,ivcur) .eq. iloft) then
         ivlst = icon(3,ivcur)
      elseif(icon(5,ivcur) .eq. iloft) then
         ivlst = icon(1,ivcur)
      elseif(icon(6,ivcur) .eq. iloft) then
         ivlst = icon(2,ivcur)
      else
          stop 1805
      endif
c
c     test current triangle
c
  200 continue
      if(icon(4,ivcur) .eq. iloft) then
         ivadj = icon(2,ivcur)
         ivfol = icon(1,ivcur)
         isi2 = icon(5,ivcur)
         isi1 = icon(6,ivcur)
      elseif(icon(5,ivcur) .eq. iloft) then
         ivadj = icon(3,ivcur)
         ivfol = icon(2,ivcur)
         isi2 = icon(6,ivcur)
         isi1 = icon(4,ivcur)
      elseif(icon(6,ivcur) .eq. iloft) then
         ivadj = icon(1,ivcur)
         ivfol = icon(3,ivcur)
         isi2 = icon(4,ivcur)
         isi1 = icon(5,ivcur)
      else
         stop 1810
      endif
c
      xtemp = x(isi2) - xdes
      ytemp = y(isi2) - ydes
      dls2 = dsqrt(xtemp**2 + ytemp**2)
      if(dls2 .ge. epf) go to 220
      call vectrd(ix, iy, ix2, iy2, ipoint, isi2, mhalf, mfull,
     *            isclp, iox, isgox, ikox, ioy, isgoy, ikoy)
      if(isgox.ne.0 .or. isgoy.ne.0) go to 220
      itype = 1
      iright = isi2
      go to 2000
  220 continue
      if(iperp.eq.1) go to 240
      xtamp = x(isi2) - x(iloft)
      ytamp = y(isi2) - y(iloft)
      dlg2 = dsqrt(xtamp**2 + ytamp**2)
      if(dlg2 .lt. epf) go to 240
      dot2 = (delx * xtamp + dely * ytamp)/dlen
      if(dot2.ge.epf .or. dot2.le.-epf) go to 260
  240 continue
c     WRITE(*,*)'GETTRI 1 DISTSN'
      call distsn(ix, iy, ix2, iy2, iloft, isi2, mhalf, mfull, isclp,
     *            iax, isgax, ikax, iay, isgay, ikay, isga)
      if(isga.gt.0) then
         dot2 = 10.0d0
      elseif(isga.lt.0) then
         dot2 =-10.0d0
      else
         dot2 = 0.0d0
      endif
  260 continue
      if(dot2 .ge. epf) go to 500
      if(dot2 .le.-epf) go to 320
      if(iperp.eq.1) go to 280
      if(dlg2 .lt. epf) go to 280
      dut2 = (dely * xtamp - delx * ytamp)/dlen
      if(dut2.ge.epf .or. dut2.le.-epf) go to 300
  280 continue
      call distsn(ix, iy, ix2, iy2, iloft, isi2, mhalf, mfull, isclp,
     *            iay, isgay, ikay, iax, nsgax, ikax, isga)
      if(isga.gt.0) then
         dut2 = 10.0d0
      elseif(isga.lt.0) then
         dut2 =-10.0d0
      else
         stop 1815
      endif
  300 continue
      if(dut2 .ge. epf) go to 400
      dot1 = -10.0
      go to 500
  320 continue
      xtump = x(isi1) - xdes
      ytump = y(isi1) - ydes
      dls1 = dsqrt(xtump**2 + ytump**2)
      if(dls1 .ge. epf) go to 340
      call vectrd(ix, iy, ix2, iy2, ipoint, isi1, mhalf, mfull,
     *            isclp, iox, isgox, ikox, ioy, isgoy, ikoy)
      if(isgox.ne.0 .or. isgoy.ne.0) go to 340
      itype = 1
      iright = isi1
      go to 2000
  340 continue
      if(iperp.eq.1) go to 360
      xtamp = x(isi1) - x(iloft)
      ytamp = y(isi1) - y(iloft)
      dlg1 = dsqrt(xtamp**2 + ytamp**2)
      if(dlg1 .lt. epf) go to 360
      dot1 = (delx * xtamp + dely * ytamp)/dlen
      if(dot1.ge.epf .or. dot1.le.-epf) go to 380
  360 continue
c     WRITE(*,*)'GETTRI 2 DISTSN'
      call distsn(ix, iy, ix2, iy2, iloft, isi1, mhalf, mfull, isclp,
     *            iax, isgax, ikax, iay, isgay, ikay, isga)
      if(isga.gt.0) then
         dot1 = 10.0d0
      elseif(isga.lt.0) then
         dot1 =-10.0d0
      else
         dot1 = 0.0d0
      endif
  380 continue
      if(dot1 .lt. epf) go to 500
      if(dot2 .gt. -epf) stop 1820
      go to 900
  400 continue
      if(iperp.eq.1) go to 420
      if(dls2 .lt. epf) go to 420
      dot3 = (-dely * xtemp + delx * ytemp)/dlen
      if(dot3.ge.epf .or. dot3.le.-epf) go to 440
  420 continue
      call distsn(ix, iy, ix2, iy2, isi2, ipoint, mhalf, mfull, isclp,
     *            iay, isgay, ikay, iax, nsgax, ikax, isga)
      if(isga.gt.0) then
         dot3 = 10.0d0
      elseif(isga.lt.0) then
         dot3 =-10.0d0
      else
         stop 1825
      endif
  440 continue
      if(dot3 .ge. epf) go to 100
      itype = 2
      iright = isi2
      ivfnd = ivcur
      go to 2000
c
  500 continue
      ivpre = ivcur
      ivcur = iabs(ivadj)
      if(ivcur .eq. ivini) stop 1830
      if(ivadj .ne. 0 .and. ivadj .ne. nvmz) go to 200
      if(ivlst .ne. 0 .and. ivlst .ne. nvmz) stop 1835
      iright = iloft
      if(dot2 .ge. epf) go to 2000
      if(dot1 .le. -epf) go to 2000
      if(dot2 .gt. -epf .or. dot1 .ge. epf) stop 1840
      if(iperp.eq.1) go to 620
      if(dls1 .lt. epf) go to 620
      dot3 = (-dely * xtump + delx * ytump)/dlen
      if(dot3.ge.epf .or. dot3.le.-epf) go to 640
  620 continue
      call distsn(ix, iy, ix2, iy2, isi1, ipoint, mhalf, mfull, isclp,
     *            iay, isgay, ikay, iax, nsgax, ikax, isga)
      if(isga.gt.0) then
         dot3 = 10.0d0
      elseif(isga.lt.0) then
         dot3 =-10.0d0
      else
         stop 1845
      endif
  640 continue
      isi2 = isi1
      if(dot3 .ge. epf) go to 100
      itype = 2
      ivfnd = ivpre
      go to 2000
c
  900 continue
      if(dls2.lt.epf) go to 920
      dalx = y(isi2) - y(isi1)
      daly = x(isi1) - x(isi2)
      daf = dsqrt(dalx**2 + daly**2)
      if(daf .lt. epf) go to 920
      dot3 = (dalx * xtemp + daly * ytemp)/daf
      if(dot3.ge.epf .or. dot3.le.-epf) go to 940
  920 continue
      call innsgn(ix, iy, ix2, iy2, isi1, isi2, ipoint, mhalf, mfull,
     *            isclp, isga)
      if(isga.gt.0) then
         dot3 = 10.0d0
      elseif(isga.lt.0) then
         dot3 =-10.0d0
      else
         dot3 = 0.0d0
      endif
  940 continue
      if(dot3 .ge. epf) go to 1100
      if(dot3 .le. -epf) go to 1000
      itype = 2
      iright = isi1
      ivfnd = ivcur
      go to 2000
 1000 continue
      itype = 3
      ivfnd = ivcur
      go to 2000
c
c     test next triangle
c
 1100 continue
      iright = isi2
      if(ivfol .eq. 0 .or. ivfol .eq. nvmz) go to 2000
      ivcur = iabs(ivfol)
      if(icon(4,ivcur) .eq. isi2) then
          isi3 = icon(5,ivcur)
          ivadj = icon(1,ivcur)
          ivpre = icon(3,ivcur)
      elseif(icon(5,ivcur) .eq. isi2) then
          isi3 = icon(6,ivcur)
          ivadj = icon(2,ivcur)
          ivpre = icon(1,ivcur)
      elseif(icon(6,ivcur) .eq. isi2) then
          isi3 = icon(4,ivcur)
          ivadj = icon(3,ivcur)
          ivpre = icon(2,ivcur)
      else
          stop 1850
      endif
      xtemp = x(isi3) - xdes
      ytemp = y(isi3) - ydes
      dls3 = dsqrt(xtemp**2 + ytemp**2)
      if(dls3 .ge. epf) go to 1200
      call vectrd(ix, iy, ix2, iy2, ipoint, isi3, mhalf, mfull,
     *            isclp, iox, isgox, ikox, ioy, isgoy, ikoy)
      if(isgox.ne.0 .or. isgoy.ne.0) go to 1200
      itype = 1
      iright = isi3
      go to 2000
 1200 continue
      if(iperp.eq.1) go to 1220
      xtamp = x(isi3) - x(iloft)
      ytamp = y(isi3) - y(iloft)
      dlg3 = dsqrt(xtamp**2 + ytamp**2)
      if(dlg3 .lt. epf) go to 1220
      dot3 = (delx * xtamp + dely * ytamp)/dlen
      if(dot3.ge.epf .or. dot3.le.-epf) go to 1240
 1220 continue
      call distsn(ix, iy, ix2, iy2, iloft, isi3, mhalf, mfull, isclp,
     *            iax, isgax, ikax, iay, isgay, ikay, isga)
      if(isga.gt.0) then
         dot3 = 10.0d0
      elseif(isga.lt.0) then
         dot3 =-10.0d0
      else
         dot3 = 0.0d0
      endif
 1240 continue
      if(dot3 .ge. epf) go to 1400
      if(dot3 .le. -epf) go to 1700
      if(iperp.eq.1) go to 1260
      if(dls3 .lt. epf) go to 1260
      dot3 = (-dely * xtemp + delx * ytemp)/dlen
      if(dot3.ge.epf .or. dot3.le.-epf) go to 1280
 1260 continue
      call distsn(ix, iy, ix2, iy2, isi3, ipoint, mhalf, mfull, isclp,
     *            iay, isgay, ikay, iax, nsgax, ikax, isga)
      if(isga.gt.0) then
         dot3 = 10.0d0
      elseif(isga.lt.0) then
         dot3 =-10.0d0
      else
         stop 1855
      endif
 1280 continue
      isi2 = isi3
      if(dot3 .ge. epf) go to 100
      itype = 3
      ivfnd = ivcur
      go to 2000
c
 1400 continue
      if(dls3.lt.epf) go to 1420
      dalx = y(isi3) - y(isi2)
      daly = x(isi2) - x(isi3)
      daf = dsqrt(dalx**2 + daly**2)
      if(daf .lt. epf) go to 1420
      dot3 = (dalx * xtemp + daly * ytemp)/daf
      if(dot3.ge.epf .or. dot3.le.-epf) go to 1440
 1420 continue
      call innsgn(ix, iy, ix2, iy2, isi2, isi3, ipoint, mhalf, mfull,
     *            isclp, isga)
      if(isga.gt.0) then
         dot3 = 10.0d0
      elseif(isga.lt.0) then
         dot3 =-10.0d0
      else
         dot3 = 0.0d0
      endif
 1440 continue
      if(dot3 .le. -epf) go to 1600
      if(dot3 .ge. epf) go to 1500
      itype = 2
      iright = isi3
      ivfnd = ivcur
      go to 2000
 1500 continue
      itype = 3
      ivfnd = ivcur
      go to 2000
 1600 continue
      isi1 = isi3
      ivfol = ivpre
      go to 1100
c
 1700 continue
      if(dls3.lt.epf) go to 1720
      dalx = y(isi3) - y(isi1)
      daly = x(isi1) - x(isi3)
      daf = dsqrt(dalx**2 + daly**2)
      if(daf .lt. epf) go to 1720
      dot3 = (dalx * xtemp + daly * ytemp)/daf
      if(dot3.ge.epf .or. dot3.le.-epf) go to 1740
 1720 continue
      call innsgn(ix, iy, ix2, iy2, isi1, isi3, ipoint, mhalf, mfull,
     *            isclp, isga)
      if(isga.gt.0) then
         dot3 = 10.0d0
      elseif(isga.lt.0) then
         dot3 =-10.0d0
      else
         dot3 = 0.0d0
      endif
 1740 continue
      if(dot3 .ge. epf) go to 1900
      if(dot3 .le. -epf) go to 1800
      itype = 2
      iright = isi1
      ivfnd = ivcur
      go to 2000
 1800 continue
      itype = 3
      ivfnd = ivcur
      go to 2000
 1900 continue
      isi2 = isi3
      ivfol = ivadj
      go to 1100
c
c     test position of point
c
 2000 continue
c     WRITE(*,*)'IPOINT=',IPOINT,' ITYPE=',ITYPE
      if(itype .eq. -1 .or. itype. eq. 1) then
         if(itype .eq. -1 .and. ileft .ne. iright) stop 1860
         if(itype .eq.  1 .and. ileft .eq. iright) stop 1865
         call vectrd(ix, iy, ix2, iy2, ipoint, iright, mhalf, mfull,
     *               isclp, iox, isgox, ikox, ioy, isgoy, ikoy)
         if(isgox.ne.0 .or. isgoy.ne.0) stop 1870
      elseif(itype .eq. 0) then
         ivcur = is(iright)
         if(icon(4,ivcur) .eq. iright) then
            ivlst = icon(3,ivcur)
            isi2 = icon(5,ivcur)
         elseif(icon(5,ivcur) .eq. iright) then
            ivlst = icon(1,ivcur)
            isi2 = icon(6,ivcur)
         elseif(icon(6,ivcur) .eq. iright) then
            ivlst = icon(2,ivcur)
            isi2 = icon(4,ivcur)
         else
            stop 1875
         endif
         if(ivlst .ne. 0 .and. ivlst .ne. nvmz) stop 1880
         iperp = 0
         delx = y(ipoint) - y(iright)
         dely = x(iright) - x(ipoint)
         dlen = dsqrt(delx**2+dely**2)
         if(dlen .lt. epf) go to 2020
         xtemp = x(isi2) - x(iright)
         ytemp = y(isi2) - y(iright)
         dls2 = dsqrt(xtemp**2+ytemp**2)
         if(dls2 .lt. epf) go to 2020
         dot2 = (delx * xtemp + dely * ytemp)/dlen
         if(dot2.ge.epf .or. dot2.le.-epf) go to 2040
 2020    continue
         iperp = 1
         call prpvec(ix, iy, ix2, iy2, ipoint, iright, mhalf, mfull,
     *               isclp, iax, isgax, ikax, iay, isgay, ikay)
         if(isgax.eq.0 .and. isgay.eq.0) stop 1885
         call distsn(ix, iy, ix2, iy2, iright, isi2, mhalf, mfull,
     *               isclp, iax, isgax, ikax, iay, isgay, ikay, isga)
         if(isga.gt.0) then
            dot2 = 10.0d0
         elseif(isga.lt.0) then
            dot2 =-10.0d0
         else
            dot2 = 0.0d0
         endif
 2040    continue
         if(dot2.le.-epf) go to 3000
c
         ivini = ivcur
 2060    continue
         if(icon(4,ivcur) .eq. iright) then
            ivadj = icon(2,ivcur)
            isi1 = icon(6,ivcur)
         elseif(icon(5,ivcur) .eq. iright) then
            ivadj = icon(3,ivcur)
            isi1 = icon(4,ivcur)
         elseif(icon(6,ivcur) .eq. iright) then
            ivadj = icon(1,ivcur)
            isi1 = icon(5,ivcur)
         else
            stop 1890
         endif
         ivcur = iabs(ivadj)
         if(ivcur.eq.ivini) stop 1895
         if(ivadj .ne. 0 .and. ivadj .ne. nvmz) go to 2060
         if(dlen .lt. epf) go to 2080
         xtemp = x(isi1) - x(iright)
         ytemp = y(isi1) - y(iright)
         dls2 = dsqrt(xtemp**2+ytemp**2)
         if(dls2 .lt. epf) go to 2080
         dot2 = (delx * xtemp + dely * ytemp)/dlen
         if(dot2.ge.epf .or. dot2.le.-epf) go to 2100
 2080    continue
         if(iperp.eq.1) go to 2090
         call prpvec(ix, iy, ix2, iy2, ipoint, iright, mhalf, mfull,
     *               isclp, iax, isgax, ikax, iay, isgay, ikay)
         if(isgax.eq.0 .and. isgay.eq.0) stop 1898
 2090    continue
         call distsn(ix, iy, ix2, iy2, iright, isi1, mhalf, mfull,
     *               isclp, iax, isgax, ikax, iay, isgay, ikay, isga)
         if(isga.gt.0) then
            dot2 = 10.0d0
         elseif(isga.lt.0) then
            dot2 =-10.0d0
         else
            dot2 = 0.0d0
         endif
 2100    continue
         if(dot2.lt.epf) stop 1900
      elseif(itype .eq. 2 .or. itype .eq. 3) then
         isi2 = icon(6,ivfnd)
         do 2300 i=4,6
            isi1 = icon(i,ivfnd)
            xtemp = x(isi2) - xdes
            ytemp = y(isi2) - ydes
            dls2 = dsqrt(xtemp**2 + ytemp**2)
            if(dls2.lt.epf) go to 2220
            dalx = y(isi2) - y(isi1)
            daly = x(isi1) - x(isi2)
            daf = dsqrt(dalx**2 + daly**2)
            if(daf .lt. epf) go to 2220
            dot3 = (dalx * xtemp + daly * ytemp)/daf
            if(dot3.ge.epf .or. dot3.le.-epf) go to 2240
 2220       continue
            call innsgn(ix, iy, ix2, iy2, isi1, isi2, ipoint, mhalf,
     *                  mfull, isclp, isga)
            if(isga.gt.0) then
               dot3 = 10.0d0
            elseif(isga.lt.0) then
               dot3 =-10.0d0
            else
               dot3 = 0.0d0
            endif
 2240       continue
            if(itype .eq. 2) then
               if(iright.eq.isi1) then
                  if(dot3.ge.epf .or. dot3.le.-epf) stop 1905
               else
                  if(dot3.gt.-epf) stop 1910
               endif
            else
               if(dot3.gt.-epf) stop 1915
            endif
            isi2 = isi1
 2300    continue
      else
         stop 1920
      endif
c
 3000 continue
      return
      end
*EDGINT
c
c     subroutine edgint to -
c
c     insert a point into current triangulation when point is in the
c     interior of an edge of a triangle
c
c     January 31, 1991
c
      subroutine edgint(icon, is, iave, ivnxt, nvmax, nvmz, ivcur,
     *                  isite, islt, iabe)
c
      integer ivnxt, nvmax, nvmz, ivcur, isite, islt, iabe
      integer icon(6,*), is(*), iave(*)
c
      integer it1, it2, it3, it4, iv1, iv2, iv3, iv4
      integer isrt, isad, ivadj, isop, ivodj
c
      it1 = ivcur
      if(icon(4,ivcur) .eq. islt) then
          isrt = icon(6,ivcur)
          isad = icon(5,ivcur)
          ivadj = icon(2,ivcur)
          iv1 = icon(3,ivcur)
          iv2 = icon(1,ivcur)
      elseif(icon(5,ivcur) .eq. islt) then
          isrt = icon(4,ivcur)
          isad = icon(6,ivcur)
          ivadj = icon(3,ivcur)
          iv1 = icon(1,ivcur)
          iv2 = icon(2,ivcur)
      elseif(icon(6,ivcur) .eq. islt) then
          isrt = icon(5,ivcur)
          isad = icon(4,ivcur)
          ivadj = icon(1,ivcur)
          iv1 = icon(2,ivcur)
          iv2 = icon(3,ivcur)
      else
          stop 1930
      endif
c
      it3 = iabs(ivadj)
      it4 = it3
      if(ivadj .eq. 0 .or. ivadj .eq. nvmz) go to 100
      if(icon(4,it4) .eq. islt) then
          isop = icon(6,it4)
          iv3 = icon(1,it4)
          iv4 = icon(2,it4)
      elseif(icon(5,it4) .eq. islt) then
          isop = icon(4,it4)
          iv3 = icon(2,it4)
          iv4 = icon(3,it4)
      elseif(icon(6,it4) .eq. islt) then
          isop = icon(5,it4)
          iv3 = icon(3,it4)
          iv4 = icon(1,it4)
      else
          stop 1935
      endif
c
  100 continue
      it2 = iave(iabe)
      if(iabe .eq. 1) then
          ivnxt = it2
          if(ivnxt .gt. nvmax) stop 1940
          iave(1) = ivnxt + 1
      else
          iabe = iabe - 1
      endif
      if(ivadj .eq. 0 .or. ivadj .eq. nvmz) go to 200
      it3 = iave(iabe)
      if(iabe .eq. 1) then
          ivnxt = it3
          if(ivnxt .gt. nvmax) stop 1945
          iave(1) = ivnxt + 1
      else
          iabe = iabe - 1
      endif
c
  200 continue
      is(isite) = it1
      if(is(isad) .eq. it1) is(isad) = it2
      if(is(isrt) .eq. it1) is(isrt) = it2
      if(is(isrt) .eq. it4) is(isrt) = it3
c
c     create new triangles
c
      icon(1,it1) = iv1
      icon(2,it1) = it2
      icon(3,it1) = it4
      icon(4,it1) = isite
      icon(5,it1) = islt
      icon(6,it1) = isad
c
      icon(1,it2) = iv2
      icon(2,it2) = it3
      icon(3,it2) = it1
      icon(4,it2) = isite
      icon(5,it2) = isad
      icon(6,it2) = isrt
c
      if(islt.eq.isad.or.isrt.eq.isad.or.isite.eq.islt.or.
     *   isite.eq.isrt.or.isite.eq.isad) stop 1950
c
      if(ivadj .lt. 0) then
          icon(3,it1) = -it4
          icon(2,it2) = -it3
      endif
c
      if(ivadj .eq. 0 .or. ivadj .eq. nvmz) go to 300
      icon(1,it3) = iv3
      icon(2,it3) = it4
      icon(3,it3) = it2
      icon(4,it3) = isite
      icon(5,it3) = isrt
      icon(6,it3) = isop
c
      icon(1,it4) = iv4
      icon(2,it4) = it1
      icon(3,it4) = it3
      icon(4,it4) = isite
      icon(5,it4) = isop
      icon(6,it4) = islt
c
      if(islt.eq.isop.or.isrt.eq.isop.or.isite.eq.isop) stop 1955
c
      if(ivadj .lt. 0) then
          icon(3,it3) = -it2
          icon(2,it4) = -it1
      endif
c
c     update neighboring triangles
c
  300 continue
      if(iv2 .eq. 0 .or. iv2 .eq. nvmz) go to 400
      ivodj = it1
      if(iv2 .lt. 0) then
          iv2 = -iv2
          it2 = -it2
          ivodj = -ivodj
      endif
c
      if(icon(1,iv2) .eq. ivodj) then
          icon(1,iv2) = it2
      elseif(icon(2,iv2) .eq. ivodj) then
          icon(2,iv2) = it2
      elseif(icon(3,iv2) .eq. ivodj) then
          icon(3,iv2) = it2
      else
          stop 1960
      endif
  400 continue
c
      if(ivadj .eq. 0 .or. ivadj .eq. nvmz) go to 500
      if(iv3 .eq. 0 .or. iv3 .eq. nvmz) go to 500
      ivodj = it4
      if(iv3 .lt. 0) then
          iv3 = -iv3
          it3 = -it3
          ivodj = -ivodj
      endif
c
      if(icon(1,iv3) .eq. ivodj) then
          icon(1,iv3) = it3
      elseif(icon(2,iv3) .eq. ivodj) then
          icon(2,iv3) = it3
      elseif(icon(3,iv3) .eq. ivodj) then
          icon(3,iv3) = it3
      else
          stop 1970
      endif
  500 continue
c
      return
      end
*INTINS
c
c     subroutine intins to -
c
c     insert a point into current triangulation when point is in
c     the interior of a triangle
c
c     January 31, 1991
c
      subroutine intins(icon, is, iave, ivnxt, nvmax, nvmz, ivcur,
     *                  isite, iabe)
c
      integer ivnxt, nvmax, nvmz, ivcur, isite, iabe
      integer icon(6,*), is(*), iave(*)
c
      integer it1, it2, it3, iv1, iv2, iv3, is4, is5, is6, ivadj
c
      it1 = ivcur
      it2 = iave(iabe)
      if(iabe .eq. 1) then
          ivnxt = it2
          if(ivnxt .gt. nvmax) stop 2000
          iave(1) = ivnxt + 1
      else
          iabe = iabe - 1
      endif
      it3 = iave(iabe)
      if(iabe .eq. 1) then
          ivnxt = it3
          if(ivnxt .gt. nvmax) stop 2010
          iave(1) = ivnxt + 1
      else
          iabe = iabe - 1
      endif
c
      iv1 = icon(1,ivcur)
      iv2 = icon(2,ivcur)
      iv3 = icon(3,ivcur)
c
      is4 = icon(4,ivcur)
      is5 = icon(5,ivcur)
      is6 = icon(6,ivcur)
c
      if(is4.eq.is5.or.is4.eq.is6.or.is5.eq.is6.or.isite.eq.is4.or.
     *   isite.eq.is5.or.isite.eq.is6) stop 2020
c
      is(isite) = it1
      if(is(is4) .eq. ivcur) is(is4) = it3
      if(is(is6) .eq. ivcur) is(is6) = it2
c
c     create new triangles
c
      icon(1,it1) = iv1
      icon(2,it1) = it2
      icon(3,it1) = it3
      icon(4,it1) = isite
      icon(5,it1) = is5
      icon(6,it1) = is6
c
      icon(1,it2) = iv2
      icon(2,it2) = it3
      icon(3,it2) = it1
      icon(4,it2) = isite
      icon(5,it2) = is6
      icon(6,it2) = is4
c
      icon(1,it3) = iv3
      icon(2,it3) = it1
      icon(3,it3) = it2
      icon(4,it3) = isite
      icon(5,it3) = is4
      icon(6,it3) = is5
c
c     update neighboring triangles
c
      if(iv2 .eq. 0 .or. iv2 .eq. nvmz) go to 100
      ivadj = ivcur
      if(iv2 .lt. 0) then
          iv2 = -iv2
          it2 = -it2
          ivadj = -ivadj
      endif
c
      if(icon(1,iv2) .eq. ivadj) then
          icon(1,iv2) = it2
      elseif(icon(2,iv2) .eq. ivadj) then
          icon(2,iv2) = it2
      elseif(icon(3,iv2) .eq. ivadj) then
          icon(3,iv2) = it2
      else
          stop 2030
      endif
  100 continue
c
      if(iv3 .eq. 0 .or. iv3 .eq. nvmz) go to 200
      ivadj = ivcur
      if(iv3 .lt. 0) then
          iv3 = -iv3
          it3 = -it3
          ivadj = -ivadj
      endif
c
      if(icon(1,iv3) .eq. ivadj) then
          icon(1,iv3) = it3
      elseif(icon(2,iv3) .eq. ivadj) then
          icon(2,iv3) = it3
      elseif(icon(3,iv3) .eq. ivadj) then
          icon(3,iv3) = it3
      else
          stop 2040
      endif
  200 continue
c
      return
      end
*OUTHUL
c
c     subroutine outhul to -
c
c     insert a point that lies outside current triangulation
c
c     February 20, 1991
c
      subroutine outhul(x, y, ix, iy, ix2, iy2, icon, is, iave,
     *                  ivnxt, nvmax, nvmz, mhalf, mfull, isclp,
     *                  ivfnd, isite, islt, iabe, epf)
c
      integer ivnxt, nvmax, nvmz, mhalf, mfull, isclp(*)
      integer ivfnd, isite, islt, iabe
      double precision epf
      double precision x(*), y(*)
      integer ix(*), iy(*), ix2(*), iy2(*)
      integer icon(6,*), is(*), iave(*)
c
      integer ivcur, isrt, ivadj, isfn, ivlst, indx, ivcar, ivfad
      integer ivfin, ivnow, ivout, ivini, isga
      double precision delx, dely, dif, xtemp, ytemp, dlst, dot
c
c     search for first visible site from isite
c
      ivcur = is(islt)
      ivini = ivcur
  100 continue
      if(icon(4,ivcur) .eq. islt) then
          isrt = icon(6,ivcur)
          ivadj = icon(2,ivcur)
      elseif(icon(5,ivcur) .eq. islt) then
          isrt = icon(4,ivcur)
          ivadj = icon(3,ivcur)
      elseif(icon(6,ivcur) .eq. islt) then
          isrt = icon(5,ivcur)
          ivadj = icon(1,ivcur)
      else
          stop 2100
      endif
c
      if(ivadj .eq. 0 .or. ivadj .eq. nvmz) go to 200
      ivcur = iabs(ivadj)
      if(ivcur.eq.ivini) stop 2105
      go to 100
c
  200 continue
      delx = y(isite) - y(isrt)
      dely = x(isrt) - x(isite)
      dif = dsqrt(delx ** 2 + dely ** 2)
      if(dif .lt. epf) go to 220
      xtemp = x(islt) - x(isite)
      ytemp = y(islt) - y(isite)
      dlst=dsqrt(xtemp**2+ytemp**2)
      if(dlst .lt. epf) go to 220
      dot = (delx * xtemp + dely * ytemp)/dif
      if(dot.ge.epf .or. dot.le.-epf) go to 240
  220 continue
c     WRITE(*,*)'OUTHUL INNSGN 1'
      call innsgn(ix, iy, ix2, iy2, isite, isrt, islt, mhalf, mfull,
     *            isclp, isga)
      if(isga.gt.0) then
         dot = 10.0d0
      elseif(isga.lt.0) then
         dot =-10.0d0
      else
         dot = 0.0d0
      endif
  240 continue
      if(dot .gt. -epf) go to 300
      islt = isrt
      ivini = ivcur
      go to 100
c
c     generate triangles with isite as a vertex
c
  300 continue
      isfn = islt
      ivlst = 0
  400 continue
      ivcur = is(islt)
      isrt = islt
      if(icon(4,ivcur) .eq. isrt) then
          indx = 3
          islt = icon(5,ivcur)
          ivadj = icon(3,ivcur)
      elseif(icon(5,ivcur) .eq. isrt) then
          indx = 1
          islt = icon(6,ivcur)
          ivadj = icon(1,ivcur)
      elseif(icon(6,ivcur) .eq. isrt) then
          indx = 2
          islt = icon(4,ivcur)
          ivadj = icon(2,ivcur)
      else
          stop 2120
      endif
      if(ivadj .ne. 0 .and. ivadj .ne. nvmz) stop 2130
c
      delx = y(isite) - y(isrt)
      dely = x(isrt) - x(isite)
      dif = dsqrt(delx ** 2 + dely ** 2)
      if(dif .lt. epf) go to 420
      xtemp = x(islt) - x(isite)
      ytemp = y(islt) - y(isite)
      dlst=dsqrt(xtemp**2+ytemp**2)
      if(dlst .lt. epf) go to 420
      dot = (delx * xtemp + dely * ytemp)/dif
c     WRITE(*,*)'2 ISITE=',ISITE,' ISRT=',ISRT,' ISLT=',ISLT,' DOT=',DOT
      if(dot.ge.epf .or. dot.le.-epf) go to 440
  420 continue
c     WRITE(*,*)'OUTHUL INNSGN 2'
      call innsgn(ix, iy, ix2, iy2, isite, isrt, islt, mhalf, mfull,
     *            isclp, isga)
      if(isga.gt.0) then
         dot = 10.0d0
      elseif(isga.lt.0) then
         dot =-10.0d0
      else
         dot = 0.0d0
      endif
  440 continue
      if(dot .gt. -epf) go to 500
c
      ivfnd = iave(iabe)
      if(iabe .eq. 1) then
          ivnxt = ivfnd
          if(ivnxt .gt. nvmax) stop 2150
          iave(1) = ivnxt + 1
      else
          iabe = iabe - 1
      endif
c
      ivcar = ivcur
      ivfad = ivfnd
      if(ivadj .eq. nvmz) then
          ivcar = -ivcar
          ivfad = -ivfad
      endif
      icon(indx,ivcur) = ivfad
      icon(1,ivfnd) = ivcar
      icon(2,ivfnd) = ivlst
      icon(4,ivfnd) = isite
      icon(5,ivfnd) = islt
      icon(6,ivfnd) = isrt
c
      if(islt.eq.isrt.or.isite.eq.islt.or.isite.eq.isrt) stop 2160
c
      if(ivlst .ne. 0) icon(3,ivlst) = ivfnd
      ivlst = ivfnd
      if(isrt .ne. isfn) then
          if(ivadj .eq. nvmz) go to 400
          ivfin = icon(2,ivfnd)
          ivnow = ivcur
  450     continue
          if(icon(4,ivnow) .eq. isrt) then
              ivout = icon(2,ivnow)
          elseif(icon(5,ivnow) .eq. isrt) then
              ivout = icon(3,ivnow)
          elseif(icon(6,ivnow) .eq. isrt) then
              ivout = icon(1,ivnow)
          else
              stop 2170
          endif
          if(ivout .lt. 0) then
              is(isrt) = -ivout
              go to 400
          endif
          if(ivout .eq. ivfin) go to 400
          ivnow = ivout
          go to 450
      else
          is(isrt) = ivfnd
          go to 400
      endif
c
  500 continue
      if(ivlst .eq. 0) stop 2180
      icon(3,ivlst) = 0
      is(isite) = ivfnd
c
      return
      end
*OPTTRI
c
c     subroutine opttri to -
c
c     optimize triangulation by edge swapping with respect to a point,
c     in the direction of a triangle that has the point as a vertex
c
c     February 15, 1991
c
      subroutine opttri(x, y, ix, iy, ix2, iy2, icon, is, epf,
     *                 ivcur, isite, nvmz, mhalf, mfull, isclp)
c
      integer ivcur, isite, nvmz, mhalf, mfull, isclp(*)
      double precision x(*), y(*)
      integer ix(*), iy(*), ix2(*), iy2(*)
      integer icon(6,*), is(*)
      double precision zlamb, dot2, dot, dat, dif1, dif2, dif3
      double precision delx1, dely1, delx2, dely2, delx3, dely3
      double precision epf, vx, vy, tdist
      integer ishat, isadj, iscur, isfin, ivadj
      integer ivout, ivcar, ivodj, itide
c
c     initialize
c
      isadj = icon(5,ivcur)
      iscur = icon(6,ivcur)
      isfin = isadj
c
      if(isite.eq.isadj.or.isite.eq.iscur.or.isadj.eq.iscur) stop 2200
c
  100 continue
      if(icon(1,ivcur) .le. 0) go to 500
      ivadj = icon(1,ivcur)
c
      if(icon(2,ivadj) .eq. ivcur) go to 150
      if(icon(1,ivadj) .eq. ivcur) then
          icon(1,ivadj) = icon(3,ivadj)
          icon(3,ivadj) = icon(2,ivadj)
          icon(5,ivadj) = icon(4,ivadj)
      elseif(icon(3,ivadj) .eq. ivcur) then
          icon(3,ivadj) = icon(1,ivadj)
          icon(1,ivadj) = icon(2,ivadj)
          icon(5,ivadj) = icon(6,ivadj)
      else
          stop 2210
      endif
      icon(2,ivadj) = ivcur
      icon(4,ivadj) = isadj
      icon(6,ivadj) = iscur
  150 continue
      ishat = icon(5,ivadj)
c
      if(ishat.eq.isite.or.ishat.eq.isadj.or.ishat.eq.iscur) stop 2220
c
c     compute center of circumcircle for triangle isite-isadj-iscur
c
      delx1 = x(ishat) - x(isite)
      dely1 = y(ishat) - y(isite)
      dif1 = dsqrt(delx1**2 + dely1**2)
      if(dif1 .lt. epf) go to 300
c
      delx2 = y(isadj) - y(iscur)
      dely2 = x(iscur) - x(isadj)
      dif2 = dsqrt(delx2**2 + dely2**2)
      if(dif2 .lt. epf) go to 300
c
      delx3 = x(isite) - x(isadj)
      dely3 = y(isite) - y(isadj)
      dif3 = dsqrt(delx3**2 + dely3**2)
      if(dif3 .lt. epf) go to 300
c
      delx1 = x(isite) - x(iscur)
      dely1 = y(isite) - y(iscur)
      dif1 = dsqrt(delx1**2 + dely1**2)
      if(dif1 .lt. epf) go to 300
c
      dot = delx2*delx1 + dely2*dely1
      if(dot.lt.epf .and. dot.gt.-epf) go to 300
      dot2 = dot/dif2
      if(dot2.lt.epf .and. dot2.gt.-epf) go to 300
      dat = delx1*delx3 + dely1*dely3
      zlamb = dat/dot
      vx = .5*(x(isadj)+x(iscur) + zlamb*delx2)
      vy = .5*(y(isadj)+y(iscur) + zlamb*dely2)
      call biscrc (x, y, ishat, isite, epf, vx, vy, tdist)
      if(tdist.le.-epf) go to 500
      if(tdist.ge.epf) go to 400
  300 continue
      call circsn(ix, iy, ix2, iy2, isite, isadj, iscur, ishat,
     *            mhalf, mfull, isclp, itide)
      if(itide .ge. 0) go to 500
c
c     switch diagonals of quadrilateral
c
  400 continue
      icon(1,ivcur) = icon(3,ivadj)
      icon(3,ivadj) = ivcur
      icon(2,ivadj) = icon(2,ivcur)
      icon(2,ivcur) = ivadj
      icon(6,ivcur) = ishat
      icon(4,ivadj) = isite
      ivout = icon(1,ivcur)
      if(ivout .eq. 0 .or. ivout .eq. nvmz) go to 470
      ivcar = ivcur
      ivodj = ivadj
      if(ivout .lt. 0) then
          ivout = -ivout
          ivcar = -ivcar
          ivodj = -ivodj
      endif
      if(icon(1,ivout) .eq. ivodj) then
          icon(1,ivout) = ivcar
      elseif(icon(2,ivout) .eq. ivodj) then
          icon(2,ivout) = ivcar
      elseif(icon(3,ivout) .eq. ivodj) then
          icon(3,ivout) = ivcar
      else
          stop 2260
      endif
  470 continue
      ivout = icon(2,ivadj)
      if(ivout .eq. 0 .or. ivout .eq. nvmz) go to 480
      ivodj = ivadj
      if(ivout .lt. 0) then
          ivout = -ivout
          ivodj = -ivodj
      endif
      icon(3,ivout) = ivodj
  480 continue
c
      if(is(isadj) .eq. ivadj) is(isadj) = ivcur
      if(is(iscur) .eq. ivcur) is(iscur) = ivadj
c
      ivcur = ivadj
      isadj = ishat
      go to 100
c
c     move to next triangle with isite as a vertex
c
  500 continue
      if(isadj .eq. isfin) go to 900
      ivadj = icon(3,ivcur)
      ivcur = ivadj
      iscur = isadj
      isadj = icon(5,ivcur)
      go to 100
c
  900 continue
      return
      end
*BISCRC
c
c This subroutine will compute the distance from center of circle to
c bisector line between vertex of triangle and ishat
c
      subroutine biscrc (x, y, ishat, ivrt, epf, xctr, yctr, tdist)
c
      double precision x(*), y(*), norm
      integer ishat, ivrt
      double precision epf, xctr, yctr, tdist
c
      double precision xm, ym, xu, yu, xd, yd, dif, dmax
c
c     find midpoint of edge from ishat to ivrt
c
      xm = (x(ishat) + x(ivrt)) / 2.0d0
      ym = (y(ishat) + y(ivrt)) / 2.0d0
c
c     find vector from ivrt to ishat
c
      xu = x(ishat) - x(ivrt)
      yu = y(ishat) - y(ivrt)
c 
      norm = dsqrt (xu**2 + yu**2)
      if (norm .lt. epf) stop 2270
c
      xd = xctr-xm
      yd = yctr-ym
      dif = dsqrt(xd**2 + yd**2)
c
      dmax = dmax1(norm,dif)
c
c     compute distance
c
      tdist = (xd*xu + yd*yu)/dmax
c
      return
      end
*RECONC
c
c     subroutine reconc to -
c
c     reconcile a line segment with current triangulation so as to
c     obtain a Delaunay triangulation constrained by the line segments
c     that have been reconciled with the triangulation so far
c
c     November 30, 1987
c
      subroutine reconc(x, y, ix, iy, ix2, iy2, is, icon, ifun, irun,
     *                  isun, izun, igun, iwun, ia, epf, mhalf, mfull,
     *                  isclp, nfmax, nhmax, nvmz, nvmp, nvmn,
     *                  ileft, iright, n)
c
      integer mhalf, mfull, isclp(*)
      integer nfmax, nhmax, nvmz, nvmp, nvmn, ileft, iright, n
      double precision epf
      double precision x(*), y(*)
      integer ix(*), iy(*), ix2(*), iy2(*)
      integer is(*), icon(6,*), ia(*)
      integer ifun(*), irun(*), isun(*), izun(*)
      integer igun(*), iwun(*)
c
      integer nkmax, iperp
      parameter (nkmax=30)
      integer iax(nkmax), iay(nkmax), isgax, isgay, ikax, ikay, isga
      double precision delx, dely, dif, xtemp, ytemp, dot1, dot2, dot
      double precision dlg, dlg1, dlg2
      integer iloft, ivini, ivcur, ivadj, ivpre, ivfol, isi1, isi2, isi3
      integer indx, indy, lf, lr, ian, ivpro, l, isam, mslop, mintr
      integer iam, iopt, j, ial, iap, iflg
c
c     test nvmn and nvmp values
c
      if(nvmn .ne. nvmz - 1 .or. nvmp .ne. -nvmn) stop 2310
c
c     compute endpoints parameters
c
      if(ileft.lt.1.or.iright.lt.1.or.ileft.gt.n.or.iright.gt.n)
     *   stop 2320
      if(ileft .eq. iright) go to 600
      iperp = 0
      delx = y(ileft) - y(iright)
      dely = x(iright) - x(ileft)
      dif = dsqrt(delx**2 + dely**2)
      if(dif .lt. epf) iperp = 1
      call prpvec(ix, iy, ix2, iy2, ileft, iright, mhalf, mfull, isclp,
     *            iax, isgax, ikax, iay, isgay, ikay)
      if(isgax.eq.0 .and. isgay.eq.0) stop 2330
      iloft = ileft
c
   50 continue
      ivini = is(iloft)
      ivcur = ivini
c
c     test vertices of current triangle
c
  100 continue
      if(icon(4,ivcur) .eq. iloft) then
          ivadj = icon(2,ivcur)
          ivpre = icon(3,ivcur)
          ivfol = icon(1,ivcur)
          isi2 = icon(5,ivcur)
          isi1 = icon(6,ivcur)
          indx = 3
          indy = 2
      elseif(icon(5,ivcur) .eq. iloft) then
          ivadj = icon(3,ivcur)
          ivpre = icon(1,ivcur)
          ivfol = icon(2,ivcur)
          isi2 = icon(6,ivcur)
          isi1 = icon(4,ivcur)
          indx = 1
          indy = 3
      elseif(icon(6,ivcur) .eq. iloft) then
          ivadj = icon(1,ivcur)
          ivpre = icon(2,ivcur)
          ivfol = icon(3,ivcur)
          isi2 = icon(4,ivcur)
          isi1 = icon(5,ivcur)
          indx = 2
          indy = 1
      else
          stop 2340
      endif
c
      if(iperp.eq.1) go to 110
      xtemp = x(isi2) - x(iloft)
      ytemp = y(isi2) - y(iloft)
      dlg2 = dsqrt(xtemp**2 + ytemp**2)
      if(dlg2 .lt. epf) go to 110
      dot2 = (delx * xtemp + dely * ytemp)/dif
      if(dot2.ge.epf .or. dot2.le.-epf) go to 120
  110 continue
      call distsn(ix, iy, ix2, iy2, iloft, isi2, mhalf, mfull, isclp,
     *            iax, isgax, ikax, iay, isgay, ikay, isga)
      if(isga.gt.0) then
         dot2 = 10.0d0
      elseif(isga.lt.0) then
         dot2 =-10.0d0
      else
         dot2 = 0.0d0
      endif
  120 continue
      if(dot2 .ge. epf) go to 200
c
      if(iperp.eq.1) go to 130
      xtemp = x(isi1) - x(iloft)
      ytemp = y(isi1) - y(iloft)
      dlg1 = dsqrt(xtemp**2 + ytemp**2)
      if(dlg1 .lt. epf) go to 130
      dot1 = (delx * xtemp + dely * ytemp)/dif
      if(dot1.ge.epf .or. dot1.le.-epf) go to 140
  130 continue
      call distsn(ix, iy, ix2, iy2, iloft, isi1, mhalf, mfull, isclp,
     *            iax, isgax, ikax, iay, isgay, ikay, isga)
      if(isga.gt.0) then
         dot1 = 10.0d0
      elseif(isga.lt.0) then
         dot1 =-10.0d0
      else
         dot1 = 0.0d0
      endif
  140 continue
      if(dot1 .lt. epf) go to 200
      if(dot2 .le. -epf) go to 150
      if(ivpre .lt. 0) go to 500
      if(ivpre .eq. 0) then
          icon(indx,ivcur) = nvmz
          go to 500
      endif
      icon(indx,ivcur) = -ivpre
      if(icon(1,ivpre) .eq. ivcur) then
          icon(1,ivpre) = -ivcur
      elseif(icon(2,ivpre) .eq. ivcur) then
          icon(2,ivpre) = -ivcur
      elseif(icon(3,ivpre) .eq. ivcur) then
          icon(3,ivpre) = -ivcur
      else
          stop 2350
      endif
      go to 500
  150 continue
      if(isi2 .eq. iright) stop 2355
      if(2 .gt. nfmax) stop 2360
      ifun(1) = iloft
      ifun(2) = isi2
      irun(1) = iloft
      irun(2) = isi1
      lf = 2
      lr = 2
      isun(1) = ivpre
      izun(1) = ivadj
      icon(4,ivcur) = -1
      ian = 1
      if(1 .gt. nhmax) stop 2370
      ia(1) = ivcur
      go to 300
c
  200 continue
      ivpro = ivcur
      ivcur = iabs(ivadj)
      if(ivcur .eq. ivini) stop 2380
      if(ivadj .ne. 0 .and. ivadj .ne. nvmz) go to 100
      if(dot2.gt.-epf .or. dot1.le.-epf) stop 2385
      isi2 = isi1
      if(ivadj .eq. nvmz) go to 500
      icon(indy,ivpro) = nvmz
      go to 500
c
c     find vertex for next triangle
c
  300 continue
      if(ivfol .le. 0) then
          write(*,*)' '
          write(*,*)' reconc: inappropriate crossing of line segments'
          write(*,*)' endpoints of first segment are'
          write(*,*)' ( ',x(isi1),' , ',y(isi1),' ) and ( ',
     *                    x(isi2),' , ',y(isi2),' )'
          write(*,*)' endpoints of second segment are'
          write(*,*)' ( ',x(ileft),' , ',y(ileft),' ) and ( ',
     *                    x(iright),' , ',y(iright),' )'
          write(*,*)' program terminated'
          write(*,*)' '
          stop 2390
      endif
      if(icon(4,ivfol) .eq. isi2) then
          isi3 = icon(5,ivfol)
          ivadj = icon(1,ivfol)
          ivpre = icon(3,ivfol)
      elseif(icon(5,ivfol) .eq. isi2) then
          isi3 = icon(6,ivfol)
          ivadj = icon(2,ivfol)
          ivpre = icon(1,ivfol)
      elseif(icon(6,ivfol) .eq. isi2) then
          isi3 = icon(4,ivfol)
          ivadj = icon(3,ivfol)
          ivpre = icon(2,ivfol)
      else
          stop 2400
      endif
c
      if(iperp.eq.1) go to 330
      xtemp = x(isi3) - x(iloft)
      ytemp = y(isi3) - y(iloft)
      dlg = dsqrt(xtemp**2 + ytemp**2)
      if(dlg .lt. epf) go to 330
      dot = (delx * xtemp + dely * ytemp)/dif
      if(dot.ge.epf .or. dot.le.-epf) go to 340
  330 continue
      call distsn(ix, iy, ix2, iy2, iloft, isi3, mhalf, mfull, isclp,
     *            iax, isgax, ikax, iay, isgay, ikay, isga)
      if(isga.gt.0) then
         dot = 10.0d0
      elseif(isga.lt.0) then
         dot =-10.0d0
      else
         dot = 0.0d0
      endif
  340 continue
      icon(4,ivfol) = -1
      ian = ian + 1
      if(ian .gt. nhmax) stop 2410
      ia(ian) = ivfol
      if((dot .le. -epf .or.  dot .ge. epf) .and. isi3 .eq. iright)
     *   stop 2415
      if(dot .gt. -epf .and. dot .lt. epf) then
          lf = lf + 1
          lr = lr + 1
          if(lf .gt. nfmax) stop 2420
          if(lr .gt. nfmax) stop 2430
          ifun(lf) = isi3
          irun(lr) = isi3
          isun(lf-1) = ivpre
          izun(lr-1) = ivadj
          isi2 = isi3
          go to 400
      elseif(dot .ge. epf) then
          lr = lr + 1
          if(lr .gt. nfmax) stop 2440
          irun(lr) = isi3
          izun(lr-1) = ivadj
          ivfol = ivpre
          isi1 = isi3
      else
          lf = lf + 1
          if(lf .gt. nfmax) stop 2450
          ifun(lf) = isi3
          isun(lf-1) = ivpre
          ivfol = ivadj
          isi2 = isi3
      endif
      go to 300
c
c     retriangulate
c
  400 continue
c
c     test for eliminated adjacent triangles
c
      do 410 l = 1, lf - 1
          isam = isun(l)
          if(isam .eq. 0 .or. isam .eq. nvmz) go to 410
          if(icon(4,iabs(isam)) .ne. -1) go to 410
          if(isam .gt. 0) then
              isun(l) = nvmp
          else
              isun(l) = nvmn
          endif
  410 continue
      do 420 l = 1, lr - 1
          isam = izun(l)
          if(isam .eq. 0 .or. isam .eq. nvmz) go to 420
          if(icon(4,iabs(isam)) .ne. -1) go to 420
          if(isam .gt. 0) then
              izun(l) = nvmp
          else
              izun(l) = nvmn
          endif
  420 continue
c
c     triangulate 'right' edge star-shaped polygon
c
      mslop = 1
      mintr = 0
      iam = 0
      iopt = 3
c
      call edgstr(x, y, ix, iy, ix2, iy2, icon, ifun, igun, iwun,
     *            ia, epf, mhalf, mfull, isclp, lf,
     *            ian, iam, mslop, mintr, j, iopt)
c
      if(j .ne. 2) stop 2460
      ial = iwun(1)
      iap = ial
c
c     reconcile new triangles with rest of triangulation
c
      iflg = 0
      call mattri(is, icon, ifun, isun, ial, lf, nvmz, nvmp, nvmn, iflg)
c
c     triangulate 'left' edge star-shaped polygon
c
      lf = lr
      do 480 l = 1, lr - 1
          ifun(l) = irun(lr + 1 - l)
          isun(l) = izun(lr - l)
  480 continue
      ifun(lr) = irun(1)
      mslop = -1
      mintr = ian + iam + 1
c
      call edgstr(x, y, ix, iy, ix2, iy2, icon, ifun, igun, iwun,
     *            ia, epf, mhalf, mfull, isclp, lf,
     *            ian, iam, mslop, mintr, j, iopt)
c
      if(j .ne. 2) stop 2470
      ial = iwun(1)
      icon(2,ial) = -iap
      icon(2,iap) = -ial
c
c     reconcile new triangles with rest of triangulation
c
      iflg = 1
      call mattri(is, icon, ifun, isun, ial, lf, nvmz, nvmp, nvmn, iflg)
c
c     reinitialize to reconcile next portion of line segment
c
  500 continue
      if(isi2 .eq. iright) go to 600
      iloft = isi2
      go to 50
c
  600 continue
      return
      end
*EDGSTR
c
c     subroutine edgsrt to -
c
c     compute a (Delaunay) triangulation for a subset of the boundary
c     of one or more adjacent edge star-shaped simple polygons
c     constrained by their boundaries
c
c     December 4, 1987
c
      subroutine edgstr(x, y, ix, iy, ix2, iy2, icon, ifun, igun, iwun,
     *                  ia, epf, mhalf, mfull, isclp, lf,
     *                  ian, iam, mslop, mintr, j, iopt)
c
      integer mhalf, mfull, isclp(*)
      integer lf, ian, iam, mslop, mintr, j, iopt
      double precision epf
      double precision x(*), y(*)
      integer ix(*), iy(*), ix2(*), iy2(*)
      integer icon(6,*), ia(*)
      integer ifun(*), igun(*), iwun(*)
c
      integer i, iq1, iq2, iq3, ial, isga
      double precision dulx, duly, dif, dlst, xtemp, ytemp, dot
c
      iwun(1) = 0
      igun(1) = ifun(1)
      igun(2) = ifun(2)
      i = 2
      j = 2
c
  100 continue
      if(i .eq. lf) go to 400
      iwun(j) = 0
      i = i + 1
      if(i .eq. lf .and. iopt .eq. 4) iopt = 1
      j = j + 1
      igun(j) = ifun(i)
      iq1 = igun(j-2)
      iq2 = igun(j-1)
      iq3 = igun(j)
c
  200 continue
      if(i .eq. lf .and. iopt .gt. 0) go to 300
      dulx = y(iq1) - y(iq2)
      duly = x(iq2) - x(iq1)
      dif = dsqrt(dulx**2 + duly**2)
      if(dif .lt. epf) go to 220
      xtemp = x(iq3) - x(iq1)
      ytemp = y(iq3) - y(iq1)
      dlst = dsqrt(xtemp**2 + ytemp**2)
      if(dlst .lt. epf) go to 220
      dot = (dulx * xtemp + duly * ytemp)/dif
      if(dot.ge.epf .or. dot.le.-epf) go to 240
  220 continue
      call innsgn(ix, iy, ix2, iy2, iq1, iq2, iq3, mhalf, mfull,
     *            isclp, isga)
      if(isga.gt.0) then
         dot = 10.0d0
      elseif(isga.lt.0) then
         dot =-10.0d0
      else
         dot = 0.0d0
      endif
  240 continue
      if(dot .lt. epf) go to 100
c
  300 continue
      iam = iam + 1
      if(iam .gt. ian) stop 2500
      ial = ia(mslop*iam+mintr)
      icon(3,ial) = iwun(j-2)
      icon(1,ial) = iwun(j-1)
      icon(2,ial) = 0
      icon(4,ial) = iq1
      icon(5,ial) = iq2
      icon(6,ial) = iq3
      if(iwun(j-2) .ne. 0) icon(2,iwun(j-2)) = ial
      if(iwun(j-1) .ne. 0) icon(2,iwun(j-1)) = ial
c
c     update triangulation to obtain a Delaunay triangulation of the
c     union of the triangles obtained so far inside the polygon
c
      if(iopt .gt. 1) call updtri(x, y, ix, iy, ix2, iy2, icon, epf,
     *                mhalf, mfull, isclp, ial, iq1, iq2, iq3)
c
      j = j - 1
      iwun(j-1) = ial
      igun(j) = iq3
c
      if(j .eq. 2) go to 100
      iq1 = igun(j-2)
      iq2 = igun(j-1)
      go to 200
c
  400 continue
c
      return
      end
*MATTRI
c
c     subroutine mattri to -
c
c     reconcile triangulation of an edge star-shaped polygon with
c     rest of triangulation
c
c     December 8, 1987
c
      subroutine mattri(is, icon, ifun, isun, ial, lf, nvmz, nvmp,
     *                  nvmn, iflg)
c
      integer ial, lf, nvmz, nvmp, nvmn, iflg
      integer is(*), icon(6,*), ifun(*), isun(*)
c
      integer ivcur, l, isite, ivadj, index, ivout, ivabs, ivpre
      integer ind1, ind2, isadj, isnxt, ivnow, l1, l2, ivini
c
c     reconcile triangulation with triangles outside polygon
c
      ivcur = ial
      do 400 l = 1, lf-1
          isite = ifun(l)
  100     continue
          if(icon(4,ivcur) .eq. isite) then
              ivadj = icon(3,ivcur)
              index = 3
          elseif(icon(5,ivcur) .eq. isite) then
              ivadj = icon(1,ivcur)
              index = 1
          elseif(icon(6,ivcur) .eq. isite) then
              ivadj = icon(2,ivcur)
              index = 2
          else
              stop 2600
          endif
          if(ivadj .eq. 0) go to 200
          ivcur = ivadj
          go to 100
c
  200     continue
          is(isite) = ivcur
          ivout = isun(l)
          if(ivout .eq. nvmp .or. ivout .eq. nvmn) go to 300
          isun(l) = 0
          icon(index,ivcur) = ivout
          if(ivout .eq. 0 .or. ivout .eq. nvmz) go to 400
          ivadj = ivcur
          if(ivout .lt. 0) ivadj = -ivcur
          ivabs = iabs(ivout)
          if(icon(4,ivabs) .eq. isite) then
              icon(2,ivabs) = ivadj
          elseif(icon(5,ivabs) .eq. isite) then
              icon(3,ivabs) = ivadj
          elseif(icon(6,ivabs) .eq. isite) then
              icon(1,ivabs) = ivadj
          else
              stop 2610
          endif
          go to 400
c
  300     continue
          isun(l) = ivcur
          if(ivout .eq. nvmn) isun(l) = -ivcur
c
  400 continue
c
c     reconcile triangulation with triangles inside polygon
c
      do 600 l = 1, lf-1
          if(isun(l) .eq. 0) go to 600
          ivout = isun(l)
          ivcur = iabs(ivout)
          isite = ifun(l)
          if(icon(4,ivcur) .eq. isite) then
              ivpre = icon(3,ivcur)
              ind1 = 3
          elseif(icon(5,ivcur) .eq. isite) then
              ivpre = icon(1,ivcur)
              ind1 = 1
          elseif(icon(6,ivcur) .eq. isite) then
              ivpre = icon(2,ivcur)
              ind1 = 2
          else
              stop 2620
          endif
          if(ivpre .ne. 0) stop 2630
          if(ivcur .ne. is(isite)) go to 600
          isadj = ifun(l+1)
c
          ivadj = ivcur
  500     continue
          ivnow = iabs(ivadj)
          if(icon(4,ivnow) .eq. isite) then
              ivadj = icon(2,ivnow)
              ind2 = 2
              isnxt = icon(6,ivnow)
          elseif(icon(5,ivnow) .eq. isite) then
              ivadj = icon(3,ivnow)
              ind2 = 3
              isnxt = icon(4,ivnow)
          elseif(icon(6,ivnow) .eq. isite) then
              ivadj = icon(1,ivnow)
              ind2 = 1
              isnxt = icon(5,ivnow)
          else
              stop 2640
          endif
          if(ivadj .ne. 0) go to 500
          if(isadj .ne. isnxt) stop 2650
          icon(ind1,ivcur) = ivnow
          if(ivout .lt. 0) icon(ind1,ivcur) = -ivnow
          icon(ind2,ivnow) = ivout
c
  600 continue
c
c     update array is
c
      if(iflg .eq. 0) then
          l1 = 2
          l2 = lf - 1
      else
          l1 = 1
          l2 = lf
      endif
      do 900 l = l1, l2
          isite = ifun(l)
          ivcur = is(isite)
          ivini = ivcur
  700     continue
          if(icon(4,ivcur) .eq. isite) then
              ivadj = icon(3,ivcur)
          elseif(icon(5,ivcur) .eq. isite) then
              ivadj = icon(1,ivcur)
          elseif(icon(6,ivcur) .eq. isite) then
              ivadj = icon(2,ivcur)
          else
              stop 2660
          endif
          if(ivadj .eq. 0 .or. ivadj .eq. nvmz) go to 800
          if(ivadj .lt. 0) is(isite) = ivcur
          ivcur = iabs(ivadj)
          if(ivcur .eq. ivini) go to 900
          go to 700
  800     continue
          is(isite) = ivcur
  900 continue
c
      return
      end
*UPDTRI
c
c     subroutine updtri to -
c
c     update a triangulation of a generalized multiply-connected
c     polygonal region to obtain a Delaunay triangulation
c     of the region
c
c     December 10, 1987
c
      subroutine updtri(x, y, ix, iy, ix2, iy2, icon, epf,
     *                  mhalf, mfull, isclp, ial, iq1, iq2, iq3)
c
      integer mhalf, mfull, isclp(*)
      integer ial, iq1, iq2, iq3
      double precision epf
      double precision x(*), y(*)
      integer ix(*), iy(*), ix2(*), iy2(*)
      integer icon(6,*)
c
      integer ivcur, isadj, iscur, isite, ivadj, ishat, itide
      double precision delx2, dely2, delx3, dely3, delx1, dely1
      double precision dif1, dif2, dif3, tdist
      double precision zlamb, vx, vy, dat, dot, dot2
c
c     initialize
c
      ivcur = ial
      isadj = iq1
      iscur = iq2
      isite = iq3
c
  100 continue
      if(icon(3,ivcur) .eq. 0) go to 500
      ivadj = icon(3,ivcur)
      ishat = icon(5,ivadj)
c
c     compute center of circumcircle for triangle isadj-iscur-iq3
c
      delx1 = x(ishat) - x(isite)
      dely1 = y(ishat) - y(isite)
      dif1 = dsqrt(delx1**2 + dely1**2)
      if(dif1 .lt. epf) go to 300
c
      delx2 = y(isadj) - y(iscur)
      dely2 = x(iscur) - x(isadj)
      dif2 = dsqrt(delx2**2 + dely2**2)
      if(dif2 .lt. epf) go to 300
c
      delx3 = x(isite) - x(isadj)
      dely3 = y(isite) - y(isadj)
      dif3 = dsqrt(delx3**2 + dely3**2)
      if(dif3 .lt. epf) go to 300
c
      delx1 = x(isite) - x(iscur)
      dely1 = y(isite) - y(iscur)
      dif1 = dsqrt(delx1**2 + dely1**2)
      if(dif1 .lt. epf) go to 300
c
      dot = delx2*delx1 + dely2*dely1
      if(dot.lt.epf .and. dot.gt.-epf) go to 300
      dot2 = dot/dif2
      if(dot2.lt.epf .and. dot2.gt.-epf) go to 300
      dat = delx1*delx3 + dely1*dely3
      zlamb = dat/dot
      vx = .5*(x(isadj)+x(iscur) + zlamb*delx2)
      vy = .5*(y(isadj)+y(iscur) + zlamb*dely2)
      call biscrc (x, y, ishat, isite, epf, vx, vy, tdist)
      if(tdist.le.-epf) go to 500
      if(tdist.ge.epf) go to 400
  300 continue
      call circsn(ix, iy, ix2, iy2, isite, isadj, iscur, ishat,
     *            mhalf, mfull, isclp, itide)
      if(itide .ge. 0) go to 500
c
c     switch diagonals of quadrilateral
c
  400 continue
      icon(3,ivcur) = icon(3,ivadj)
      icon(3,ivadj) = icon(1,ivadj)
      icon(1,ivadj) = icon(1,ivcur)
      icon(1,ivcur) = ivadj
      icon(5,ivcur) = ishat
      icon(4,ivadj) = ishat
      icon(5,ivadj) = iscur
      icon(6,ivadj) = iq3
      if(icon(3,ivcur) .ne. 0) icon(2,icon(3,ivcur)) = ivcur
      if(icon(1,ivadj) .ne. 0) icon(2,icon(1,ivadj)) = ivadj
c
      ivcur = ivadj
      isadj = ishat
      go to 100
c
c     move to next triangle with iq3 as a vertex
c
  500 continue
      if(isadj .eq. iq1) go to 600
      ivadj = icon(2,ivcur)
      ivcur = ivadj
      iscur = isadj
      isadj = icon(4,ivcur)
      go to 100
c
  600 continue
c
      return
      end
*DELELI
c
c     subroutine deleli to -
c
c     eliminate a vertex from a Delaunay triangulation so as to obtain
c     a Delaunay triangulation that does not include the vertex
c     (vertex can not be in an inserted edge)
c
c     February 28, 1991
c
      subroutine deleli(x, y, ix, iy, ix2, iy2, is, icon, ifun, irun,
     *                  isun, izun, igun, iwun, ia, iave, iaze, epf,
     *                  mhalf, mfull, isclp, nfmax, nhmax, njmax,
     *                  nkmax, nvmz, nvmp, nvmn, isite, n, isucs,
     *                  isnew, iabe, iase)
c
      integer mhalf, mfull, isclp(*)
      integer nfmax, nhmax, njmax, nkmax, nvmz, nvmp, nvmn
      integer isite, n, isucs, isnew, iabe, iase
      double precision epf
      double precision x(*), y(*)
      integer ix(*), iy(*), ix2(*), iy2(*)
      integer is(*), icon(6,*), ia(*)
      integer ifun(*), irun(*), isun(*), izun(*)
      integer igun(*), iwun(*), iave(*), iaze(*)
c
      integer mkmax, iperp
      parameter (mkmax=30)
      integer iax(mkmax), iay(mkmax), isgax, isgay, ikax, ikay, isga
      integer ivini, ivpre, ispre, ivout, ivadj, isadj, lf, lr, ian, i
      integer ivcur, ivtro, l, isam, mslop, mintr, iam, iopt, jr, ial
      integer iflg, ivabs, iscor, indx, lg, j, jf, isnow, iscur, iel
      double precision delx, dely, dif, dlg, xtemp, ytemp, dot
c
c     test nvmn and nvmp values
c
      if(nvmn .ne. nvmz - 1 .or. nvmp .ne. -nvmn) stop 2800
c
c     obtain 'right' side information
c
      isucs = 1
      if(n - iase + 1 .le. 3) then
          isucs = -2
          go to 2300
      endif
      if(isite .lt. 1 .or. isite .gt. n) stop 2810
      ivini = is(isite)
      if(icon(4,ivini) .eq. isite) then
          ivpre = icon(1,ivini)
          ispre = icon(5,ivini)
          ivout = icon(3,ivini)
          ivadj = icon(2,ivini)
          isadj = icon(6,ivini)
      elseif(icon(5,ivini) .eq. isite) then
          ivpre = icon(2,ivini)
          ispre = icon(6,ivini)
          ivout = icon(1,ivini)
          ivadj = icon(3,ivini)
          isadj = icon(4,ivini)
      elseif(icon(6,ivini) .eq. isite) then
          ivpre = icon(3,ivini)
          ispre = icon(4,ivini)
          ivout = icon(2,ivini)
          ivadj = icon(1,ivini)
          isadj = icon(5,ivini)
      else
          stop 2820
      endif
c
      lf = 0
      if(ivout .lt. 0) then
          isucs = 0
          go to 2300
      endif
c
c     compute 'splitting' vector
c
      iperp = 0
      delx = y(ispre) - y(isite)
      dely = x(isite) - x(ispre)
      dif = dsqrt(delx ** 2 + dely ** 2)
      if(dif .lt. epf) iperp = 1
      call prpvec(ix, iy, ix2, iy2, ispre, isite, mhalf, mfull, isclp,
     *            iax, isgax, ikax, iay, isgay, ikay)
      if(isgax.eq.0 .and. isgay.eq.0) stop 2830
c
c     eliminate first triangle
c
      lr = 2
      if(lr .gt. nfmax) stop 2840
      irun(1) = ispre
      irun(2) = isadj
      izun(1) = ivpre
      icon(4,ivini) = -1
      ian = 1
      if(ian .gt. nhmax) stop 2850
      ia(1) = ivini
c
c     analyze next triangle if any
c
  200 continue
      if(ivadj .lt. 0) then
          isucs = 0
          go to 2300
      endif
      if(ivadj .eq. 0) go to 700
      ivcur = ivadj
      if(icon(4,ivcur) .eq. isite) then
          ivpre = icon(1,ivcur)
          ivadj = icon(2,ivcur)
          isadj = icon(6,ivcur)
      elseif(icon(5,ivcur) .eq. isite) then
          ivpre = icon(2,ivcur)
          ivadj = icon(3,ivcur)
          isadj = icon(4,ivcur)
      elseif(icon(6,ivcur) .eq. isite) then
          ivpre = icon(3,ivcur)
          ivadj = icon(1,ivcur)
          isadj = icon(5,ivcur)
      else
          stop 2860
      endif
c
c     ascertain which side new vertex is on
c
      if(isadj .eq. ispre) go to 400
      if(iperp.eq.1) go to 210
      xtemp = x(isadj) - x(ispre)
      ytemp = y(isadj) - y(ispre)
      dlg = dsqrt(xtemp**2 + ytemp**2)
      if(dlg .lt. epf) go to 210
      dot = (delx * xtemp + dely * ytemp)/dif
      if(dot.ge.epf .or. dot.le.-epf) go to 220
  210 continue
      call distsn(ix, iy, ix2, iy2, ispre, isadj, mhalf, mfull, isclp,
     *            iax, isgax, ikax, iay, isgay, ikay, isga)
      if(isga.gt.0) then
         dot = 10.0d0
      elseif(isga.lt.0) then
         dot =-10.0d0
      else
         dot = 0.0d0
      endif
  220 continue
      if(dot .ge. epf) go to 400
c
c     add next vertex for 'right' side
c
      if(lf .ne. 0) stop 2865
      lr = lr + 1
      if(lr .gt. nfmax) stop 2870
      irun(lr) = isadj
      izun(lr-1) = ivpre
      icon(4,ivcur) = -1
      ian = ian + 1
      if(ian .gt. nhmax) stop 2880
      ia(ian) = ivcur
      go to 200
c
c     obtain 'left' side information
c
  400 continue
      if(ivout.eq.0) stop 2890
      lf = lf + 1
      if(lf .gt. nfmax) stop 2900
      ifun(lf) = isadj
      isun(lf) = ivpre
      icon(4,ivcur) = -1
      ian = ian + 1
      if(ian .gt. nhmax) stop 2910
      ia(ian) = ivcur
      if(isadj .eq. ispre) go to 500
      go to 200
  500 continue
      if(lf .lt. 2) stop 2920
      ivtro = isun(1)
      do 600 l = 2, lf
          isun(l-1) = isun(l)
  600 continue
c
c     test for eliminated adjacent triangles
c
      do 900 l = 1, lf -1
          isam = isun(l)
          if(isam .eq. 0 .or. isam .eq. nvmz) go to 900
          if(icon(4,iabs(isam)) .ne. -1) go to 900
          if(isam .gt. 0) then
              isun(l) = nvmp
          else
              isun(l) = nvmn
          endif
  900 continue
c
  700 continue
      do 800 l = 1, lr - 1
          isam = izun(l)
          if(isam .eq. 0 .or. isam .eq. nvmz) go to 800
          if(icon(4,iabs(isam)) .ne. -1) go to 800
          if(isam .gt. 0) then
              izun(l) = nvmp
          else
              izun(l) = nvmn
          endif
  800 continue
c
      mslop = 1
      mintr = 0
      iam = 0
      isnew = irun(lr)
      if(ivout .eq. 0) go to 1000
      if(lf.eq.0) stop 2930
      go to 1200
c
c     point to be eliminated is on convex hull
c
 1000 continue
      lr = lr + 1
      if(lr .gt. nfmax) stop 2940
      irun(lr) = isite
      izun(lr-1) = 0
c
c     triangulate 'right' side in a Delaunay way
c
      iopt = 4
c
      call edgstr(x, y, ix, iy, ix2, iy2, icon, irun, igun, iwun,
     *            ia, epf, mhalf, mfull, isclp, lr,
     *            ian, iam, mslop, mintr, jr, iopt)
c
      ial = iwun(1)
      is(isite) = ial
c
c     reconcile with rest of triangulation
c
      iflg = 1
      call mattri(is, icon, irun, izun, ial, lr, nvmz, nvmp, nvmn, iflg)
c
c     eliminate superfluous triangles
c
 1100 continue
      ivcur = icon(3,ial)
      if(ivcur.eq.0) then
         write(*,*)' '
         write(*,*)'After deletion of a point other points are colinear'
         write(*,*)'Program does not handle this unlikely event'
         write(*,*)'Program terminated'
         write(*,*)' '
         stop
      endif
      ivabs = iabs(ivcur)
      iscor = icon(5,ial)
      if(icon(4,ivabs) .eq. iscor) then
          indx = 3
      elseif(icon(5,ivabs) .eq. iscor) then
          indx = 1
      elseif(icon(6,ivabs) .eq. iscor) then
          indx = 2
      else
          stop 2950
      endif
      icon(indx,ivabs) = 0
      if(ivcur .lt. 0) icon(indx,ivabs) = nvmz
      is(iscor) = ivabs
      icon(4,ial) = -1
      iabe = iabe + 1
      if(iabe .gt. njmax) stop 2960
      iave(iabe) = ial
      ial = icon(1,ial)
      if(ial .ne. 0) go to 1100
      go to 2200
c
c     point to be eliminated is not on convex hull
c
 1200 continue
c
c     triangulate 'right' side (not necessarily Delaunay)
c
      iopt = 0
c
      call edgstr(x, y, ix, iy, ix2, iy2, icon, irun, igun, iwun,
     *            ia, epf, mhalf, mfull, isclp, lr,
     *            ian, iam, mslop, mintr, jr, iopt)
c
c     update 'left' side information
c
      lg = lf
      do 1300 j = 2, jr
          lf = lf + 1
          if(lf .gt. nfmax) stop 2970
          ifun(lf) = igun(j)
          isun(lf-1) = iwun(j-1)
 1300 continue
c
c     triangulate 'left' side in a Delaunay fashion
c
      iopt = 2
c
      call edgstr(x, y, ix, iy, ix2, iy2, icon, ifun, igun, iwun,
     *            ia, epf, mhalf, mfull, isclp, lf,
     *            ian, iam, mslop, mintr, jf, iopt)
c
      if(jf .ne. 2) stop 2980
      ial = iwun(1)
c
c     merge the two sides
c
      lf = lg
      do 1400 j = 1, jr - 1
          iwun(j) = isun(lf+j-1)
 1400 continue
c
      j = jr
      ivcur = ial
 1500 continue
      ivadj = icon(1,ivcur)
      if(ivadj .eq. 0) go to 1600
      ivcur = ivadj
      go to 1500
 1600 continue
      j = j - 1
      ivadj = iwun(j)
      icon(1,ivcur) = ivadj
      if(ivadj .ne. 0) icon(2,ivadj) = ivcur
      ivcur = icon(3,ivcur)
      if(j .ne. 1) go to 1500
c
c     add each triangle to 'left' polygon and optimize
c
      j = jr - 1
 1700 continue
      if(j .eq. 0) go to 1900
      ivcur = iwun(j)
      if(ivcur .eq. 0) go to 1800
      iwun(j) = icon(3,ivcur)
      j = j + 1
      iwun(j) = icon(1,ivcur)
      indx = 2
      isnow = icon(5,ivcur)
      isadj = icon(6,ivcur)
      iscur = icon(4,ivcur)
c
      call addtri(x, y, ix, iy, ix2, iy2, icon, epf, mhalf, mfull,
     *            isclp, ivcur, indx, isnow, isadj, iscur)
c
      go to 1700
 1800 continue
      j = j - 1
      go to 1700
 1900 continue
c
c     reconcile everything with rest of triangulation
c
      if(icon(4,ial) .ne. ifun(1) .or. icon(6,ial) .ne. isnew) 
     *   stop 2990
      icon(2,ial) = ivtro
      if(ivtro .eq. 0 .or. ivtro .eq. nvmz) go to 2000
      iel = ial
      if(ivtro .lt. 0) then
          ivtro = -ivtro
          iel = -iel
      endif
      if(icon(4,ivtro) .eq. isnew) then
          icon(2,ivtro) = iel
      elseif(icon(5,ivtro) .eq. isnew) then
          icon(3,ivtro) = iel
      elseif(icon(6,ivtro) .eq. isnew) then
          icon(1,ivtro) = iel
      else
          stop 3000
      endif
c
 2000 continue
      is(isnew) = ial
      do 2100 l = 2, lr
          lf = lf + 1
          if(lf .gt. nfmax) stop 3010
          ifun(lf) = irun(l)
          isun(lf-1) = izun(l-1)
 2100 continue
c
      iflg = 1
      call mattri(is, icon, ifun, isun, ial, lf, nvmz, nvmp, nvmn, iflg)
c
c     record eliminated triangles
c
      if(iam .ge. ian) stop 3020
      iam = iam + 1
      do 2150 i = iam, ian
          iabe = iabe + 1
          if(iabe .gt. njmax) stop 3030
          iave(iabe) = ia(i)
 2150 continue
c
c     record eliminated point
c
 2200 continue
      is(isite) = 0
      iase = iase + 1
      if(iase .gt. nkmax) stop 3040
      iaze(iase) = isite
c
 2300 continue
      return
      end
*ADDTRI
c
c     subroutine addtri to -
c
c     add a triangle to a Delaunay triangulation of a polygon and get
c     a new Delaunay triangulation of the new polygon if the new
c     triangle has only one side in common with the polygon
c
c     March 6, 1991
c
      subroutine addtri(x, y, ix, iy, ix2, iy2, icon, epf, mhalf,
     *                  mfull, isclp, ivcur, indx, isite, isadj, iscur)
c
      integer mhalf, mfull, isclp(*)
      integer ivcur, indx, isite, isadj, iscur
      double precision epf
      double precision x(*), y(*)
      integer ix(*), iy(*), ix2(*), iy2(*)
      integer icon(6,*)
      integer isfin, ivadj, indy, ishat, itide
      double precision delx1, dely1, delx2, dely2, delx3, dely3
      double precision dif1, dif2, dif3, tdist
      double precision dat, dot, dot2, zlamb, vx, vy
      integer iv1, iv2
c
c     initialize
c
      isfin = isadj
c
  100 continue
      ivadj = icon(indx,ivcur)
      if(ivadj .eq. 0) go to 500
c
      if(icon(1,ivadj) .eq. ivcur) then
          indy = 4
      elseif(icon(2,ivadj) .eq. ivcur) then
          indy = 5
      elseif(icon(3,ivadj) .eq. ivcur) then
          indy = 6
      else
          stop 3050
      endif
      ishat = icon(indy,ivadj)
c
c     compute center of circumcircle for triangle isite-isadj-iscur
c
      delx1 = x(ishat) - x(isite)
      dely1 = y(ishat) - y(isite)
      dif1 = dsqrt(delx1**2 + dely1**2)
      if(dif1 .lt. epf) go to 300
c
      delx2 = y(isadj) - y(iscur)
      dely2 = x(iscur) - x(isadj)
      dif2 = dsqrt(delx2**2 + dely2**2)
      if(dif2 .lt. epf) go to 300
c
      delx3 = x(isite) - x(isadj)
      dely3 = y(isite) - y(isadj)
      dif3 = dsqrt(delx3**2 + dely3**2)
      if(dif3 .lt. epf) go to 300
c
      delx1 = x(isite) - x(iscur)
      dely1 = y(isite) - y(iscur)
      dif1 = dsqrt(delx1**2 + dely1**2)
      if(dif1 .lt. epf) go to 300
c
      dot = delx2*delx1 + dely2*dely1
      if(dot.lt.epf .and. dot.gt.-epf) go to 300
      dot2 = dot/dif2
      if(dot2.lt.epf .and. dot2.gt.-epf) go to 300
      dat = delx1*delx3 + dely1*dely3
      zlamb = dat/dot
      vx = .5*(x(isadj)+x(iscur) + zlamb*delx2)
      vy = .5*(y(isadj)+y(iscur) + zlamb*dely2)
      call biscrc (x, y, ishat, isite, epf, vx, vy, tdist)
      if(tdist.le.-epf) go to 500
      if(tdist.ge.epf) go to 400
  300 continue
      call circsn(ix, iy, ix2, iy2, isite, isadj, iscur, ishat,
     *            mhalf, mfull, isclp, itide)
      if(itide .ge. 0) go to 500
c
c     switch diagonals of quadrilateral
c
  400 continue
      iv1 = ivadj
      iv2 = ivcur
      if(indy .eq. 5) then
          iv1 = ivcur
          iv2 = ivadj
      endif
      if(indx .ne. 1 .and. indy .ne. 4) then
          icon(3,iv1) = icon(3,iv2)
          icon(3,iv2) = icon(1,iv2)
          icon(1,iv2) = icon(1,iv1)
          icon(1,iv1) = iv2
          icon(5,iv1) = icon(5,iv2)
          icon(4,iv2) = icon(5,iv2)
          icon(5,iv2) = icon(6,iv2)
          icon(6,iv2) = icon(6,iv1)
          if(icon(3,iv1) .ne. 0) icon(2,icon(3,iv1)) = iv1
          if(icon(1,iv2) .ne. 0) icon(2,icon(1,iv2)) = iv2
          ivcur = ivadj
      else
          icon(1,iv1) = icon(1,iv2)
          icon(1,iv2) = icon(3,iv2)
          icon(3,iv2) = icon(3,iv1)
          icon(3,iv1) = iv2
          icon(5,iv1) = icon(5,iv2)
          icon(6,iv2) = icon(5,iv2)
          icon(5,iv2) = icon(4,iv2)
          icon(4,iv2) = icon(4,iv1)
          if(icon(1,iv1) .ne. 0) icon(2,icon(1,iv1)) = iv1
          if(icon(3,iv2) .ne. 0) icon(2,icon(3,iv2)) = iv2
      endif
      if(indy .eq. 4) indx = 3
      isadj = ishat
      go to 100
c
c     move to next triangle with isite as a vertex
c
  500 continue
      if(isadj .eq. isfin) go to 600
      if(indx .eq. 1) then
          ivadj = icon(3,ivcur)
      elseif(indx .eq. 2) then
          ivadj = icon(1,ivcur)
      else
          ivadj = icon(2,ivcur)
      endif
      ivcur = ivadj
      iscur = isadj
      if(icon(4,ivcur) .eq. isite) then
          indx = 1
          isadj = icon(5,ivcur)
      elseif(icon(5,ivcur) .eq. isite) then
          indx = 2
          isadj = icon(6,ivcur)
      elseif(icon(6,ivcur) .eq. isite) then
          indx = 3
          isadj = icon(4,ivcur)
      else
          stop 3080
      endif
      go to 100
c
  600 continue
      return
      end
*COMPRS
c
c     subroutine comprs to -
c
c     compress array icon so that there are no gaps in it
c
c     March 15, 1991
c
      subroutine comprs(is, icon, id, n, ivnxt, nmax, nvmax, nvmz)
c
      integer n, ivnxt, nmax, nvmax, nvmz
      integer is(*), icon(6,*), id(*)
      integer ivnot, i, j, ivcur
c
c     test n, ivnxt, nvmz
c
      if(n .gt. nmax .or. ivnxt .gt. nvmax) stop 3090
c
      if(nvmz .ne. -nvmax - 1) stop 3100
c
c     compress icon
c
      ivnot = 0
      do 200 i = 1, ivnxt
          if(icon(4,i) .eq. -1) go to 200
          ivnot = ivnot + 1
          id(i) = ivnot
          do 100 j = 1, 6
              icon(j,ivnot) = icon(j,i)
  100     continue
  200 continue
      ivnxt = ivnot
c
c     update icon for triangles
c
      do 400 i = 1, ivnxt
          do 300 j = 1, 3
              ivcur = icon(j,i)
              if(ivcur .eq. 0 .or. ivcur .eq. nvmz) go to 300
              if(ivcur .gt. 0) then
                  icon(j,i) = id(ivcur)
              else
                  ivcur = -ivcur
                  icon(j,i) = -id(ivcur)
              endif
  300     continue
  400 continue
c
c     update is for triangles
c
      do 500 i = 1, n
          if(is(i) .le. 0) go to 500
          is(i) = id(is(i))
  500 continue
c
      return
      end
*CMPRS2
c
c     subroutine cmprs2 to -
c
c     compress arrays icon, vx, vy, so that only Voronoi vertices
c     inside or on boundary of rectangle are kept
c
c     April 2, 2013
c
      subroutine cmprs2(is, icon, id, vx, vy, ifl, nr, n, ivnxt,
     *                  nmax, nvmax)
c
      double precision vx(*), vy(*)
      integer nr, n, ivnxt, nmax, nvmax
      integer is(*), icon(6,*), id(*), ifl(*)
      integer ivnot, i, j, ivcur, ivini, ivadj
c
c     test n, nr, ivnxt, nvmz
c
      if(n .ge. nr) stop 3105
      if(nr .gt. nmax .or. ivnxt .gt. nvmax) stop 3110
c
c     redefine is array properly: if the boundary of cell for i
c     does not intersect edge of rectangle it remains the same;
c     if it does and has a vertex in the interior of rectangle
c     it becomes vertex before first vertex in the interior
c     of rectangle in a counterclockwise direction; otherwise
c     it remains the same, a vertex on edge of rectangle
c
      do 180 i = 1, n
         ivini = is(i)
         if(ivini.le.0) go to 180
         ivcur = ivini
c
c     test whether boundary of cell intersects edge of rectangle
c
   20    continue
         if(icon(4,ivcur).gt.n .or. icon(5,ivcur).gt.n .or.
     *      icon(6,ivcur).gt.n) go to 40
         if(icon(4,ivcur).eq.i) then
            ivadj = icon(2,ivcur)
         elseif(icon(5,ivcur).eq.i) then
            ivadj = icon(3,ivcur)
         elseif(icon(6,ivcur).eq.i) then
            ivadj = icon(1,ivcur)
         else
            stop 3115
         endif
         if(ivadj.eq.ivini) go to 180
         if(ivadj.le.0) stop 3120
         ivcur = ivadj
         go to 20
c
c     boundary of cell intersects edge of rectangle
c     test whether it has a vertex in interior of rectangle
c
   40    continue
         ivini = ivcur
   50    continue
         if(icon(4,ivcur).eq.i) then
            ivadj = icon(2,ivcur)
         elseif(icon(5,ivcur).eq.i) then
            ivadj = icon(3,ivcur)
         elseif(icon(6,ivcur).eq.i) then
            ivadj = icon(1,ivcur)
         else
            stop 3125
         endif
         if(ivadj.eq.ivini) go to 180
         if(ivadj.le.0) stop 3130
         if(icon(4,ivadj).le.n .and. icon(5,ivadj).le.n .and.
     *      icon(6,ivadj).le.n) go to 60
         ivcur = ivadj
         go to 50
   60    continue
         is(i) = ivcur
  180 continue
c
c     mark Voronoi vertices outside rectangle
c
      do 300 i = 1, ivnxt
         if(icon(4,i).le.n .or. icon(5,i).le.n .or.
     *      icon(6,i).le.n) go to 300
         icon(4,i) = -1
  300 continue
c
c     for each Voronoi vertex not outside rectangle eliminate
c     Voronoi vertices outside rectangle that are its neighbors;
c     for corresponding Delaunay triangle eliminate vertices
c     of triangle that are reflection points
c
      do 600 i = 1, ivnxt
         if(icon(4,i).eq.-1) go to 600
         do 400 j = 1, 3
            ivcur = icon(j,i)
            if(ivcur.le.0) go to 400
            if(icon(4,ivcur).eq.-1) icon(j,i) = 0
  400    continue
         do 500 j = 4, 6
            if(icon(j,i).gt.n) icon(j,i) = 0
  500    continue
  600 continue
c
c     compress icon, vx, vy, ifl
c
      ivnot = 0
      do 1250 i = 1, ivnxt
          if(icon(4,i) .eq. -1) go to 1250
          ivnot = ivnot + 1
          id(i) = ivnot
          do 1200 j = 1, 6
              icon(j,ivnot) = icon(j,i)
 1200     continue
          vx(ivnot) = vx(i)
          vy(ivnot) = vy(i)
          ifl(ivnot) = ifl(i)
 1250 continue
      ivnxt = ivnot
c
c     update icon for triangles
c
      do 1400 i = 1, ivnxt
          do 1300 j = 1, 3
              ivcur = icon(j,i)
              if(ivcur .eq. 0) go to 1300
              icon(j,i) = id(ivcur)
 1300     continue
 1400 continue
c
c     update is for triangles
c
      do 1500 i = 1, n
          if(is(i) .le. 0) go to 1500
          is(i) = id(is(i))
 1500 continue
c
      return
      end
*CONTST
c
c     subroutine contst to -
c
c     test consistency of triangulation
c
c     November 9, 1993
c
      subroutine contst(icon, is, id, nv, ivnxt, nvmz)
c
      integer nv, ivnxt, nvmz
      integer icon(6,*), is(*), id(*)
      integer site1, site2, i, iscur, ivini, itrct
      integer ivcur, ivcur2, ivabs2, ivlst, isit1, isit2, itemp
c
c     test initial triangle for each site
c
      do 50 i = 1, nv
          iscur = is(i)
          if(iscur .le. 0) go to 50
          if(icon(4,iscur) .ne. i .and. icon(5,iscur) .ne. i .and.
     *       icon(6,iscur) .ne. i) stop 3200
   50 continue
c
c     initialize
c
      do 60 i = 1, ivnxt
          id(i) = 1
   60 continue
c
      do 70 i = 1, nv
          if(is(i) .gt. 0) go to 80
   70 continue
      stop 3220
   80 continue
      ivini = is(i)
      itrct = 1
      id(ivini) = 0
      ivcur = ivini
  100 continue
      ivcur2 = icon(2,ivcur)
      if(ivcur2 .eq. 0 .or. ivcur2 .eq. nvmz) go to 200
      ivabs2 = iabs(ivcur2)
      ivlst = ivcur
      if(ivcur2 .lt. 0) ivlst = -ivcur
      site1 = icon(6,ivcur)
      site2 = icon(4,ivcur)
      if(id(ivabs2) .eq. 1) go to 500
      if(icon(2,ivabs2) .eq. ivlst) then
          isit1 = icon(4,ivabs2)
          isit2 = icon(6,ivabs2)
      elseif(icon(3,ivabs2) .eq. ivlst) then
          isit1 = icon(5,ivabs2)
          isit2 = icon(4,ivabs2)
      elseif(icon(1,ivabs2) .eq. ivlst .and. ivabs2 .eq. ivini) then
          isit1 = icon(6,ivabs2)
          isit2 = icon(5,ivabs2)
      else
          stop 3240
      endif
      if(site1 .ne. isit1 .or. site2 .ne. isit2) stop 3260
  200 continue
      ivcur2 = icon(3,ivcur)
      if(ivcur2 .eq. 0 .or. ivcur2 .eq. nvmz) go to 300
      ivabs2 = iabs(ivcur2)
      ivlst = ivcur
      if(ivcur2 .lt. 0) ivlst = -ivcur
      site1 = icon(4,ivcur)
      site2 = icon(5,ivcur)
      if(id(ivabs2) .eq. 1) go to 500
      if(icon(2,ivabs2) .eq. ivlst) then
          isit1 = icon(4,ivabs2)
          isit2 = icon(6,ivabs2)
      elseif(icon(3,ivabs2) .eq. ivlst) then
          isit1 = icon(5,ivabs2)
          isit2 = icon(4,ivabs2)
      elseif(icon(1,ivabs2) .eq. ivlst .and. ivabs2 .eq. ivini) then
          isit1 = icon(6,ivabs2)
          isit2 = icon(5,ivabs2)
      else
          stop 3280
      endif
      if(site1 .ne. isit1 .or. site2 .ne. isit2) stop 3300
  300 continue
      ivcur2 = icon(1,ivcur)
      if(ivcur .eq. ivini) go to 400
      if(ivcur2 .eq. 0 .or. ivcur2 .eq. nvmz) stop 3320
      ivabs2 = iabs(ivcur2)
      ivlst = ivcur
      if(ivcur2 .lt. 0) ivlst = -ivcur
      if(icon(3,ivabs2) .eq. ivlst) then
          ivcur = ivabs2
          go to 300
      elseif(icon(2,ivabs2) .eq. ivlst) then
          ivcur = ivabs2
          go to 200
      elseif(icon(1,ivabs2) .eq. ivlst .and. ivabs2 .eq. ivini) then
          go to 900
      else
          stop 3340
      endif
  400 continue
      if(ivcur2 .eq. 0 .or. ivcur2 .eq. nvmz) go to 900
      ivabs2 = iabs(ivcur2)
      ivlst = ivcur
      if(ivcur2 .lt. 0) ivlst = -ivcur
      site1 = icon(5,ivcur)
      site2 = icon(6,ivcur)
      if(id(ivabs2) .eq. 1) go to 500
      if(icon(2,ivabs2) .eq. ivlst) then
          isit1 = icon(4,ivabs2)
          isit2 = icon(6,ivabs2)
      elseif(icon(3,ivabs2) .eq. ivlst) then
          isit1 = icon(5,ivabs2)
          isit2 = icon(4,ivabs2)
      else
          stop 3360
      endif
      if(site1 .ne. isit1 .or. site2 .ne. isit2) stop 3380
      go to 900
c
c     shift icon for ivabs2
c
  500 continue
      if(icon(1,ivabs2) .eq. ivlst) go to 600
      if(icon(2,ivabs2) .eq. ivlst) then
          itemp = icon(1,ivabs2)
          icon(1,ivabs2) = icon(2,ivabs2)
          icon(2,ivabs2) = icon(3,ivabs2)
          icon(3,ivabs2) = itemp
          itemp = icon(4,ivabs2)
          icon(4,ivabs2) = icon(5,ivabs2)
          icon(5,ivabs2) = icon(6,ivabs2)
          icon(6,ivabs2) = itemp
      elseif(icon(3,ivabs2) .eq. ivlst) then
          itemp = icon(1,ivabs2)
          icon(1,ivabs2) = icon(3,ivabs2)
          icon(3,ivabs2) = icon(2,ivabs2)
          icon(2,ivabs2) = itemp
          itemp = icon(4,ivabs2)
          icon(4,ivabs2) = icon(6,ivabs2)
          icon(6,ivabs2) = icon(5,ivabs2)
          icon(5,ivabs2) = itemp
      else
          stop 3390
      endif
c
      if(site1 .ne. icon(6,ivabs2) .or.
     *   site2 .ne. icon(5,ivabs2)) stop 3400
c
  600 continue
      ivcur = ivabs2
      id(ivcur) = 0
      itrct = itrct + 1
      go to 100
c
  900 continue
      if(itrct .ne. ivnxt) stop 3420
c
      do 1000 i = 1, ivnxt
          if(id(i) .ne. 0) stop 3440
 1000 continue
c
      write(*,*)'****************************************'
      write(*,*)'consistency check satisfied'
      write(*,*)'****************************************'
      return
      end
*CIRCRI
c
c     this subroutine tests how well the circle criterion is satisfied
c     by the triangles
c
      subroutine circri(x, y, ix, iy, ix2, iy2, icon, mhalf, mfull,
     *                  isclp, ivnxt, epf)
c
      integer mhalf, mfull, isclp(*)
      integer ivnxt
      double precision epf
      double precision x(*), y(*)
      integer ix(*), iy(*), ix2(*), iy2(*)
      integer icon(6,*)
      double precision dif1, dif2, dif3, delx1, dely1, delx2, dely2
      double precision delx3, dely3, tdist
      double precision vx, vy, zlamb, dat, dot, dot2
      integer isite, isadj, iscur, i, ishat, k, ivadj, j, insti
      integer icerr, itide
c
c     initialize
c
      icerr = 0
c
c     test all triangles
c
      do 200 i = 1, ivnxt
         isite = icon(4,i)
         isadj = icon(5,i)
         iscur = icon(6,i)
c
c    compute center for current triangle
c
         insti = 0
         delx2 = y(isadj) - y(iscur)
         dely2 = x(iscur) - x(isadj)
         dif2 = dsqrt(delx2**2 + dely2**2)
         if(dif2 .lt. epf) go to 30
c
         delx3 = x(isite) - x(isadj)
         dely3 = y(isite) - y(isadj)
         dif3 = dsqrt(delx3**2 + dely3**2)
         if(dif3 .lt. epf) go to 30
c
         delx1 = x(isite) - x(iscur)
         dely1 = y(isite) - y(iscur)
         dif1 = dsqrt(delx1**2 + dely1**2)
         if(dif1 .lt. epf) go to 30
c
         dot = delx2*delx1 + dely2*dely1
         if(dot.lt.epf .and. dot.gt.-epf) go to 30
         dot2 = dot/dif2
         if(dot2.lt.epf .and. dot2.gt.-epf) go to 30
         dat = delx1*delx3 + dely1*dely3
         zlamb = dat/dot
         vx = .5*(x(isadj)+x(iscur) + zlamb*delx2)
         vy = .5*(y(isadj)+y(iscur) + zlamb*dely2)
         go to 40
   30    continue
         insti = 1
   40    continue
c
         do 100 j = 1, 3
            ivadj = icon(j,i)
            if(ivadj.le.0) go to 100
            do 50 k = 1, 3
                 if(icon(k,ivadj).eq.i) go to 75
   50       continue
            stop 3560
   75       continue
            ishat = icon(k+3,ivadj)
            if(insti.eq.1) go to 90
            delx1 = x(ishat) - x(isite)
            dely1 = y(ishat) - y(isite)
            dif1 = dsqrt(delx1**2 + dely1**2)
            if(dif1 .lt. epf) go to 90
            call biscrc(x,y,ishat,isite,epf,vx,vy,tdist)
            if(tdist.le.-epf) go to 100
            if(tdist.lt.epf) go to 90
            icerr = icerr+1
            go to 100
   90       continue
            call circsn(ix, iy, ix2, iy2, isite, isadj, iscur, ishat,
     *                  mhalf, mfull, isclp, itide)
            if(itide .ge. 0) go to 100
            icerr = icerr+1
  100    continue
  200 continue
c
      write(*,*) '************************************************'
      write(*,*) 'Number of circle criterion violations = ',icerr
      write(*,*) '(This number should equal 0)'
      write(*,*) '************************************************'
      write(*,*) ' '
c
      return
      end
*ORIENT
c
c     This subroutine will test the orientation of the triangles
c
      subroutine orient(x, y, ix, iy, ix2, iy2, icon, nt,
     *                  mhalf, mfull, isclp, epf)
c
      double precision x(*), y(*)
      integer ix(*), iy(*), ix2(*), iy2(*)
      integer icon(6,*)
      double precision epf
      integer isclp(*), mhalf, mfull, nt
      integer a, b, c, i, isga, idmin
      double precision delx, dely, dif, xtemp, ytemp, dlst, dot
c
c     test all triangles
c
      idmin = 0
      do 300 i=1,nt
         a=icon(4,i)
         b=icon(5,i)
         c=icon(6,i)
         if(a.le.0 .or. b.le.0 .or. c.le.0) stop 3600
         delx = y(a) - y(b)
         dely = x(b) - x(a)
         dif = dsqrt(delx ** 2 + dely ** 2)
         if(dif .lt. epf) go to 220
         xtemp = x(c) - x(a)
         ytemp = y(c) - y(a)
         dlst=dsqrt(xtemp**2+ytemp**2)
         if(dlst .lt. epf) go to 220
         dot = (delx * xtemp + dely * ytemp)/dif
         if(dot.ge.epf .or. dot.le.-epf) go to 240
  220    continue
         call innsgn(ix, iy, ix2, iy2, a, b, c, mhalf, mfull, isclp,
     *               isga)
         if(isga.gt.0) then
            dot = 10.0d0
         elseif(isga.lt.0) then
            dot =-10.0d0
         else
            dot = 0.0d0
         endif
  240    continue
         if(dot .lt. epf) idmin = idmin+1
  300 continue
c
      write(*,*) '************************************************'
      write(*,*) 'Number of orientation violations = ',idmin
      write(*,*) '(This number should equal 0)'
      write(*,*) '************************************************'
      write(*,*) ' '
c
      return
      end
*CONVEX
c
c  Subroutine to verify convexity of union of triangles
c
      subroutine convex(x, y, ix, iy, ix2, iy2, icon, is, nv, nt,
     *                  nvmz, mhalf, mfull, isclp, epf)
c
      double precision x(*), y(*)
      integer ix(*), iy(*), ix2(*), iy2(*)
      integer icon(6,*), is(*), nv, nt, nvmz
      double precision epf
      integer isclp(*), mhalf, mfull, isga
      integer idmin, i, ivcur, isite, iscur, isadj, isini, ivadj, itot
      double precision delx, dely, dif, xtemp, ytemp, dlst, dot
c
      idmin = 0
c
c     identify a triangle at boundary of convex hull
c
      do 300 i = 1, nt
         if(icon(1,i).eq.0 .or. icon(1,i).eq.nvmz .or.
     *      icon(2,i).eq.0 .or. icon(2,i).eq.nvmz .or.
     *      icon(3,i).eq.0 .or. icon(3,i).eq.nvmz) go to 400
  300 continue
      stop 3700
  400 continue
      ivcur = i
c
c     identify first two consecutive vertices at boundary of convex hull
c
      if(icon(1,ivcur).eq.0 .or. icon(1,ivcur).eq.nvmz) then
         isite = icon(6,ivcur)
         iscur = icon(5,ivcur)
      elseif(icon(2,ivcur).eq.0 .or. icon(2,ivcur).eq.nvmz) then
         isite = icon(4,ivcur)
         iscur = icon(6,ivcur)
      elseif(icon(3,ivcur).eq.0 .or. icon(3,ivcur).eq.nvmz) then
         isite = icon(5,ivcur)
         iscur = icon(4,ivcur)
      else
         stop 3720
      endif
      isini = isite
c
c     test current vertex isite at boundary of convex hull
c
      itot = 0
  500 continue
      itot = itot+1
      if(itot.gt.nv) stop 3730
      ivcur = is(isite)
      if(ivcur.le.0 .or. ivcur.gt.nt) stop 3740
      if(icon(4,ivcur).eq.isite) then
         ivadj = icon(3,ivcur)
         isadj = icon(5,ivcur)
      elseif(icon(5,ivcur).eq.isite) then
         ivadj = icon(1,ivcur)
         isadj = icon(6,ivcur)
      elseif(icon(6,ivcur).eq.isite) then
         ivadj = icon(2,ivcur)
         isadj = icon(4,ivcur)
      else
         stop 3760
      endif
      if(ivadj.ne.0 .and. ivadj.ne.nvmz) stop 3780
c
      delx = y(isite) - y(iscur)
      dely = x(iscur) - x(isite)
      dif = dsqrt(delx ** 2 + dely ** 2)
      if(dif .lt. epf) go to 520
      xtemp = x(isadj) - x(isite)
      ytemp = y(isadj) - y(isite)
      dlst=dsqrt(xtemp**2+ytemp**2)
      if(dlst .lt. epf) go to 520
      dot = (delx * xtemp + dely * ytemp)/dif
      if(dot.ge.epf .or. dot.le.-epf) go to 540
  520 continue
      call innsgn(ix, iy, ix2, iy2, isite, iscur, isadj, mhalf, mfull,
     *            isclp, isga)
      if(isga.gt.0) then
         dot = 10.0d0
      elseif(isga.lt.0) then
         dot =-10.0d0
      else
         dot = 0.0d0
      endif
  540 continue
      if(dot .ge. epf) idmin = idmin+1
      iscur = isite
      isite = isadj
      if(isite.ne.isini) go to 500
c
      write(*,*) '************************************************'
      write(*,*) 'Number of convexity violations = ',idmin
      write(*,*) '(This number should equal 0)'
      write(*,*) '************************************************'
      write(*,*) ' '
c
      return
      end
*INNSGN
c
c     Routine for determining sign of distance from point ifou
c     to line that contains point ifir and isec. Distance is
c     positive if points ifir, isec, ifou are the vertices of
c     a triangle with nonempty interior and appear in this order
c     in a counterclockwise direction around the boundary of the
c     triangle.
c
      subroutine innsgn(x, y, x2, y2, ifir, isec, ifou,
     *                  mhalf, mfull, isclp, isga)
c
      integer x(*), y(*), x2(*), y2(*)
      integer ifir, isec, ifou
      integer isclp(*), mhalf, mfull, nkmax
      parameter (nkmax = 30)
      integer ia(nkmax), io(nkmax), iax(nkmax), iay(nkmax)
      integer iu(nkmax), iw(nkmax)
      integer ix4(nkmax), iy4(nkmax)
      integer ixf(nkmax), iyf(nkmax)
      integer ixfiw, iyfiw, ixsew, iysew
      integer ixfow, iyfow
      integer ixfi2, iyfi2, ixse2, iyse2
      integer ixfo2, iyfo2
      integer isgxf, isgyf, ikxf, ikyf
      integer isgx4, isgy4, ikx4, iky4
      integer isga, ika, isgo, iko
      integer isgax, ikax, isgay, ikay
      integer isgu, isgw, iku, ikw
c
      ixfiw = x(ifir)
      iyfiw = y(ifir)
      ixsew = x(isec)
      iysew = y(isec)
      ixfow = x(ifou)
      iyfow = y(ifou)
c
      ixfi2 = x2(ifir)
      iyfi2 = y2(ifir)
      ixse2 = x2(isec)
      iyse2 = y2(isec)
      ixfo2 = x2(ifou)
      iyfo2 = y2(ifou)
c
      call decmp2(ixf, isgxf, ikxf, ixfiw, ixfi2, mhalf, mfull, isclp)
      call decmp2(iyf, isgyf, ikyf, iyfiw, iyfi2, mhalf, mfull, isclp)
c
      call decmp2(io, isgo, iko, ixsew, ixse2, mhalf, mfull, isclp)
      call muldif(io, ixf, iay, isgo, isgxf, isgay, iko, ikxf, ikay,
     *              nkmax, mhalf)
      call decmp2(iu, isgu, iku, ixfow, ixfo2, mhalf, mfull, isclp)
      call muldif(iu, ixf, ix4, isgu, isgxf, isgx4, iku, ikxf, ikx4,
     *              nkmax, mhalf)
      call decmp2(io, isgo, iko, iysew, iyse2, mhalf, mfull, isclp)
      call muldif(io, iyf, iax, isgo, isgyf, isgax, iko, ikyf, ikax,
     *              nkmax, mhalf)
      isgax=-isgax
      call decmp2(iu, isgu, iku, iyfow, iyfo2, mhalf, mfull, isclp)
      call muldif(iu, iyf, iy4, isgu, isgyf, isgy4, iku, ikyf, iky4,
     *              nkmax, mhalf)
c
      call mulmul(iax, ix4, iw, isgax, isgx4, isgw, ikax, ikx4, ikw,
     *              nkmax, mhalf)
c
      call mulmul(iay, iy4, iu, isgay, isgy4, isgu, ikay, iky4, iku,
     *              nkmax, mhalf)
      isgu =-isgu
      call muldif(iw, iu, ia, isgw, isgu, isga, ikw, iku, ika,
     *              nkmax, mhalf)
c
      return
      end
*DISTSN
c
c     Routine for determining sign of distance from point ifou
c     to line that contains point ifir and that is perpendicular
c     to vector (iax,iay) with origin point ifir. The distance
c     is positive if point ifou is in halfplane that contains
c     vector (iax,iay).
c
      subroutine distsn(x, y, x2, y2, ifir, ifou, mhalf, mfull, isclp,
     *                  iax, isgax, ikax, iay, isgay, ikay, isga)
c
      integer x(*), y(*), x2(*), y2(*)
      integer ifir, ifou
      integer isclp(*), mhalf, mfull, nkmax
      parameter (nkmax = 30)
      integer iax(*), iay(*)
      integer ia(nkmax), iu(nkmax), iw(nkmax)
      integer ix4(nkmax), iy4(nkmax)
      integer ixf(nkmax), iyf(nkmax)
      integer ixfiw, iyfiw, ixfow, iyfow
      integer ixfi2, iyfi2, ixfo2, iyfo2
      integer isgxf, isgyf, ikxf, ikyf
      integer isgx4, isgy4, ikx4, iky4
      integer isga, ika
      integer isgax, ikax, isgay, ikay
      integer isgu, isgw, iku, ikw
c
      ixfiw = x(ifir)
      iyfiw = y(ifir)
      ixfow = x(ifou)
      iyfow = y(ifou)
c
      ixfi2 = x2(ifir)
      iyfi2 = y2(ifir)
      ixfo2 = x2(ifou)
      iyfo2 = y2(ifou)
c
      call decmp2(ixf, isgxf, ikxf, ixfiw, ixfi2, mhalf, mfull, isclp)
      call decmp2(iyf, isgyf, ikyf, iyfiw, iyfi2, mhalf, mfull, isclp)
c
      call decmp2(iu, isgu, iku, ixfow, ixfo2, mhalf, mfull, isclp)
      call muldif(iu, ixf, ix4, isgu, isgxf, isgx4, iku, ikxf, ikx4,
     *              nkmax, mhalf)
      call decmp2(iu, isgu, iku, iyfow, iyfo2, mhalf, mfull, isclp)
      call muldif(iu, iyf, iy4, isgu, isgyf, isgy4, iku, ikyf, iky4,
     *              nkmax, mhalf)
c
      call mulmul(iax, ix4, iw, isgax, isgx4, isgw, ikax, ikx4, ikw,
     *              nkmax, mhalf)
c
      call mulmul(iay, iy4, iu, isgay, isgy4, isgu, ikay, iky4, iku,
     *              nkmax, mhalf)
      isgu =-isgu
      call muldif(iw, iu, ia, isgw, isgu, isga, ikw, iku, ika,
     *              nkmax, mhalf)
c
      return
      end
*VECTRD
c
c     Routine for determining vector <ifir,isec>
c
      subroutine vectrd(x, y, x2, y2, ifir, isec, mhalf, mfull, isclp,
     *                  iax, isgax, ikax, iay, isgay, ikay)
c
      integer x(*), y(*), x2(*), y2(*)
      integer iax(*),iay(*)
      integer ifir, isec
      integer isclp(*), mhalf, mfull, nkmax
      parameter (nkmax = 30)
      integer io(nkmax)
      integer ixf(nkmax), iyf(nkmax)
      integer ixfiw, iyfiw, ixsew, iysew
      integer ixfi2, iyfi2, ixse2, iyse2
      integer isgxf, isgyf, ikxf, ikyf
      integer isgax, ikax, isgay, ikay
      integer isgo, iko
c
      ixfiw = x(ifir)
      iyfiw = y(ifir)
      ixsew = x(isec)
      iysew = y(isec)
c
      ixfi2 = x2(ifir)
      iyfi2 = y2(ifir)
      ixse2 = x2(isec)
      iyse2 = y2(isec)
c
      call decmp2(ixf, isgxf, ikxf, ixfiw, ixfi2, mhalf, mfull, isclp)
      call decmp2(iyf, isgyf, ikyf, iyfiw, iyfi2, mhalf, mfull, isclp)
c
      call decmp2(io, isgo, iko, ixsew, ixse2, mhalf, mfull, isclp)
      call muldif(io, ixf, iax, isgo, isgxf, isgax, iko, ikxf, ikax,
     *              nkmax, mhalf)
      call decmp2(io, isgo, iko, iysew, iyse2, mhalf, mfull, isclp)
      call muldif(io, iyf, iay, isgo, isgyf, isgay, iko, ikyf, ikay,
     *              nkmax, mhalf)
c
      return
      end
*PRPVEC
c
c     Routine for determining vector perpendicular to <ifir,isec>
c     that makes a right angle to <ifir,isec> in a counterclockwise
c     direction from it
c
      subroutine prpvec(x, y, x2, y2, ifir, isec, mhalf, mfull, isclp,
     *                  iax, isgax, ikax, iay, isgay, ikay)
c
      integer x(*), y(*), x2(*), y2(*)
      integer iax(*),iay(*)
      integer ifir, isec
      integer isclp(*), mhalf, mfull, nkmax
      parameter (nkmax = 30)
      integer io(nkmax)
      integer ixf(nkmax), iyf(nkmax)
      integer ixfiw, iyfiw, ixsew, iysew
      integer ixfi2, iyfi2, ixse2, iyse2
      integer isgxf, isgyf, ikxf, ikyf
      integer isgax, ikax, isgay, ikay
      integer isgo, iko
c
      ixfiw = x(ifir)
      iyfiw = y(ifir)
      ixsew = x(isec)
      iysew = y(isec)
c
      ixfi2 = x2(ifir)
      iyfi2 = y2(ifir)
      ixse2 = x2(isec)
      iyse2 = y2(isec)
c
      call decmp2(ixf, isgxf, ikxf, ixfiw, ixfi2, mhalf, mfull, isclp)
      call decmp2(iyf, isgyf, ikyf, iyfiw, iyfi2, mhalf, mfull, isclp)
c
      call decmp2(io, isgo, iko, ixsew, ixse2, mhalf, mfull, isclp)
      call muldif(io, ixf, iay, isgo, isgxf, isgay, iko, ikxf, ikay,
     *              nkmax, mhalf)
      call decmp2(io, isgo, iko, iysew, iyse2, mhalf, mfull, isclp)
      call muldif(io, iyf, iax, isgo, isgyf, isgax, iko, ikyf, ikax,
     *              nkmax, mhalf)
      isgax=-isgax
c
      return
      end
*CIRCSN
c
c     subroutine for determining position of point ifou with
c     respect to circle determined by points ifir, isec, ithi
c     if positive then ifou is outside the circle
c     if negative then ifou is inside the circle
c     if zero then ifou is on the boundary of the circle
c
      subroutine circsn(x, y, x2, y2, ifir, isec, ithi, ifou,
     *                  mhalf, mfull, isclp, ipout)
c
      integer x(*), y(*), x2(*), y2(*)
      integer ifir, isec, ithi, ifou, mhalf, mfull, nkmax, ipout
      integer isclp(*)
      parameter (nkmax = 30)
      integer iu(nkmax)
      integer iq2(nkmax), iq3(nkmax), iq4(nkmax)
      integer ix2(nkmax), iy2(nkmax)
      integer ix3(nkmax), iy3(nkmax)
      integer ix4(nkmax), iy4(nkmax)
      integer ixf(nkmax), iyf(nkmax)
      integer ixf2(nkmax), iyf2(nkmax)
      integer ixthw, iythw, ixfow, iyfow
      integer ixsew, iysew, ixfuw, iyfuw
      integer ixth2, iyth2, ixfo2, iyfo2
      integer ixse2, iyse2, ixfu2, iyfu2
      integer isgq2, isgq3, isgq4, ikq2, ikq3, ikq4
      integer isgxf, isgyf, ikxf, ikyf
      integer isgxf2, isgyf2, ikxf2, ikyf2
      integer isgx2, isgy2, ikx2, iky2
      integer isgx3, isgy3, ikx3, iky3 
      integer isgx4, isgy4, ikx4, iky4
      integer isgu, iku
c
      ixfuw = x(ifir)
      iyfuw = y(ifir)
      ixsew = x(isec)
      iysew = y(isec)
      ixthw = x(ithi)
      iythw = y(ithi)
      ixfow = x(ifou)
      iyfow = y(ifou)
c
      ixfu2 = x2(ifir)
      iyfu2 = y2(ifir)
      ixse2 = x2(isec)
      iyse2 = y2(isec)
      ixth2 = x2(ithi)
      iyth2 = y2(ithi)
      ixfo2 = x2(ifou)
      iyfo2 = y2(ifou)
c
      call decmp2(ixf, isgxf, ikxf, ixfuw, ixfu2, mhalf, mfull, isclp)
      call decmp2(iyf, isgyf, ikyf, iyfuw, iyfu2, mhalf, mfull, isclp)
      call mulmul(ixf, ixf, ixf2, isgxf, isgxf, isgxf2, ikxf, ikxf,
     *              ikxf2, nkmax, mhalf)
      call mulmul(iyf, iyf, iyf2, isgyf, isgyf, isgyf2, ikyf, ikyf,
     *              ikyf2, nkmax, mhalf)
      if(isgxf2.lt.0 .or. isgyf2.lt.0) stop 4600
c
      call frterm(ixsew, iysew, ixf, iyf, isgxf, isgyf, ikxf, ikyf,
     *            ixf2, iyf2, isgxf2, isgyf2, ikxf2, ikyf2, ix2, iy2,
     *            iq2, isgx2, isgy2, isgq2, ikx2, iky2, ikq2,
     *            mhalf, mfull, ixse2, iyse2, isclp)
c
      call frterm(ixthw, iythw, ixf, iyf, isgxf, isgyf, ikxf, ikyf,
     *            ixf2, iyf2, isgxf2, isgyf2, ikxf2, ikyf2, ix3, iy3,
     *            iq3, isgx3, isgy3, isgq3, ikx3, iky3, ikq3,
     *            mhalf, mfull, ixth2, iyth2, isclp)
c
      call frterm(ixfow, iyfow, ixf, iyf, isgxf, isgyf, ikxf, ikyf,
     *            ixf2, iyf2, isgxf2, isgyf2, ikxf2, ikyf2, ix4, iy4,
     *            iq4, isgx4, isgy4, isgq4, ikx4, iky4, ikq4,
     *            mhalf, mfull, ixfo2, iyfo2, isclp)
c
      call detrm3(iq2, ix2, iy2, isgq2, isgx2, isgy2,
     *            iq3, ix3, iy3, isgq3, isgx3, isgy3,
     *            iq4, ix4, iy4, isgq4, isgx4, isgy4,
     *            ikq2, ikq3, ikq4, ikx2, ikx3, ikx4,
     *            iky2, iky3, iky4, iu, isgu, iku, mhalf)
c
      ipout = isgu
c
      return
      end
*NUMDEN
c
c     subroutine for determining numerator and denominator used in
c     formulas for computing x-,y-coordinates of center of circle
c     determined by points ifir, isec, ithi
c
      subroutine numden(x, y, x2, y2, ifir, isec, ithi, mhalf, mfull,
     *                  isclp, io, isgo, iko, inx, iny, isgnx, isgny,
     *                  iknx, ikny)
c
      integer x(*), y(*), x2(*), y2(*)
      integer io(*), isgo, iko
      integer inx(*), iny(*), isgnx, isgny, iknx, ikny
      integer ifir, isec, ithi, mhalf, mfull, nkmax
      integer isclp(*)
      parameter (nkmax = 30)
      integer iq2(nkmax), iq3(nkmax)
      integer ix2(nkmax), iy2(nkmax)
      integer ix3(nkmax), iy3(nkmax)
      integer ixf(nkmax), iyf(nkmax)
      integer ixf2(nkmax), iyf2(nkmax)
      integer ixthw, iythw, ixsew, iysew, ixfuw, iyfuw
      integer ixth2, iyth2, ixse2, iyse2, ixfu2, iyfu2
      integer isgq2, isgq3, ikq2, ikq3
      integer isgxf, isgyf, ikxf, ikyf
      integer isgxf2, isgyf2, ikxf2, ikyf2
      integer isgx2, isgy2, ikx2, iky2
      integer isgx3, isgy3, ikx3, iky3 
c
      ixfuw = x(ifir)
      iyfuw = y(ifir)
      ixsew = x(isec)
      iysew = y(isec)
      ixthw = x(ithi)
      iythw = y(ithi)
c
      ixfu2 = x2(ifir)
      iyfu2 = y2(ifir)
      ixse2 = x2(isec)
      iyse2 = y2(isec)
      ixth2 = x2(ithi)
      iyth2 = y2(ithi)
c
      call decmp2(ixf, isgxf, ikxf, ixfuw, ixfu2, mhalf, mfull, isclp)
      call decmp2(iyf, isgyf, ikyf, iyfuw, iyfu2, mhalf, mfull, isclp)
      call mulmul(ixf, ixf, ixf2, isgxf, isgxf, isgxf2, ikxf, ikxf,
     *              ikxf2, nkmax, mhalf)
      call mulmul(iyf, iyf, iyf2, isgyf, isgyf, isgyf2, ikyf, ikyf,
     *              ikyf2, nkmax, mhalf)
      if(isgxf2.lt.0 .or. isgyf2.lt.0) stop 4610
c
      call frterm(ixsew, iysew, ixf, iyf, isgxf, isgyf, ikxf, ikyf,
     *            ixf2, iyf2, isgxf2, isgyf2, ikxf2, ikyf2, ix2, iy2,
     *            iq2, isgx2, isgy2, isgq2, ikx2, iky2, ikq2,
     *            mhalf, mfull, ixse2, iyse2, isclp)
c
      call frterm(ixthw, iythw, ixf, iyf, isgxf, isgyf, ikxf, ikyf,
     *            ixf2, iyf2, isgxf2, isgyf2, ikxf2, ikyf2, ix3, iy3,
     *            iq3, isgx3, isgy3, isgq3, ikx3, iky3, ikq3,
     *            mhalf, mfull, ixth2, iyth2, isclp)
c
      call dterm2(ix2, iy2, isgx2, isgy2, ix3, iy3, isgx3, isgy3,
     *            ikx2, ikx3, iky2, iky3, io, isgo, iko, mhalf)
c
      call dterm2(iq2, iy2, isgq2, isgy2, iq3, iy3, isgq3, isgy3,
     *            ikq2, ikq3, iky2, iky3, inx, isgnx, iknx, mhalf)
c
      call dterm2(ix2, iq2, isgx2, isgq2, ix3, iq3, isgx3, isgq3,
     *            ikx2, ikx3, ikq2, ikq3, iny, isgny, ikny, mhalf)
c
      return
      end
*FRTERM
c
      subroutine frterm(ixsew, iysew, ixf, iyf, isgxf, isgyf, ikxf,
     *                  ikyf, ixf2, iyf2, isgxf2, isgyf2, ikxf2, ikyf2,
     *                  ix2, iy2, iq2, isgx2, isgy2, isgq2, ikx2, iky2,
     *                  ikq2, mhalf, mfull, ixse2, iyse2, isclp)
c
      integer ixsew, iysew, isgxf, isgyf, isgxf2, isgyf2
      integer ikxf, ikyf, ikxf2, ikyf2
      integer isgx2, isgy2, isgq2
      integer ikx2, iky2, ikq2, mhalf, mfull
      integer ixse2, iyse2, isclp(*)
      integer isgo, isgu, isgv, isgp, iko, iku, ikv, ikp
      integer ixf(*), iyf(*), ixf2(*), iyf2(*)
      integer ix2(*), iy2(*), iq2(*)
      integer nkmax
      parameter (nkmax = 30)
      integer io(nkmax), iu(nkmax), iv(nkmax), ip(nkmax)
c
      call decmp2(io, isgo, iko, ixsew, ixse2, mhalf, mfull, isclp)
      call muldif(io, ixf, ix2, isgo, isgxf, isgx2, iko, ikxf, ikx2,
     *              nkmax, mhalf)
      call mulmul(io, io, iu, isgo, isgo, isgu, iko, iko, iku,
     *              nkmax, mhalf)
      call muldif(iu, ixf2, ip, isgu, isgxf2, isgp, iku, ikxf2, ikp,
     *              nkmax, mhalf)
c
      call decmp2(io, isgo, iko, iysew, iyse2, mhalf, mfull, isclp)
      call muldif(io, iyf, iy2, isgo, isgyf, isgy2, iko, ikyf, iky2,
     *              nkmax, mhalf)
      call mulmul(io, io, iu, isgo, isgo, isgu, iko, iko, iku,
     *              nkmax, mhalf)
      call muldif(iu, iyf2, iv, isgu, isgyf2, isgv, iku, ikyf2, ikv,
     *              nkmax, mhalf)
      isgv=-isgv
      call muldif(ip, iv, iq2, isgp, isgv, isgq2, ikp, ikv, ikq2,
     *              nkmax, mhalf)
c
      return
      end
*DETRM3
c
      subroutine detrm3(ix2, iy2, iz2, isgx2, isgy2, isgz2,
     *                  ix3, iy3, iz3, isgx3, isgy3, isgz3,
     *                  ix4, iy4, iz4, isgx4, isgy4, isgz4,
     *                  ikx2, ikx3, ikx4, iky2, iky3, iky4,
     *                  ikz2, ikz3, ikz4, io, isgo, iko, mhalf)
c
      integer nkmax
      parameter (nkmax = 30)
      integer io(*), iu(nkmax), iv(nkmax), iw(nkmax)
      integer ix2(*), iy2(*), iz2(*), ix3(*), iy3(*), iz3(*)
      integer ix4(*), iy4(*), iz4(*)
      integer isgx2, isgy2, isgz2, isgx3, isgy3, isgz3
      integer isgx4, isgy4, isgz4, ikx2, ikx3, ikx4, iky2, iky3, iky4
      integer ikz2, ikz3, ikz4, isgo, iko, mhalf
      integer isgu, isgv, isgw, iku, ikv, ikw
c
      call mulmul(ix3, iy4, iv, isgx3, isgy4, isgv, ikx3, iky4, ikv,
     *              nkmax, mhalf)
      call mulmul(ix4, iy3, iu, isgx4, isgy3, isgu, ikx4, iky3, iku,
     *              nkmax, mhalf)
      call muldif(iv, iu, iw, isgv, isgu, isgw, ikv, iku, ikw,
     *              nkmax, mhalf)
      call mulmul(iw, iz2, io, isgw, isgz2, isgo, ikw, ikz2, iko,
     *              nkmax, mhalf)
c
      call mulmul(ix2, iy4, iv, isgx2, isgy4, isgv, ikx2, iky4, ikv,
     *              nkmax, mhalf)
      call mulmul(ix4, iy2, iu, isgx4, isgy2, isgu, ikx4, iky2, iku,
     *              nkmax, mhalf)
      call muldif(iv, iu, iw, isgv, isgu, isgw, ikv, iku, ikw,
     *              nkmax, mhalf)
      call mulmul(iw, iz3, iu, isgw, isgz3, isgu, ikw, ikz3, iku,
     *              nkmax, mhalf)
      call muldif(io, iu, iw, isgo, isgu, isgw, iko, iku, ikw,
     *              nkmax, mhalf)
c
      call mulmul(ix3, iy2, iv, isgx3, isgy2, isgv, ikx3, iky2, ikv,
     *              nkmax, mhalf)
      call mulmul(ix2, iy3, iu, isgx2, isgy3, isgu, ikx2, iky3, iku,
     *              nkmax, mhalf)
      call muldif(iv, iu, io, isgv, isgu, isgo, ikv, iku, iko,
     *              nkmax, mhalf)
      call mulmul(io, iz4, iu, isgo, isgz4, isgu, iko, ikz4, iku,
     *              nkmax, mhalf)
      call muldif(iw, iu, io, isgw, isgu, isgo, ikw, iku, iko,
     *              nkmax, mhalf)
c
      return
      end
*DTERM2
c
      subroutine dterm2(ix2, iy2, isgx2, isgy2, ix3, iy3, isgx3, isgy3,
     *                  ikx2, ikx3, iky2, iky3, io, isgo, iko, mhalf)
c
      integer nkmax
      parameter (nkmax = 30)
      integer io(*), iu(nkmax), iv(nkmax)
      integer ix2(*), iy2(*), ix3(*), iy3(*)
      integer isgx2, isgy2, isgx3, isgy3, ikx2, ikx3, iky2, iky3
      integer isgo, iko, mhalf
      integer isgu, isgv, iku, ikv
c
      call mulmul(ix2, iy3, iu, isgx2, isgy3, isgu, ikx2, iky3, iku,
     *              nkmax, mhalf)
      call mulmul(ix3, iy2, iv, isgx3, isgy2, isgv, ikx3, iky2, ikv,
     *              nkmax, mhalf)
      call muldif(iu, iv, io, isgu, isgv, isgo, iku, ikv, iko,
     *              nkmax, mhalf)
c
      return
      end
*DECMP2
c
c     subroutine decmp2 
c
c     to decompose a regular or non-regular length integer
c
      subroutine decmp2(ia, isga, ika, iwi, iwi2, mhalf, mfull, isclp)
c
      integer ia(*), isga, ika, iwi, iwi2, mhalf, mfull, isclp(*)
      integer nkmax
      parameter (nkmax = 30)
      integer iu(nkmax), io(nkmax), isgu, isgo, iku, iko, isgcl, ikcl
c
      call decomp(ia, isga, iwi, mhalf)
      ika = 2
      if(iwi2.ne.0) then
         isgcl = 1
         ikcl = 2
         call mulmul(ia, isclp, iu, isga, isgcl, isgu, ika, ikcl,
     *                 iku, nkmax, mhalf)
         if(iwi2.eq.mfull) iwi2 = 0
         call decomp(io, isgo, iwi2, mhalf)
         isgo = -isgo
         iko = 2
         call muldif(iu, io, ia, isgu, isgo, isga, iku, iko, ika,
     *                 nkmax, mhalf)
      endif
c
      return
      end
*DECOMP
c
c     subroutine decomp
c
c     to decompose a regular length integer
c
c     iwi = isga*(ia(1) + ia(2) * mhalf)
c
c     iwi  is a regular length integer
c     isga is a sign integer (-1, 0, 1)
c     ia(1) and ia(2) are integers less than mhalf
c
      subroutine decomp(ia, isga, iwi, mhalf)
c
      integer ia(*), isga, iwi, mhalf, ivi
c
      if(iwi.gt.0) then
         isga = 1
         ivi = iwi
      elseif(iwi.lt.0) then
         isga =-1
         ivi = -iwi
      else
         isga = 0
         ia(1) = 0
         ia(2) = 0
         return
      endif
      ia(2) = ivi/mhalf
      ia(1) = ivi - ia(2)*mhalf
c
      return
      end
*MULMUL
c
c     subroutine mulmul
c
c     to perform a multiple precision integer multiplication
c     (for multiplying 2 or more integers)
c
c     io = ia * ib
c
c     ia represents a decomposed integer
c     ib represents a decomposed integer
c     io is the product of ia and ib in its decomposed form
c
      subroutine mulmul(ia, ib, io, isga, isgb, isgo, ika, ikb, iko,
     *                    nkmax, mhalf)
c
      integer ia(*), ib(*), io(*)
      integer isga, isgb, isgo, ika, ikb, iko, nkmax, mhalf
      integer i, ipt, ipr, iko1, k, j
c
      if(isga.eq.0.or.isgb.eq.0)then
         isgo=0
         iko = 2
         io(1) = 0
         io(2) = 0
         return
      endif
c
      iko = ika + ikb
      if(iko.gt.nkmax) stop 4710
c
      if(isga.gt.0)then
         if(isgb.gt.0)then
            isgo = 1
         else
            isgo =-1
         endif
      else
         if(isgb.gt.0)then
            isgo =-1
         else
            isgo = 1
         endif
      endif
c
      iko1 = iko - 1
      ipr = 0
c
      do 200 i = 1, iko1
         ipt = ipr
         k = i
         do 180 j = 1, ikb
            if(k .lt. 1) go to 190
            if(k .gt. ika) go to 150
            ipt = ipt + ia(k)*ib(j)
  150       continue
            k = k - 1
  180    continue
  190    continue
         ipr = ipt/mhalf
         io(i) = ipt - ipr*mhalf
  200 continue
c
      io(iko) = ipr
      if(ipr.ge.mhalf) stop 4720
c
      iko1 = iko
      do 300 i = iko1, ika+1, -1
         if(io(i) .ne. 0) go to 400
         iko = iko - 1
  300 continue
  400 continue
c
      return
      end
*MULDIF
c
c     subroutine muldif
c
c     to perform a multiple precision integer subtraction
c     (for subtracting a decomposed product from another)
c
c     io = ia - ib
c
c     ia represents a decomposed regular length integer or the 
c        decomposed product of two or more regular length integers
c     ib is similarly described
c     io is a decomposed integer which represents ia - ib
c
      subroutine muldif(ia, ib, io, isga, isgb, isgo, ika, ikb, iko,
     *                    nkmax, mhalf)
c
      integer ia(*), ib(*), io(*)
      integer isga, isgb, isgo, ika, ikb, iko, nkmax, mhalf
      integer i, iko1, irel
c
      if(isgb.eq.0)then
         if(isga.eq.0)then
            isgo=0
            iko = 2
            io(1) = 0
            io(2) = 0
            return
         endif
         isgo = isga
         iko = ika
         do 100 i=1,iko
            io(i) = ia(i)
  100    continue
      elseif(isga.eq.0)then
         isgo =-isgb
         iko = ikb
         do 200 i=1,iko
            io(i) = ib(i)
  200    continue
      else
         iko = ika
         if(ikb.lt.ika) then
            do 300 i=ikb+1,ika
               ib(i) = 0
  300       continue
         elseif(ika.lt.ikb) then
            iko = ikb
            do 400 i=ika+1,ikb
               ia(i) = 0
  400       continue
         endif
         if(isga*isgb.gt.0)then
            irel = 0
            do 500 i = iko, 1, -1
               if(ia(i).gt.ib(i))then
                  irel = 1
                  go to 600
               elseif(ia(i).lt.ib(i))then
                  irel = -1
                  go to 600
               endif
  500       continue
  600       continue
            if(irel.eq.0)then
               isgo = 0
               do 700 i=1,iko
                  io(i) = 0
  700          continue
            else
               isgo=isga*irel
               io(1) = (ia(1)-ib(1))*irel
               do 800 i=2,iko
                  if(io(i-1).lt.0) then
                     io(i) =-1
                     io(i-1) = io(i-1) + mhalf
                  else
                     io(i) = 0
                  endif
                  io(i) = io(i) + (ia(i)-ib(i))*irel
  800          continue
               if(io(iko).lt.0) stop 4810
            endif
         else
            isgo=isga
            io(1) = ia(1)+ib(1)
            do 900 i=2,iko
               if(io(i-1).ge.mhalf) then
                  io(i) = 1
                  io(i-1) = io(i-1) - mhalf
               else
                  io(i) = 0
               endif
               io(i) = io(i) + ia(i)+ib(i)
  900       continue
            if(io(iko).ge.mhalf) then
               iko = iko+1
               if(iko.gt.nkmax) stop 4820
               io(iko) = 1
               io(iko-1) = io(iko-1) - mhalf
            endif
         endif
      endif
c
      if(iko .eq. 2) go to 1400
      iko1 = iko
      do 1300 i = iko1, 3, -1
         if(io(i) .ne. 0) go to 1400
         iko = iko - 1
 1300 continue
 1400 continue
c
      return
      end
*DOUBNM
c
      subroutine doubnm(io, isgo, iko, r215, dnum)
c
      integer io(*)
      double precision dnum, r215, rpwr
      integer isgo, iko, i
c
      if(isgo.eq.0) then
         dnum = dble(0)
         go to 1000
      else
         if(iko .lt. 2) stop 5710
         if(iko .gt. 68) stop 5720
         rpwr = dble(1)
         dnum = dble(io(1))
         do 100 i = 2, iko
            rpwr = rpwr*r215
            dnum = dnum + dble(io(i))*rpwr
  100    continue
      endif
      if(isgo.lt.0) dnum = -dnum
c
 1000 continue
      return
      end
