c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12
*MAIN
c
c     main program to -
c
c     compute Delaunay triangulations & Voronoi diagrams of dynamic
c     data, i. e. of a set of moving points in 2-d space.
c
c     Currently program is set up to handle 100 or fewer input files
c     making up dynamic data. Each input data set contains the same
c     number of distinct data points in the plane. Each time an input
c     data set is read the current Delaunay triangulation and Voronoi
c     diagram are updated on a local basis, i. e. taking into account
c     only those points that have moved. If the percentage of points
c     that have moved is too high then the whole Delaunay triangulation
c     and Voronoi diagram are computed.
c
c     The input files with data points have names pn1, pn2, ...., pn100.
c     The output files with triangulation information have names
c     tn1, tn2, ..., tn100.
c     The output files with information about Voronoi vertices have
c     names vn1, vn2, ..., vn100.
c
c     Unless the whole Delaunay triangulation and Voronoi diagram
c     are computed, triangulation information and information about
c     Voronoi vertices are saved in the ouput files on a local basis,
c     i. e. triangulation information and information about Voronoi
c     vertices that have not been affected by the moving points are
c     not saved.
c
c     After an input data file is read, if the whole Delaunay
c     triangulation and Voronoi diagram have been computed,
c     triangulation information is saved as follows. Arrays and
c     variables appearing here are described below. Here iflag is
c     always equal to zero signifying that the whole Delaunay
c     triangulation and Voronoi diagram have been computed.
c
c     open(unit=16,file=tn(idcr))
c     write(16,*) iflag
c     write(16,*) n, ivnxt, nvmz, icsfig, idyn
c     write(16,130)((icon(i,j),i=1,6),j=1,ivnxt)
c     write(16,140)(is(i),i=1,n)
c 130 format(6i10)
c 140 format(5i10)
c     close(unit=16)
c
c     Information about Voronoi vertices is saved as follows. Arrays
c     and variables appearing here are described below.
c
c     open(unit=18,file=vn(idcr))
c     write(18,*) n, ivnxt, icsfig, larg
c     write(18,190) (id(i),i=1,ivnxt)
c     do 180 i=1,ivnxt
c        write(18,*) vx(i),vy(i)
c 180 continue
c 190 format(40(1x,i1))
c     close(unit=18)
c
c     If the Delaunay triangulation and Voronoi diagram have been
c     updated on a local basis, then triangulation information is
c     saved as follows for the purpose of minimizing the size of
c     the output file. Arrays and variables appearing here are
c     described below. Here iflag is always equal to one signifying
c     that the Delaunay triangulation and Voronoi diagram have been
c     updated on a local basis.
c
c     open(unit=16,file=tn(idcr))
c     write(16,*) iflag
c     write(16,*) n, ivnxt, nvmz, icsfig
c     write(16,*) nd, izn, ian, ibn, icn, icnxt
c     write(16,140)(idl(i),i=1,nd)
c     write(16,140)(iz(i),i=1,izn)
c     write(16,140)(ia(i),i=1,ian)
c     write(16,140)(ib(i),i=1,ibn)
c     if(icn.ne.0) write(16,140)(ic(i),i=1,icn)
c     write(16,130)((icon(i,ib(j)),i=1,6),j=1,ibn)
c     if(icn.ne.0) write(16,130)((icon(i,ic(j)),i=1,6),j=1,icn)
c     write(16,140)(is(iz(i)),i=1,izn)
c     close(unit=16)
c
c     On a local basis vertices of Voronoi diagram are saved as follows.
c     Arrays and variables appearing here are described below.
c
c     open(unit=18,file=vn(idcr))
c     write(18,*) n, ivnxt, icsfig, larg
c     write(18,1050) (id(ib(i)),i=1,ibn)
c     do 1040 i=1,ibn
c        write(18,*) vx(ib(i)),vy(ib(i))
c1040 continue
c1050 format(40(1x,i1))
c     close(unit=18)
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            Note: Even though in what follows and throughtout the
c                  program the possibility of having inserted edges
c                  is treated as a real possibility, in this program
c                  that possibility does not exist. The comments and
c                  programming addressing that possibility are still
c                  valid even though the possibility never occurs.
c                  They should have no effect on either the correctness
c                  of the documention or the efficiency and correctness
c                  of the program.
c
c     idyn   Number of files or data sets that contain input points
c            making up dynamic data. The number of points in a file
c            is the same for all files. Points are not allowed to
c            be duplicated in a file.
c
c     idcr   Index of current input points file in dynamic data.
c    
c     nv     Number of input points in each data set in dynamic
c            data for which Delaunay triangulations and corresponding
c            Voronoi diagrams are computed. Value should be the same
c            for all data sets in dynamic data.
c
c     nd     Number of points or vertices to be deleted (because they
c            have moved) from a given Delaunay triangulation. The
c            same points with their new coordinates are then added
c            again to the Delaunay triangulation.
c
c     ndlim  Upper bound on number of vertices that can be deleted
c            and added again to the current triangulation before
c            program switches to just computing the whole next
c            triangulation from scratch. One fourth of the total
c            number of vertices may be a good value for it.
c
c     iflag  Flag that indicates whether all or part of current
c            Delaunay triangulation or Voronoi diagram was computed.
c            Equals 0 if all was computed. Equals 1 if only a part.
c
c     icnxt  Number of triangles in current Delaunay triangulation
c            or number of Voronoi vertices in the corresponding
c            Voronoi diagram.
c
c     ivnxt  Number of current and previous triangles whose information
c            appears in array icon, i. e. length of array icon.
c
c     iabe   Number of previously used triangles available for use in
c            the current triangulation.
c
c     ian    Number of triangles that have been eliminated from
c            previous triangulation.
c
c     ibn    Number of triangles that have been added to current
c            triangulation.
c
c     icn    Number of triangles that are in both the previous and
c            current triangulation adjacent to a triangle that has
c            been added to current triangulation.
c
c     izn    Number of points or vertices that are vertices of
c            triangles that have been added to current 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     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 a 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 id, it, vx,
c            vy (described below). Value is the maximum allowed number
c            of current and previous triangles whose information
c            appears in array icon. Value at least that of ivnxt.
c            Value should be at least two times that of nmax plus
c            some more to account for eliminated Delaunay triangles
c            (Voronoi vertices) in the previous Delaunay triangulation
c            and Voronoi diagram.
c
c     ndmax  Parameter which is the dimension of array idl (described
c            below). Value is the maximum allowed number of vertices
c            to be deleted (and then added) from a given Delaunay
c            triangulation. Value at least that of nd (described above).
c
c     nfmax  Parameter which is the dimension of several arrays used
c            internally during the deletion of a point (see below).
c            Value at some multiple of the square root of nmax
c            (described above) seems to work.
c
c     njmax  Parameter which is the dimension of arrays iave, ia, ib
c            (see below). Value at an integral multiple of ndmax
c            (described above) seems to work.
c
c     nkmax  Parameter which is the dimension of array ic (see below).
c            Value at an integral multiple of ndmax (described above)
c            seems to work.
c
c     nlmax  Parameter which is the dimension of array iz (see below).
c            Value at an integral multiple of ndmax (described above)
c            seems to work.
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     x1     Array similar to array x (described above) for previous
c            data set of input points.
c
c     y1     Array similar to array y (described above) for previous
c            data set of input points.
c
c     x2     Array similar to array x (described above) for current
c            data set of input points.
c
c     y2     Array similar to array y (described above) for current
c            data set of input points.
c
c     vx     Array that contains on output x-coordinates of
c            circumcenters (Voronoi vertices) of triangles 1 through
c            ivnxt. nvmax is the 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. nvmax is the 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, then
c            icon(i,j), i = 1, ..., 6, provides information about
c            triangle j that is a triangle in either a previous
c            triangulation or in the current one. Assuming it is
c            in the current one, icon(i,j), i = 1, 2, 3, are the
c            triangle numbers or indices for the three triangles
c            adjacent to triangle j in the order in which they appear
c            in a counter-clockwise direction around triangle j. Given
c            i, i=1, 2, or 3, if icon(i,j) is negative and not equal
c            to nvmz, then -icon(i,j) is the actual index for the ith
c            adjacent triangle to triangle j and the negativity of
c            icon(i,j) simply implies that the edge shared by triangle
c            j and triangle -icon(i,j) is part of an inserted line
c            segment. Given i, i=1, 2, or 3, if icon(i,j) equals either
c            zero or nvmz then the ith adjcent triangle to triangle j
c            does not exist and the edge of triangle j that triangle j
c            would 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. Finally, icon(i,j), i = 4, 5, 6, are the vertex
c            numbers or indices for the three vertices of triangle j in
c            the order in which they appear in a counterclockwise
c            direction around triangle j. In addition, assuming without
c            any loss of generality that for i=1, 2, 3, icon(i,j) is
c            positive, for i=4, 5, 6, icon(i,j) is the vertex number
c            for the vertex of triangle j that is not a vertex for
c            triangle icon(i-3,j).
c
c     id     Array with dimension nvmax. For i = 1,..., ivnxt, id(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 id(i) will equal 0.
c
c     it     Array with dimension nvmax. For i = 1,..., ivnxt, it(i)
c            equals 1 if triangle i is in both the current and previous
c            triangulations. Equals -1 if triangle i is in previous
c            triangulation but not in current one. Equals 2 if it is
c            in current triangulation but not in the previous one.
c            Equals 3 (instead of 1) if it is in both the current and
c            previous triangulations but it is adjacent to a triangle j
c            in the current triangulation with it(j) equal to 2.
c            Equals -2 if it is in neither one but at one point it was
c            in the current one. Equals 0 if it is not and has not
c            been in either one.
c
c     idl    Array with dimension ndmax. For i = 1,..., nd, idl(i)
c            is the ith vertex to be deleted (and then added) from a
c            given Delaunay triangulation.
c
c     ia     Array of dimension njmax. For i = 1,..., ian, ia(i)
c            is the ith triangle that has been eliminated from
c            previous triangulation. icon(j,ia(i)), j=1, 2, 3,
c            may not be what they were in previous triangulation.
c
c     ib     Array of dimension njmax. For i = 1,..., ibn, ib(i)
c            is the ith triangle that has been added to current
c            triangulation. icon(j,ib(i)), j=1,...,6, are what
c            they are supposed to be in current triangulation.
c
c     ic     Array of dimension nkmax. For i = 1,..., icn, ic(i)
c            is the ith triangle that is in both the previous and
c            current triangulation adjacent to a triangle that has
c            been added to current triangulation. icon(j,ic(i)),
c            j=1, 2, 3, may not be what they were in previous
c            triangulation, but icon(j,ic(i)), j =1,...,6, are what
c            they are supposed to be in current triangulation.
c
c     iz     Array of dimension nlmax. For i = 1,...,izn, iz(i)
c            is the ith point or vertex that is a vertex of a
c            triangle that has been added to current triangulation.
c
c     iave   Array of dimension njmax. For i = 1,..., iabe, iave(i)
c            is the ith previously used triangle available for use
c            in the current triangulation.
c
c     ix
c     iy
c     ix2
c     iy2    Arrays of dimension nmax used internally.
c
c     ifun
c     irun
c     isun
c     izun
c     igun
c     iwun   Arrays of dimension nfmax used internally.
c
c
c***********************************************************************
c
      program main
c
      integer nmax, nvmax, ndmax
      integer nfmax, njmax, nkmax, nlmax
      parameter (nmax = 160000, nvmax = 2.1*nmax, ndmax = 40000)
      parameter (nfmax = 5*400, njmax = 3*ndmax)
      parameter (nkmax = ndmax, nlmax = 2*ndmax)
c
      double precision x(nmax), y(nmax), vx(nvmax), vy(nvmax)
      double precision x1(nmax), y1(nmax), x2(nmax), y2(nmax)
      integer ix(nmax), iy(nmax), ix2(nmax), iy2(nmax)
      integer is(nmax), icon(6,nvmax), id(nvmax), it(nvmax), iz(nlmax)
      integer idl(ndmax), iave(njmax), ia(njmax), ib(njmax), ic(nkmax)
      integer ifun(nfmax), irun(nfmax), isun(nfmax)
      integer izun(nfmax), igun(nfmax), iwun(nfmax)
c
      integer idyn, idcr, n, nv, icnxt, ivnxt, nvmz, nd, ndlim, iflag
      integer iabe, ian, ibn, icn, izn, icsfig, mhalf, mfull, isclp(2)
      double precision epf
      double precision r215, deps, dscle, dfull, dfill
c
      character*3 nu(100)
      character*10 pn(100), tn(100), vn(100)
c
      integer i, j, k, nvmn, nvmp, ileft, idupl, inpro, ivi
      integer ndl, nde, ndv, nds
      double precision xcor, ycor, xmax, xmin, ymax, ymin
c
      integer ng, ifir, isec, ithi, larg, izero
      double precision derro, dist1, dist2, dist3
      double precision dista, disti, diffd, dnom
c
      izero = 0
      derro = 0.0
c
      write(*,*)' '
      write(*,*)'This program is for computing Voronoi diagrams ',
     *          'of dynamic data.'
      write(*,*)' '
      write(*,*)'Currently program is set up to handle 100 or fewer ',
     *          'input files'
      write(*,*)'making up dynamic data.'
      write(*,*)' '
      write(*,*)'Enter number of input points files that make up ',
     *          'dynamic data:'
      read(5,*) idyn
      if(idyn.lt.2 .or. idyn.gt.100) stop 5
      if(idyn.eq.izero) izero = 0
c
c     initialize for names of dynamic input data files and output files
c
      do 5 j = 1, 100
        do 3 i =1,3
           nu(j)(i:i) = ' '
    3   continue
    5 continue

      nu(1) = '1'
      nu(2) = '2'
      nu(3) = '3'
      nu(4) = '4'
      nu(5) = '5'
      nu(6) = '6'
      nu(7) = '7'
      nu(8) = '8'
      nu(9) = '9'
c
      do 10 j = 0, 9
         i=j+10
         nu(i)(1:1) = '1'
         i=j+20
         nu(i)(1:1) = '2'
         i=j+30
         nu(i)(1:1) = '3'
         i=j+40
         nu(i)(1:1) = '4'
         i=j+50
         nu(i)(1:1) = '5'
         i=j+60
         nu(i)(1:1) = '6'
         i=j+70
         nu(i)(1:1) = '7'
         i=j+80
         nu(i)(1:1) = '8'
         i=j+90
         nu(i)(1:1) = '9'
   10 continue
c
      do 15 j = 1, 9
         i=j*10
         nu(i+0)(2:2) = '0'
         nu(i+1)(2:2) = '1'
         nu(i+2)(2:2) = '2'
         nu(i+3)(2:2) = '3'
         nu(i+4)(2:2) = '4'
         nu(i+5)(2:2) = '5'
         nu(i+6)(2:2) = '6'
         nu(i+7)(2:2) = '7'
         nu(i+8)(2:2) = '8'
         nu(i+9)(2:2) = '9'
   15 continue
c
      nu(100) = '100'
c
      do 25 j = 1, 100
         do 20 i = 1, 10
            pn(j)(i:i) = ' '
            tn(j)(i:i) = ' '
            vn(j)(i:i) = ' '
   20    continue
c
         pn(j)(1:2) = 'pn'
         tn(j)(1:2) = 'tn'
         vn(j)(1:2) = 'vn'
c
         pn(j)(3:5) = nu(j)
         tn(j)(3:5) = nu(j)
         vn(j)(3:5) = nu(j)
c
   25 continue
c
      write(*,*)' '
      write(*,*)'Enter icsfig, i. e. number 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

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
c
      epf = 0.01d0

c
c*****************************************************************
c
c     read first points file
c
      idcr = 1
      write(*,*)' '
      write(*,*)'Index of current points file = ',idcr
      write(*,*)' '
      write(*,*)'Reading first points file (x-,y-coordinates) ...'

      open(unit=11,file=pn(idcr))
      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
      x1(nv) = xcor
      y1(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 15
      write(*,*)' '
      write(*,*)'Number of points read =',nv
      n=nv
      ndlim = nv/2

      write(*,*)' '
      write(*,*)' xmin = ',xmin,' xmax = ',xmax
      write(*,*)' ymin = ',ymin,' ymax = ',ymax
c
c
c*****************************************************************
c
c     transform coordinates into integer decomposition
c
  120 continue
c     write(*,*)' '
c     write(*,*)'Transforming all point coordinates from current ',
c    *          'input file into integer decompositions ...'

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

c
c*****************************************************************
c
c     initialize triangle pointer and triangle status array
c
      do 124 i = 1, n
         is(i) = 0
  124 continue
c
      do 126 i = 1, nvmax
         it(i) = 0
  126 continue
c
      write(*,*)' '
      write(*,*)'Computing entire Delaunay triangulation for ',
     *          'current input file ...'
c
      ian = 0
      ibn = 0
      call tricmp(x, y, ix, iy, ix2, iy2, is, icon, it, iave, ia, ib,
     *            epf, mhalf, mfull, isclp, nvmz, nv, iabe, ian, ibn,
     *            njmax, nvmax, ivnxt, idupl, inpro, r215, deps, ileft)
c
      if(ian.ne.0) stop 20
      if(ibn.ne.ivnxt) stop 25
      if(iabe. ne.1) stop 30
      if(idupl.ne.0) stop 35
      if(inpro.ne.0) stop 40
c
      iflag = 0
      icnxt = ivnxt
      do 128 i = 1, ivnxt
         it(i) = 1
  128 continue
c
c**********************************************************************
c
c     compress array icon
c
c     call comprs(is, icon, id, it, n, ivnxt, nvmax, nvmz)
c
c**********************************************************************
c
c     save information about current Delaunay triangulation
c
c     write(*,*)' '
c     write(*,*)'Saving information about current Delaunay ',
c    *          'triangulation ... '

      open(unit=16,file=tn(idcr))
      write(16,*) iflag
      write(16,*) n, ivnxt, nvmz, icsfig, idyn
      write(16,130)((icon(i,j),i=1,6),j=1,ivnxt)
      write(16,140)(is(i),i=1,n)
  130 format(6i10)
  140 format(5i10)
      close(unit=16)

c
c********************************************************************
c
c     compute and save circumcenters (x-,y-coordinates) of triangles,
c     i. e. Voronoi vertices in first Voronoi diagram
c
c     write(*,*)' '
c     write(*,*)'Computing circumcenters (x-,y-coordinates) ',
c    *          'of triangles,'
c     write(*,*)'i. e. Voronoi vertices in current Voronoi diagram ...'

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

      do 150 i = 1, ivnxt
         if(it(i).eq.0) stop 42
         if(id(i).eq.0) go to 150
         go to 160
  150 continue
c     write(*,*)' '
c     write(*,*)'All Voronoi vertices marked as too large.'
      go to 175
  160 continue
      ng = 0
      derro = 0.0d0
      larg = 0
      do 170 i = 1, ivnxt
         if(it(i).eq.0) stop 43
         if(id(i).eq.0) then
            larg = larg+1
            go to 170
         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 45
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
  170 continue
c     write(*,*)' '
c     write(*,*)'Number of Voronoi vertices marked as too ',
c    *          'large = ',larg
      write(*,*)' '
      write(*,*)'Voronoi distance error = ',derro
      write(*,*)'(This number should be close to zero).'
  175 continue

c     write(*,*)' '
c     write(*,*)'Saving circumcenters/Voronoi vertices ... '

      open(unit=18,file=vn(idcr))
      write(18,*) n, ivnxt, icsfig, larg
      write(18,190) (id(i),i=1,ivnxt)
      do 180 i=1,ivnxt
         write(18,*) vx(i),vy(i)
  180 continue
  190 format(40(1x,i1))
      close(unit=18)

c
c*********************************************************************
c
c     test consistency of triangulation
c
c     write(*,*)' '
c     write(*,*)'Testing consistency of current Delaunay ',
c    *          'triangulation ...'
      call contst(icon, it, is, id, nv, ivnxt, icnxt, nvmz)
c
c*********************************************************************
c
c     test circle criterion
c
c     write(*,*)' '
c     write(*,*)'Testing circle criterion ... '
      call circri(x, y, ix, iy, ix2, iy2, icon, it, ib, ibn, mhalf,
     *            mfull, isclp, epf)
c
c*********************************************************************
c
c     test orientation of triangles
c
c     write(*,*)' '
c     write(*,*)'Testing orientation of triangles ... '
      call orient(x, y, ix, iy, ix2, iy2, icon, it, ib, ibn, mhalf,
     *            mfull, isclp, epf)
c
c*********************************************************************
c
c     test convexity of union of triangles
c
c     write(*,*)' '
c     write(*,*)'Testing convexity of union of triangles ... '
      call convex(x, y, ix, iy, ix2, iy2, icon, it, is, nv, ivnxt,
     *             nvmz, mhalf, mfull, isclp, epf)
c
c********************************************************************
c
c     computations are completed
c
      write(*,*)' '
      write(*,*)'Computations completed for current data set.'
      write(*,*)' '
c     write(*,*)'Number of vertices in current Delaunay ',
c    *          'triangulation: ', nv
      write(*,*)'Number of Delaunay triangles or Voronoi vertices: ',
     *           ivnxt
      write(*,*)' '
c
c*****************************************************************
c
c     read next points file
c
  195 continue
      if(idcr.ge.idyn) go to 1060
      idcr = idcr+1
      write(*,*)' '
      write(*,*)'Index of current points file = ',idcr
      write(*,*)' '
      write(*,*)'Reading current points file (x-,y-coordinates) ...'
      open(unit=11,file=pn(idcr))
      nv = 0
  200 continue
      read(11,*,end=210) xcor, ycor
      nv = nv + 1
      if(nv .gt. nmax) stop 50
      x2(nv) = xcor
      y2(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 200
  210 continue
      close(unit=11)
      if(n.ne.nv) stop 55

      write(*,*)' '
      write(*,*)' xmin = ',xmin,' xmax = ',xmax
      write(*,*)' ymin = ',ymin,' ymax = ',ymax
c
c**********************************************************************
c
c     identify points to be deleted and then added
c
      nd = 0
      do 250 i = 1, n
         if(x1(i).eq.x2(i) .and. y1(i).eq.y2(i)) go to 250
         nd = nd + 1
         if(nd.gt.ndmax) stop 60
         if(nd.gt.ndlim) go to 260
         idl(nd) = i
  250 continue
      if(nv-nd.lt.3) go to 260
      go to 270
c
  260 continue
      write(*,*)' '
      write(*,*)'Too many points that have moved as identified by ',
     *          'comparing previous'
      write(*,*)'and current points files in dynamic data '
      write(*,*)' '
      write(*,*)'Entire Delaunay triangulation for current points ',
     *          'file will be'
      write(*,*)'computed from scratch'
c
      do 265 i = 1, n
         x(i) = x2(i)
         y(i) = y2(i)
         x1(i) = x(i)
         y1(i) = y(i)
  265 continue
      go to 120
c
c**********************************************************************
c
c     delete set of points from triangulation before they move
c
  270 continue
      write(*,*)' '
      write(*,*)'Updating last Delaunay triangulation for ',
     *          'current input file ...'
      write(*,*)' '
      write(*,*)'Number of points that have moved as identified by ',
     *          'comparing previous'
      write(*,*)'and current points files in dynamic data = ',nd
      if(nd.eq.0) stop 65
c     write(*,*)' '
c     write(*,*)'Deleting vertices from triangulation ',
c    *          'before they move ...'

      ian = 0
      ibn = 0
      call pntdel(x, y, idl, ix, iy, ix2, iy2, is, icon, it, iave, ia,
     *            ib, ifun, irun, isun, izun, igun, iwun, epf, mhalf,
     *            mfull, isclp, nvmz, nv, nd, iabe, ian, ibn, ivnxt,
     *            nfmax, njmax, nvmp, nvmn, ndl, nde, ndv, nds,
     *            ileft, nvmax)
      if(ndl.ne. nd. or. nde.ne.0 .or. ndv.ne.0 .or. nds.ne.0) stop 70
      if(ian.eq.0 .or. ibn.eq.0) stop 75
c
c     WRITE(*,*)'270 IABE=',IABE
c
c
c**********************************************************************
c
c     add set of points again into triangulation as they have moved
c
c     write(*,*)' '
c     write(*,*)'Adding vertices into triangulation ',
c    *          'that have moved ...'
c
      ndl = nd
      nd = ndl
      do 300 j = 1, ndl
         i = idl(j)
         x(i) = x2(i)
         y(i) = y2(i)
         x1(i) = x(i)
         y1(i) = y(i)
  300 continue
c
      call triadd(x, y, idl, ix, iy, ix2, iy2, is, icon, it, iave, ia,
     *            ib, epf, icsfig, mhalf, mfull, isclp, nvmz, nv, nd,
     *            iabe, ian, ibn, njmax, nvmax, ivnxt, idupl, inpro,
     *            ileft, dscle, dfull, dfill)
c
      if(idupl.ne.0) stop 80
      if(inpro.ne.0) stop 85
      if(ian.eq.0 .or. ibn.eq.0) stop 90
c
c**********************************************************************
c
c     compress array ib, update array it accordingly, test arrays, etc.
c
c     WRITE(*,*)'IVNXT =',IVNXT
c     WRITE(*,*)'IT =',(IT(I),I=1,IVNXT)
      j = 0
      do 400 i = 1, ibn
         k = ib(i)
         if(it(k).eq.-2) then
            it(k) = 0
         elseif(it(k).eq.2) then
            j = j+1
            ib(j) = k
         else
            stop 95
         endif
  400 continue
      ibn = j
      if(ibn.eq.0) stop 100
c
      icn = 0
      izn = 0
      do 450 i = 1, ibn
         k = ib(i)
         if(it(k).ne.2) stop 105
         do 420 j = 1, 3
            ivi = icon(j,k)
            if(ivi .eq. 0 .or. ivi .eq. nvmz) go to 420
            ivi = iabs(ivi)
            if(it(ivi).eq.2 .or. it(ivi).eq.3) go to 420
            if(it(ivi).ne.1) stop 110
            it(ivi) = 3
            icn = icn + 1
            if(icn .gt. nkmax) stop 115
            ic(icn) = ivi
  420    continue
         do 440 j = 4, 6
            ivi = icon(j,k)
            if(is(ivi).lt.0) go to 440
            if(is(ivi).eq.0) stop 120
            is(ivi) = -is(ivi)
            izn = izn + 1
            if(izn .gt. nlmax) stop 125
            iz(izn) = ivi
  440    continue
  450 continue
      if(izn.eq.0) stop 130
      ivi = icn
      icn = ivi
      ivi = izn
      izn = ivi
c     WRITE(*,*)'IS ='
c     WRITE(*,140)(IS(I),I=1,N)
c     WRITE(*,*)'IZN =',IZN
c     WRITE(*,*)'IZ =',(IZ(I),I=1,IZN)
c
      do 470 i = 1, izn
         k = iz(i)
         if(is(k).ge.0) stop 135
         is(k) = -is(k)
  470 continue
c
      if(iabe.eq.1) go to 550
      do 500 i = 2, iabe
         k = iave(i)
         if(it(k).ne.0) stop 140
  500 continue
  550 continue
c
c     WRITE(*,*)'IABE =',IABE
c     WRITE(*,*)'IAVE =',(IAVE(I),I=1,IABE)
c     WRITE(*,*)'IVNXT =',IVNXT
c     WRITE(*,*)'IT =',(IT(I),I=1,IVNXT)
c
      do 600 i = 1, ian
         k = ia(i)
         if(it(k).ne.-1) stop 145
         it(k) = 0
         iabe = iabe + 1
         if(iabe .gt. njmax) stop 150
         iave(iabe) = k
  600 continue
c
      do 700 i = 1, ibn
         k = ib(i)
         if(it(k).ne.2) stop 155
         it(k) = 1
  700 continue
c
      if(icn.eq.0) go to 850
      do 800 i = 1, icn
         k = ic(i)
         if(it(k).ne.3) stop  160
         it(k) = 1
  800 continue
  850 continue
c
c     WRITE(*,*)'OLD ICNXT =',ICNXT
      iflag = 1
      if(icnxt + ibn + iabe - ian - 1 .ne. ivnxt) stop 165
      icnxt = icnxt - ian + ibn
c     WRITE(*,*)'NEW ICNXT =',ICNXT
c
c**********************************************************************
c
c     WRITE(*,*)'IABE =',IABE
c     WRITE(*,*)'IAVE =',(IAVE(I),I=1,IABE)
c     WRITE(*,*)'IVNXT =',IVNXT
c     WRITE(*,*)'IT =',(IT(I),I=1,IVNXT)
c     WRITE(*,*)'IAN =',IAN
c     WRITE(*,*)'IA =',(IA(I),I=1,IAN)
c     WRITE(*,*)'IBN =',IBN
c     WRITE(*,*)'IB =',(IB(I),I=1,IBN)
c     WRITE(*,*)'ICN =',ICN
c     IF(ICN.NE.0) WRITE(*,*)'IC =',(IC(I),I=1,ICN)
c     WRITE(*,*)'ICON ='
c     WRITE(*,130)((ICON(I,J),I=1,6),J=1,IVNXT)
c     WRITE(*,*)'IS ='
c     WRITE(*,140)(IS(I),I=1,N)
c
c**********************************************************************
c
c     compress array icon
c
c     call comprs(is, icon, id, it, n, ivnxt, nvmax, nvmz)
c
c**********************************************************************
c
c     save information about current Delaunay triangulation
c
c     write(*,*)' '
c     write(*,*)'Saving partial information about current Delaunay ',
c    *          'triangulation ... '

      open(unit=16,file=tn(idcr))
      write(16,*) iflag
      write(16,*) n, ivnxt, nvmz, icsfig
      write(16,*) nd, izn, ian, ibn, icn, icnxt
      write(16,140)(idl(i),i=1,nd)
      write(16,140)(iz(i),i=1,izn)
      write(16,140)(ia(i),i=1,ian)
      write(16,140)(ib(i),i=1,ibn)
      if(icn.ne.0) write(16,140)(ic(i),i=1,icn)
      write(16,130)((icon(i,ib(j)),i=1,6),j=1,ibn)
      if(icn.ne.0) write(16,130)((icon(i,ic(j)),i=1,6),j=1,icn)
      write(16,140)(is(iz(i)),i=1,izn)
      close(unit=16)
c
c********************************************************************
c
c     compute and save circumcenters (x-,y-coordinates) of triangles,
c     i. e. Voronoi vertices in current triangulation
c
c     write(*,*)' '
c     write(*,*)'Computing circumcenters (x-,y-coordinates) ',
c    *          'of triangles,'
c     write(*,*)'i. e. Voronoi vertices in current Voronoi diagram ...'

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

      do 1010 j = 1, ibn
         i = ib(j)
         if(it(i).eq.0) stop 167
         if(id(i).eq.0) go to 1010
         go to 1015
 1010 continue
c     write(*,*)' '
c     write(*,*)'All Voronoi vertices marked as too large.'
      go to 1030
 1015 continue
      ng = 0
      derro = 0.0d0
      larg = 0
      do 1020 j = 1, ibn
         i = ib(j)
         if(it(i).eq.0) stop 168
         if(id(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 170
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
c     write(*,*)' '
c     write(*,*)'Number of Voronoi vertices marked as too ',
c    *          'large = ',larg
      write(*,*)' '
      write(*,*)'Voronoi distance error = ',derro
      write(*,*)'(This number should be close to zero).'
 1030 continue

c     write(*,*)' '
c     write(*,*)'Saving circumcenters/Voronoi vertices ... '

      open(unit=18,file=vn(idcr))
      write(18,*) n, ivnxt, icsfig, larg
      write(18,1050) (id(ib(i)),i=1,ibn)
      do 1040 i=1,ibn
         write(18,*) vx(ib(i)),vy(ib(i))
 1040 continue
 1050 format(40(1x,i1))
      close(unit=18)

c
c*********************************************************************
c
c     test consistency of triangulation
c
c     write(*,*)' '
c     write(*,*)'Testing consistency of current Delaunay ',
c    *          'triangulation ...'
c     call contst(icon, it, is, id, nv, ivnxt, icnxt, nvmz)
      call cnstst(icon, it, is, ib, iz, ibn, izn, nvmz)
c
c*********************************************************************
c
c     test circle criterion
c
c     write(*,*)' '
c     write(*,*)'Testing circle criterion ... '
      call circri(x, y, ix, iy, ix2, iy2, icon, it, ib, ibn, mhalf,
     *            mfull, isclp, epf)
c
c*********************************************************************
c
c     test orientation of triangles
c
c     write(*,*)' '
c     write(*,*)'Testing orientation of triangles ... '
      call orient(x, y, ix, iy, ix2, iy2, icon, it, ib, ibn, mhalf,
     *            mfull, isclp, epf)
c
c*********************************************************************
c
c     test convexity of union of triangles
c
c     write(*,*)' '
c     write(*,*)'Testing convexity of union of triangles ... '
c     call convex(x, y, ix, iy, ix2, iy2, icon, it, is, nv, ivnxt,
c    *             nvmz, mhalf, mfull, isclp, epf)
      call cnvtst(x, y, ix, iy, ix2, iy2, icon, is, iz, izn, nvmz,
     *            mhalf, mfull, isclp, epf)
c
c********************************************************************
c
c     computations are completed
c
      write(*,*)' '
      write(*,*)'Computations completed for current data set.'
      write(*,*)' '
c     write(*,*)'Number of vertices in current Delaunay ',
c    *          'triangulation: ',nv
      write(*,*)'Number of Delaunay triangles or Voronoi vertices: ',
     *           icnxt
      write(*,*)' '
      go to 195
c
 1060 continue
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 1100 i = 1, nv
c         write(*,*)' pnt: ',i,' pnter: ',is(i),
c    *               ' x: ',x(i),' y: ',y(i)
c1100 continue
c     write(*,*)'Number of triangles: '
c     write(*,*) ivnxt
c     write(*,*)'Triangles: '
c     do 1200 i = 1, ivnxt
c         write(*,1400) i, (icon(j,i),j=1,6)
c1200 continue
c1400 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 280
      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(*,*)'Number of significant figures of largest ',
     *             'nondecimal part of'
         write(*,*)'a point coordinate appears to be greater than 9.'
         write(*,*)'Program is terminated.'
         stop 285
      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 290
      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 295
      call decomp(isclp, isgcl, isclu, mhalf)
      if(isgcl.ne.1) stop 297
c
c     test lengths of x, y-coordinates, shift and make them integers
c
      dfull = dble(mfull)
      if(dscle.lt.deps) stop 300
      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 305
         ix(i) = idint(x(i))
         if(iabs(ix(i)).ge.mfull) stop 310
         decml = (x(i) - dint(x(i)))*dscle
         if(dabs(decml).ge.dfull) stop 312
         ix2(i) = idnint(decml)
         if(iabs(ix2(i)).ge.mfull) stop 315
         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 320
         iy(i) = idint(y(i))
         if(iabs(iy(i)).ge.mfull) stop 325
         decml = (y(i) - dint(y(i)))*dscle
         if(dabs(decml).ge.dfull) stop 327
         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(*,*)'Number of significant figures of nondecimal',
     *             'part of a'
         write(*,*)'new point coordinate appears to be greater than 9.'
         write(*,*)'Program is terminated.'
         stop 335
      endif
      itsfig = ibsfig + icsfig
c     WRITE(*,*)'ITSFIG=',ITSFIG,' IBSFIG=',IBSFIG,' ICSFIG=',ICSFIG
      if(itsfig.gt.14) then
         write(*,*)'Total number of significant figures of a new ',
     *             'point coordinate'
         write(*,*)'appears to be greater than 14.'
         stop 340
      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 350
      x(n) = idint(xi(n))
      if(iabs(x(n)).ge.mfull) stop 355
      decml = (xi(n) - dint(xi(n)))*dscle
      if(dabs(decml).ge.dfull) stop 360
      x2(n) = idnint(decml)
      if(iabs(x2(n)).ge.mfull) stop 365
      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 370
      y(n) = idint(yi(n))
      if(iabs(y(n)).ge.mfull) stop 375
      decml = (yi(n) - dint(yi(n)))*dscle
      if(dabs(decml).ge.dfull) stop 380
      y2(n) = idnint(decml)
      if(iabs(y2(n)).ge.mfull) stop 385
      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, it, iave,
     *                  ia, ib, epf, mhalf, mfull, isclp, nvmz, nv,
     *                  iabe, ian, ibn, njmax, nvmax, ivnxt, idupl,
     *                  inpro, r215, deps, ileft)
c
      double precision x(*), y(*)
      integer ix(*), iy(*), ix2(*), iy2(*)
      integer is(*), icon(6,*), it(*), iave(*), ia(*), ib(*)
      double precision epf, r215, deps
      integer mhalf, mfull, nvmz, nv
      integer iabe, ian, ibn, njmax, nvmax, ivnxt, 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
      double precision xtemp, ytemp, delx, dely
      double precision xmax, xmin, dotx, dnux, dot
      integer i, i1, i2, i3, ip, 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 400
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 410
c
      call prpvec(ix, iy, ix2, iy2, i1, i2, mhalf, mfull, isclp,
     *            iax, isgax, ikax, iay, isgay, ikay)
      if(isgay.le.0) stop 420
      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 430
      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) stop 440
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 450
      if(isga .gt. 0) go to 400
      ip = i2
      i2 = i3
      i3 = ip
c
  400 continue
c     write(*,*)' '
c     write(*,*)'Initial triangle has vertices:',i1,i2,i3
      ivnxt = 1
      if(ivnxt.gt.nvmax) stop 460
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
      if(iabe.gt.njmax) stop 470
      iave(1)=ivnxt+1
      it(1) = 2
      ibn = 1
      ib(1) = 1
      ileft=i1
      insrt=1
      idupl=0
      inpro=0
      if(nv.eq.3) go to 480
c
      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, it, iave, ia,
     *               ib, epf, mhalf, mfull, isclp, nvmz, ileft, iright,
     *               nv, iabe, ian, ibn, njmax, isite, 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/10000)*10000)
     *      write(*,*)'Number of points inserted =',isite
  450 continue
c
  480 continue
      return
      end
*TRIADD
c
c     driver subroutine triadd to -
c
c     add points into current Delaunay triangulation
c
      subroutine triadd(x, y, idl, ix, iy, ix2, iy2, is, icon, it,
     *                  iave, ia, ib, epf, icsfig, mhalf, mfull, isclp,
     *                  nvmz, nv, nd, iabe, ian, ibn, njmax, nvmax,
     *                  ivnxt, idupl, inpro, ileft, dscle, dfull, dfill)
c
      double precision x(*), y(*)
      integer idl(*), ix(*), iy(*), ix2(*), iy2(*)
      integer is(*), icon(6,*), it(*), iave(*), ia(*), ib(*)
      double precision epf, dscle, dfull, dfill
      integer icsfig, mhalf, mfull, isclp(*)
      integer nvmz, nv, nd, iabe, ian, ibn, njmax, nvmax, ivnxt
c
      integer ileft, iright, insrt, isst
      integer idupl, inpro, isite, itype, ipoint
c
      insrt=1
      idupl=0
      inpro=0
c
      write(*,*)' '
      do 580 isst = 1, nd
         isite = idl(isst)
         if(is(isite).eq.nvmz) then
            inpro=inpro+1
            go to 580
         endif
         if(is(isite).ne.0) stop 500
         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, it, iave, ia,
     *               ib, epf, mhalf, mfull, isclp, nvmz, ileft, iright,
     *               nv, iabe, ian, ibn, njmax, isite, 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(isst.le.(isst/1000)*1000)
     *      write(*,*)'Number of points inserted =',isst
  580 continue
      return
      end
*PNTDEL
c
c     driver subroutine pntdel to -
c
c     delete points or vertices in current Delaunay triangulation
c
      subroutine pntdel(x, y, idl, ix, iy, ix2, iy2, is, icon, it,
     *                  iave, ia, ib, ifun, irun, isun, izun, igun,
     *                  iwun, epf, mhalf, mfull, isclp, nvmz,
     *                  nv, nd, iabe, ian, ibn, ivnxt, nfmax, njmax,
     *                  nvmp, nvmn, ndl, nde, ndv, nds, ileft, nvmax)
c
      double precision x(*), y(*)
      integer idl(*), ix(*), iy(*), ix2(*), iy2(*)
      integer is(*), icon(6,*), it(*), iave(*), ia(*), ib(*)
      integer ifun(*), irun(*), isun(*), izun(*)
      integer igun(*), iwun(*)
      double precision epf
      integer mhalf, mfull, isclp(*)
      integer nvmz, nv, nd, iabe, ian, ibn, ivnxt, nvmax
      integer nfmax, njmax
      integer nvmp, nvmn, ndl, nde, ndv, nds
c
      integer ileft, iright
      integer i, isucs, isnew
c
      ndl=0
      nde=0
      ndv=0
      nds=0
c
      write(*,*)' '
      do 720 i = 1, nd
         isucs = -1
         iright = idl(i)
         call deleli(x, y, ix, iy, ix2, iy2, is, icon, it, ifun, irun,
     *               isun, izun, igun, iwun, ia, ib, iave, epf,
     *               mhalf, mfull, isclp, nfmax, njmax,
     *               nvmz, nvmp, nvmn, iright, nv, isucs,
     *               isnew, ian, ibn, iabe, ivnxt, nvmax)
c
         if(isucs .eq. 1) then
            ndl=ndl+1
         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/1000)*1000)
     *      write(*,*)'Number of points deleted =',i
  720 continue
      ileft = isnew
      if(is(isnew).eq.0) stop 520
c
      return
      end
*VORCMP
c
c     This subroutine computes Voronoi vertices.
c 
c
      subroutine vorcmp(ix, iy, ix2, iy2, xp, yp, icon, it, ifl, ib,
     *                  ibn, dscle, mhalf, mfull, isclp)
c
      integer ix(*), iy(*), ix2(*), iy2(*)
      double precision xp(*), yp(*)
      integer icon(6,*), it(*), ifl(*), ib(*)
      integer ibn, 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, j, ifir, isec, ithi, isgo, iko
      integer isgnx, isgny, iknx, ikny
c
c     compute vertices
c
      ng = 0
      dfull = dble(mfull)
      r215 = dble(mhalf)
      deps = dble(0.9)
      if(dscle.lt.deps) stop 620
      dscl2 = 2.0d0*dscle
      write(*,*)' '
      do 1000 j = 1, ibn
         i = ib(j)
         if(j.le.(j/10000)*10000)
     *   write(*,*)'Number of computed vertices =',j
         if(it(i).eq.0) stop 630
         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 640
         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 660
         call doubnm(io, isgo, iko, r215, dnom)
         if(dnom.lt.deps) stop 680
         call doubnm(inx, isgnx, iknx, r215, xnum)
         call doubnm(iny, isgny, ikny, r215, ynum)
         xp(i) = (xnum/dnom)/dscl2
         yp(i) = (ynum/dnom)/dscl2
 1000 continue
c
c     identify vertices whose sizes are too big
c
      do 1600 j = 1, ibn
         i = ib(j)
         ifl(i) = 1
         dnom = dmax1(dabs(xp(i)),dabs(yp(i)))
         if(dnom.gt.dfull) ifl(i) = 0
 1600 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, it, iave,
     *                  ia, ib, epf, mhalf, mfull, isclp, nvmz, ileft,
     *                  iright, n, iabe, ian, ibn, njmax, isite, nvmax,
     *                  ivnxt, insrt, itype, ipoint)
c
      integer mhalf, mfull, isclp(*)
      integer nvmz, ileft, iright, n, iabe, ian, ibn, njmax, isite
      integer nvmax, ivnxt, insrt, itype, ipoint
      double precision epf
      double precision x(*), y(*)
      integer ix(*), iy(*), ix2(*), iy2(*)
      integer is(*), icon(6,*), it(*), iave(*), ia(*), ib(*)
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
c     WRITE(*,*)'ISITE=',ISITE,' ITYPE=',ITYPE
      if(itype .eq. -1 .or. itype .eq. 1) go to 400
      if(insrt .eq. 0) go to 400
      if(isite.le.0) stop 700
      if(itype .eq. 2) then
          call edgint(icon, it, is, iave, ia, ib, ivnxt, nvmax, nvmz,
     *                ivfnd, isite, iright, iabe, ian, ibn, njmax)
      elseif(itype .eq. 3) then
          call intins(icon, it, is, iave, ia, ib, ivnxt, nvmax, nvmz,
     *                ivfnd, isite, iabe, ian, ibn, njmax)
      else
          call outhul(x, y, ix, iy, ix2, iy2, icon, it, is, iave, ib,
     *                ivnxt, nvmax, nvmz, mhalf, mfull, isclp,
     *                ivfnd, isite, iright, iabe, ibn, njmax, 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, it, is, iave, ia, ib,
     *            epf, ivcur, isite, nvmz, mhalf, mfull, iabe, ian,
     *            ibn, ivnxt, njmax, nvmax, 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 800
      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 805
      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 810
      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 815
      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 820
      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 825
      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 830
      if(ivadj .ne. 0 .and. ivadj .ne. nvmz) go to 200
      if(ivlst .ne. 0 .and. ivlst .ne. nvmz) stop 835
      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 840
      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 845
      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 850
      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 855
      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 860
         if(itype .eq.  1 .and. ileft .eq. iright) stop 865
         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 870
      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 875
         endif
         if(ivlst .ne. 0 .and. ivlst .ne. nvmz) stop 880
         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 885
         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 890
         endif
         ivcur = iabs(ivadj)
         if(ivcur.eq.ivini) stop 895
         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 898
 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 900
      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 905
               else
                  if(dot3.gt.-epf) stop 910
               endif
            else
               if(dot3.gt.-epf) stop 915
            endif
            isi2 = isi1
 2300    continue
      else
         stop 920
      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, it, is, iave, ia, ib, ivnxt, nvmax, nvmz,
     *                  ivcur, isite, islt, iabe, ian, ibn, njmax)
c
      integer ivnxt, nvmax, nvmz, ivcur, isite, islt
      integer iabe, ian, ibn, njmax
      integer icon(6,*), it(*), is(*), iave(*), ia(*), ib(*)
c
      integer it1, it2, it3, it4, iv1, iv2, iv3, iv4
      integer isrt, isad, ivadj, isop, ivodj
c
      if(it(ivcur).eq.1) then
         it(ivcur) = -1
         ian = ian + 1
         if(ian .gt. njmax) stop 930
         ia(ian) = ivcur
         it1 = iave(iabe)
         if(iabe .eq. 1) then
            ivnxt = it1
            if(ivnxt .gt. nvmax) stop 935
            iave(1) = ivnxt + 1
         else
            iabe = iabe - 1
         endif
         if(it(it1).ne.0 .and. it(it1).ne.-2) stop 940
         if(it(it1).eq.0) then
            ibn = ibn+1
            if(ibn .gt. njmax) stop 945
            ib(ibn) = it1
         endif
         it(it1) = 2
      elseif(it(ivcur).eq.2) then
         it1 = ivcur
      else
         stop 950
      endif
c
      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 955
      endif
c
      ivodj = iabs(ivadj)
      it3 = ivodj
      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 960
      endif
c
  100 continue
      it2 = iave(iabe)
      if(iabe .eq. 1) then
          ivnxt = it2
          if(ivnxt .gt. nvmax) stop 965
          iave(1) = ivnxt + 1
      else
          iabe = iabe - 1
      endif
      if(it(it2).ne.0 .and. it(it2).ne.-2) stop 970
      if(it(it2).eq.0) then
         ibn = ibn+1
         if(ibn .gt. njmax) stop 975
         ib(ibn) = it2
      endif
      it(it2) = 2
      if(ivadj .eq. 0 .or. ivadj .eq. nvmz) go to 200
c
      it3 = iave(iabe)
      if(iabe .eq. 1) then
          ivnxt = it3
          if(ivnxt .gt. nvmax) stop 980
          iave(1) = ivnxt + 1
      else
          iabe = iabe - 1
      endif
      if(it(it3).ne.0 .and. it(it3).ne.-2) stop 985
      if(it(it3).eq.0) then
         ibn = ibn+1
         if(ibn .gt. njmax) stop 990
         ib(ibn) = it3
      endif
      it(it3) = 2
c
      if(it(it4).eq.1) then
         it(it4) = -1
         ian = ian + 1
         if(ian .gt. njmax) stop 995
         ia(ian) = it4
         it4 = iave(iabe)
         if(iabe .eq. 1) then
            ivnxt = it4
            if(ivnxt .gt. nvmax) stop 1000
            iave(1) = ivnxt + 1
         else
            iabe = iabe - 1
         endif
         if(it(it4).ne.0 .and. it(it4).ne.-2) stop 1005
         if(it(it4).eq.0) then
            ibn = ibn+1
            if(ibn .gt. njmax) stop 1010
            ib(ibn) = it4
         endif
         it(it4) = 2
      elseif(it(it4).ne.2) then
         stop 1015
      endif
c
  200 continue
      is(isite) = it1
      if(is(islt) .eq. ivcur) is(islt) = it1
      if(is(isad) .eq. ivcur) is(isad) = it2
      if(is(isrt) .eq. ivcur) is(isrt) = it2
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 1020
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
      if(is(isrt) .eq. ivodj) is(isrt) = it3
      if(is(islt) .eq. ivodj) is(islt) = it4
      if(is(isop) .eq. ivodj) is(isop) = it4
c
      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 1025
c
      if(ivadj .lt. 0) then
          icon(3,it3) = -it2
          icon(2,it4) = -it1
      endif
c
c     update neighboring triangles
c
  300 continue
      if(iv1 .eq. 0 .or. iv1 .eq. nvmz) go to 350
      ivodj = ivcur
      if(iv1 .lt. 0) then
          iv1 = -iv1
          it1 = -it1
          ivodj = -ivodj
      endif
c
      if(icon(1,iv1) .eq. ivodj) then
          icon(1,iv1) = it1
      elseif(icon(2,iv1) .eq. ivodj) then
          icon(2,iv1) = it1
      elseif(icon(3,iv1) .eq. ivodj) then
          icon(3,iv1) = it1
      else
          stop 1030
      endif
  350 continue
c
      if(iv2 .eq. 0 .or. iv2 .eq. nvmz) go to 400
      ivodj = ivcur
      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 1035
      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 450
      ivodj = iabs(ivadj)
      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 1040
      endif
  450 continue
c
      if(iv4 .eq. 0 .or. iv4 .eq. nvmz) go to 500
      ivodj = iabs(ivadj)
      if(iv4 .lt. 0) then
          iv4 = -iv4
          it4 = -it4
          ivodj = -ivodj
      endif
c
      if(icon(1,iv4) .eq. ivodj) then
          icon(1,iv4) = it4
      elseif(icon(2,iv4) .eq. ivodj) then
          icon(2,iv4) = it4
      elseif(icon(3,iv4) .eq. ivodj) then
          icon(3,iv4) = it4
      else
          stop 1045
      endif
  500 continue
c
      ivcur = iabs(it1)
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, it, is, iave, ia, ib, ivnxt, nvmax, nvmz,
     *                  ivcur, isite, iabe, ian, ibn, njmax)
c
      integer ivnxt, nvmax, nvmz, ivcur, isite, iabe, ian, ibn, njmax
      integer icon(6,*), it(*), is(*), iave(*), ia(*), ib(*)
c
      integer it1, it2, it3, iv1, iv2, iv3, is4, is5, is6, ivadj
c
      if(it(ivcur).eq.1) then
         it(ivcur) = -1
         ian = ian + 1
         if(ian .gt. njmax) stop 1050
         ia(ian) = ivcur
         it1 = iave(iabe)
         if(iabe .eq. 1) then
            ivnxt = it1
            if(ivnxt .gt. nvmax) stop 1055
            iave(1) = ivnxt + 1
         else
            iabe = iabe - 1
         endif
         if(it(it1).ne.0 .and. it(it1).ne.-2) stop 1060
         if(it(it1).eq.0) then
            ibn = ibn+1
            if(ibn .gt. njmax) stop 1065
            ib(ibn) = it1
         endif
         it(it1) = 2
      elseif(it(ivcur).eq.2) then
         it1 = ivcur
      else
         stop 1070
      endif
c
      it2 = iave(iabe)
      if(iabe .eq. 1) then
          ivnxt = it2
          if(ivnxt .gt. nvmax) stop 1075
          iave(1) = ivnxt + 1
      else
          iabe = iabe - 1
      endif
      if(it(it2).ne.0 .and. it(it2).ne.-2) stop 1080
      if(it(it2).eq.0) then
         ibn = ibn+1
         if(ibn .gt. njmax) stop 1085
         ib(ibn) = it2
      endif
      it(it2) = 2
c
      it3 = iave(iabe)
      if(iabe .eq. 1) then
          ivnxt = it3
          if(ivnxt .gt. nvmax) stop 1090
          iave(1) = ivnxt + 1
      else
          iabe = iabe - 1
      endif
      if(it(it3).ne.0 .and. it(it3).ne.-2) stop 1095
      if(it(it3).eq.0) then
         ibn = ibn+1
         if(ibn .gt. njmax) stop 1100
         ib(ibn) = it3
      endif
      it(it3) = 2
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 1105
c
      is(isite) = it1
      if(is(is4) .eq. ivcur) is(is4) = it3
      if(is(is5) .eq. ivcur) is(is5) = it1
      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(iv1 .eq. 0 .or. iv1 .eq. nvmz) go to 50
      ivadj = ivcur
      if(iv1 .lt. 0) then
          iv1 = -iv1
          it1 = -it1
          ivadj = -ivadj
      endif
c
      if(icon(1,iv1) .eq. ivadj) then
          icon(1,iv1) = it1
      elseif(icon(2,iv1) .eq. ivadj) then
          icon(2,iv1) = it1
      elseif(icon(3,iv1) .eq. ivadj) then
          icon(3,iv1) = it1
      else
          stop 1110
      endif
   50 continue
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 1115
      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 1120
      endif
  200 continue
c
      ivcur = iabs(it1)
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, it, is, iave, ib,
     *                  ivnxt, nvmax, nvmz, mhalf, mfull, isclp,
     *                  ivfnd, isite, islt, iabe, ibn, njmax, epf)
c
      integer ivnxt, nvmax, nvmz, mhalf, mfull, isclp(*)
      integer ivfnd, isite, islt, iabe, ibn, njmax
      double precision epf
      double precision x(*), y(*)
      integer ix(*), iy(*), ix2(*), iy2(*)
      integer icon(6,*), it(*), is(*), iave(*), ib(*)
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 1125
      endif
c
      if(ivadj .eq. 0 .or. ivadj .eq. nvmz) go to 200
      ivcur = iabs(ivadj)
      if(ivcur.eq.ivini) stop 1130
      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
c     WRITE(*,*)'1 ISITE=',ISITE,' ISRT=',ISRT,' ISLT=',ISLT,' DOT=',DOT
      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
      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 1135
      endif
      if(ivadj .ne. 0 .and. ivadj .ne. nvmz) stop 1140
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 1145
          iave(1) = ivnxt + 1
      else
          iabe = iabe - 1
      endif
      if(it(ivfnd).ne.0 .and. it(ivfnd).ne.-2) stop 1150
      if(it(ivfnd).eq.0) then
         ibn = ibn+1
         if(ibn .gt. njmax) stop 1155
         ib(ibn) = ivfnd
      endif
      it(ivfnd) = 2
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 1160
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 1170
          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 1180
      if(ivlst.ne.ivfnd) stop 1190
      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, it, is, iave, ia,
     *                 ib, epf, ivcur, isite, nvmz, mhalf, mfull, iabe,
     *                 ian, ibn, ivnxt, njmax, nvmax, isclp)
c
      integer ivcur, isite, nvmz, mhalf, mfull, isclp(*)
      double precision x(*), y(*)
      integer ix(*), iy(*), ix2(*), iy2(*)
      integer icon(6,*), it(*), is(*), iave(*), ia(*), ib(*)
      integer iabe, ian, ibn, ivnxt, njmax, nvmax
      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, it1, it2
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 1200
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 1205
      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 1210
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
      if(it(ivcur).eq.1) then
         it(ivcur) = -1
         ian = ian + 1
         if(ian .gt. njmax) stop 1215
         ia(ian) = ivcur
         it1 = iave(iabe)
         if(iabe .eq. 1) then
            ivnxt = it1
            if(ivnxt .gt. nvmax) stop 1220
            iave(1) = ivnxt + 1
         else
            iabe = iabe - 1
         endif
         if(it(it1).ne.0 .and. it(it1).ne.-2) stop 1225
         if(it(it1).eq.0) then
            ibn = ibn+1
            if(ibn .gt. njmax) stop 1230
            ib(ibn) = it1
         endif
         it(it1) = 2
      elseif(it(ivcur).eq.2) then
         it1 = ivcur
      else
         stop 1235
      endif
      if(it(ivadj).eq.1) then
         it(ivadj) = -1
         ian = ian + 1
         if(ian .gt. njmax) stop 1240
         ia(ian) = ivadj
         it2 = iave(iabe)
         if(iabe .eq. 1) then
            ivnxt = it2
            if(ivnxt .gt. nvmax) stop 1245
            iave(1) = ivnxt + 1
         else
            iabe = iabe - 1
         endif
         if(it(it2).ne.0 .and. it(it2).ne.-2) stop 1250
         if(it(it2).eq.0) then
            ibn = ibn+1
            if(ibn .gt. njmax) stop 1255
            ib(ibn) = it2
         endif
         it(it2) = 2
      elseif(it(ivadj).eq.2) then
         it2 = ivadj
      else
         stop 1260
      endif
      icon(1,it1) = icon(3,ivadj)
      icon(3,it2) = it1
      icon(2,it2) = icon(2,ivcur)
      icon(2,it1) = it2
      icon(3,it1) = icon(3,ivcur)
      icon(1,it2) = icon(1,ivadj)
      icon(4,it1) = isite
      icon(5,it1) = isadj
      icon(6,it1) = ishat
      icon(4,it2) = isite
      icon(5,it2) = ishat
      icon(6,it2) = iscur
      ivout = icon(1,it1)
      if(ivout .eq. 0 .or. ivout .eq. nvmz) go to 420
      ivcar = it1
      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 1265
      endif
  420 continue
      ivout = icon(1,it2)
      if(ivout .eq. 0 .or. ivout .eq. nvmz) go to 440
      ivcar = it2
      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 1270
      endif
  440 continue
      ivout = icon(3,it1)
      if(ivout .eq. 0 .or. ivout .eq. nvmz) go to 460
      ivodj = it1
      if(ivout .lt. 0) then
          ivout = -ivout
          ivodj = -ivodj
      endif
      icon(2,ivout) = ivodj
  460 continue
      ivout = icon(2,it2)
      if(ivout .eq. 0 .or. ivout .eq. nvmz) go to 480
      ivodj = it2
      if(ivout .lt. 0) then
          ivout = -ivout
          ivodj = -ivodj
      endif
      icon(3,ivout) = ivodj
  480 continue
c
      if(is(isite) .eq. ivcur) is(isite) = it1
      if(is(isadj) .eq. ivcur) is(isadj) = it1
      if(is(iscur) .eq. ivcur) is(iscur) = it2
      if(is(isadj) .eq. ivadj) is(isadj) = it1
      if(is(ishat) .eq. ivadj) is(ishat) = it2
      if(is(iscur) .eq. ivadj) is(iscur) = it2
c
      ivcur = it2
      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 1275
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
*EDGSTR
c
c     subroutine edgstr 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, it, ifun, igun,
     *                  iwun, iave, ib, epf, mhalf, mfull, isclp, lf,
     *                  njmax, iabe, ibn, j, iopt, ivnxt, nvmax)
c
      integer mhalf, mfull, isclp(*)
      integer lf, iabe, ibn, njmax, j, iopt, ivnxt, nvmax
      double precision epf
      double precision x(*), y(*)
      integer ix(*), iy(*), ix2(*), iy2(*)
      integer icon(6,*), it(*), iave(*), ib(*)
      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
      ial = iave(iabe)
      if(iabe .eq. 1) then
          ivnxt = ial
          if(ivnxt .gt. nvmax) stop 1300
          iave(1) = ivnxt + 1
      else
          iabe = iabe - 1
      endif
      if(it(ial).ne.0 .and. it(ial).ne.-2) stop 1310
      if(it(ial).eq.0) then
         ibn = ibn+1
         if(ibn .gt. njmax) stop 1320
         ib(ibn) = ial
      endif
      it(ial) = 2
      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 1600
          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 1610
          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 1620
          endif
          if(ivpre .ne. 0) stop 1630
          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 1640
          endif
          if(ivadj .ne. 0) go to 500
          if(isadj .ne. isnxt) stop 1650
          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 1660
          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, it, ifun,
     *                  irun, isun, izun, igun, iwun, ia, ib, iave,
     *                  epf, mhalf, mfull, isclp, nfmax, njmax, nvmz,
     *                  nvmp, nvmn, isite, n, isucs, isnew, ian, ibn,
     *                  iabe, ivnxt, nvmax)
c
      integer mhalf, mfull, isclp(*)
      integer nfmax, njmax, nvmz, nvmp, nvmn
      integer isite, n, isucs, isnew, ian, ibn, iabe, ivnxt, nvmax
      double precision epf
      double precision x(*), y(*)
      integer ix(*), iy(*), ix2(*), iy2(*)
      integer is(*), icon(6,*), it(*), ia(*), ib(*)
      integer ifun(*), irun(*), isun(*), izun(*)
      integer igun(*), iwun(*), iave(*)
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
      integer ivcur, ivtro, l, isam, 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 1800
c
c     obtain 'right' side information
c
      isucs = 1
      if(isite .lt. 1 .or. isite .gt. n) stop 1810
      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 1820
      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 1830
c
c     eliminate first triangle
c
      lr = 2
      if(lr .gt. nfmax) stop 1840
      irun(1) = ispre
      irun(2) = isadj
      izun(1) = ivpre
c     icon(4,ivini) = -1
      if(it(ivini).eq.1) then
         it(ivini) = -1
         ian = ian+1
         if(ian .gt. njmax) stop 1850
         ia(ian) = ivini
      elseif(it(ivini).eq.2) then
         it(ivini) = -2
         iabe = iabe + 1
         if(iabe .gt. njmax) stop 1860
         iave(iabe) = ivini
      else
         stop 1870
      endif
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 1880
      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 1885
      lr = lr + 1
      if(lr .gt. nfmax) stop 1890
      irun(lr) = isadj
      izun(lr-1) = ivpre
c     icon(4,ivcur) = -1
      if(it(ivcur).eq.1) then
         it(ivcur) = -1
         ian = ian + 1
         if(ian .gt. njmax) stop 1900
         ia(ian) = ivcur
      elseif(it(ivcur).eq.2) then
         it(ivcur) = -2
         iabe = iabe + 1
         if(iabe .gt. njmax) stop 1910
         iave(iabe) = ivcur
      else
         stop 1915
      endif
      go to 200
c
c     obtain 'left' side information
c
  400 continue
      if(ivout.eq.0) stop 1920
      lf = lf + 1
      if(lf .gt. nfmax) stop 1930
      ifun(lf) = isadj
      isun(lf) = ivpre
c     icon(4,ivcur) = -1
      if(it(ivcur).eq.1) then
         it(ivcur) = -1
         ian = ian + 1
         if(ian .gt. njmax) stop 1940
         ia(ian) = ivcur
      elseif(it(ivcur).eq.2) then
         it(ivcur) = -2
         iabe = iabe + 1
         if(iabe .gt. njmax) stop 1950
         iave(iabe) = ivcur
      else
         stop 1955
      endif
      if(isadj .eq. ispre) go to 500
      go to 200
  500 continue
      if(lf .lt. 2) stop 1960
      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
c         if(icon(4,iabs(isam)) .ne. -1) go to 900
          if(it(iabs(isam)).ne.-1 .and. it(iabs(isam)).ne.-2) 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
c         if(icon(4,iabs(isam)) .ne. -1) go to 800
          if(it(iabs(isam)).ne.-1 .and. it(iabs(isam)).ne.-2) go to 800
          if(isam .gt. 0) then
              izun(l) = nvmp
          else
              izun(l) = nvmn
          endif
  800 continue
c
      isnew = irun(lr)
      if(ivout .eq. 0) go to 1000
      if(lf.eq.0) stop 1970
      go to 1200
c
c     point to be eliminated is on convex hull
c
 1000 continue
      lr = lr + 1
      if(lr .gt. nfmax) stop 1980
      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, it, irun, igun, iwun,
     *            iave, ib, epf, mhalf, mfull, isclp, lr, njmax,
     *            iabe, ibn, jr, iopt, ivnxt, nvmax)
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(*,*)'After deletion of a point other points are colinear'
         write(*,*)'Program does not handle this unlikely event'
         write(*,*)'Program terminated'
         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 1990
      endif
      icon(indx,ivabs) = 0
      if(ivcur .lt. 0) icon(indx,ivabs) = nvmz
      is(iscor) = ivabs
c     icon(4,ial) = -1
      if(it(ial).eq.2) then
         it(ial) = -2
         iabe = iabe + 1
         if(iabe .gt. njmax) stop 2000
         iave(iabe) = ial
      else
         stop 2010
      endif
      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, it, irun, igun, iwun,
     *            iave, ib, epf, mhalf, mfull, isclp, lr, njmax,
     *            iabe, ibn, jr, iopt, ivnxt, nvmax)
c
c     update 'left' side information
c
      lg = lf
      do 1300 j = 2, jr
          lf = lf + 1
          if(lf .gt. nfmax) stop 2020
          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, it, ifun, igun, iwun,
     *            iave, ib, epf, mhalf, mfull, isclp, lf, njmax,
     *            iabe, ibn, jf, iopt, ivnxt, nvmax)
c
      if(jf .ne. 2) stop 2030
      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 2040
      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 2050
      endif
c
 2000 continue
      is(isnew) = ial
      do 2100 l = 2, lr
          lf = lf + 1
          if(lf .gt. nfmax) stop 2060
          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 point
c
 2200 continue
      is(isite) = 0
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 3000
      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
c*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
c     subroutine comprs(is, icon, id, it, n, ivnxt, nvmax, nvmz)
c
c     integer n, ivnxt, nvmax, nvmz
c     integer is(*), icon(6,*), id(*), it(*)
c     integer ivnot, i, j, ivcur
c
c     test n, ivnxt, nvmz
c
c     if(n .gt. nvmax .or. ivnxt .gt. nvmax) stop 3100
c
c     if(nvmz .ne. -nvmax - 1) stop 3120
c
c     compress icon
c
c     ivnot = 0
c     do 200 i = 1, ivnxt
c         if(icon(4,i) .eq. -1) go to 200
c         if(it(i).le.0) go to 200
c         ivnot = ivnot + 1
c         id(i) = ivnot
c         do 100 j = 1, 6
c             icon(j,ivnot) = icon(j,i)
c 100     continue
c 200 continue
c     ivnxt = ivnot
c
c     update icon for triangles
c
c     do 400 i = 1, ivnxt
c         do 300 j = 1, 3
c             ivcur = icon(j,i)
c             if(ivcur .eq. 0 .or. ivcur .eq. nvmz) go to 300
c             if(ivcur .gt. 0) then
c                 icon(j,i) = id(ivcur)
c             else
c                 ivcur = -ivcur
c                 icon(j,i) = -id(ivcur)
c             endif
c 300     continue
c 400 continue
c
c     update is for triangles
c
c     do 500 i = 1, n
c         if(is(i) .le. 0) go to 500
c         is(i) = id(is(i))
c 500 continue
c
c     return
c     end
*CONTST
c
c     subroutine contst to -
c
c     test consistency of triangulation
c
c     November 9, 1993
c
      subroutine contst(icon, it, is, id, nv, ivnxt, icnxt, nvmz)
c
      integer nv, ivnxt, icnxt, nvmz
      integer icon(6,*), it(*), 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. icnxt) stop 3410
c
      do 1000 i = 1, ivnxt
          if(id(i).ne.0 .and. it(i).ne.0) stop 3420
 1000 continue
c
c     write(*,*)'****************************************'
c     write(*,*)'consistency check satisfied'
c     write(*,*)'****************************************'
c
      return
      end
*CNSTST
c
c     subroutine contst to -
c
c     test consistency of triangulation
c
c     November 9, 1993
c
      subroutine cnstst(icon, it, is, ib, iz, ibn, izn, nvmz)
c
      integer nvmz, ibn, izn
      integer icon(6,*), it(*), is(*), ib(*), iz(*)
      integer site1, site2, i, iscur
      integer ivcur, ivcur2, ivabs2, ivlst, isit1, isit2
c
c     test initial triangle for each site
c
      do 50 i = 1, izn
          site1 = iz(i)
          iscur = is(site1)
          if(iscur .le. 0) stop 3440
          if(icon(4,iscur).ne.site1 .and. icon(5,iscur).ne.site1 .and.
     *       icon(6,iscur) .ne. site1) stop 3450
   50 continue
c
c     test each triangle for neighbor connections
c
      do 300 i = 1, ibn
         ivcur = ib(i)
         if(it(ivcur).eq.0) stop 3460
c
         ivcur2 = icon(1,ivcur)
         if(ivcur2 .eq. 0 .or. ivcur2 .eq. nvmz) go to 100
         ivabs2 = iabs(ivcur2)
         ivlst = ivcur
         if(ivcur2 .lt. 0) ivlst = -ivcur
         site1 = icon(5,ivcur)
         site2 = icon(6,ivcur)
         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) then
            isit1 = icon(6,ivabs2)
            isit2 = icon(5,ivabs2)
         else
            stop 3470
         endif
         if(site1 .ne. isit1 .or. site2 .ne. isit2) stop 3480
c
  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(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) then
            isit1 = icon(6,ivabs2)
            isit2 = icon(5,ivabs2)
         else
            stop 3490
         endif
         if(site1 .ne. isit1 .or. site2 .ne. isit2) stop 3500
c
  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(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) then
            isit1 = icon(6,ivabs2)
            isit2 = icon(5,ivabs2)
         else
            stop 3510
         endif
         if(site1 .ne. isit1 .or. site2 .ne. isit2) stop 3520
c
  300 continue
c
c     write(*,*)'****************************************'
c     write(*,*)'consistency check satisfied'
c     write(*,*)'****************************************'
c
      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, it, ib, ibn,
     *                  mhalf, mfull, isclp, epf)
c
      integer mhalf, mfull, isclp(*), ibn
      double precision epf
      double precision x(*), y(*)
      integer ix(*), iy(*), ix2(*), iy2(*)
      integer icon(6,*), it(*), ib(*)
      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, l
      integer icerr, itide
c
c     initialize
c
      icerr = 0
c
c     test all triangles
c
      do 200 l = 1, ibn
         i = ib(l)
         if(it(i).eq.0) stop 3540
         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
      if(icerr.ne.0) stop 3570
c     write(*,*) '************************************************'
c     write(*,*) 'Number of circle criterion violations = ',icerr
c     write(*,*) '(This number should equal 0)'
c     write(*,*) '************************************************'
c     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, it, ib, ibn,
     *                  mhalf, mfull, isclp, epf)
c
      double precision x(*), y(*)
      integer ix(*), iy(*), ix2(*), iy2(*)
      integer icon(6,*), it(*), ib(*)
      double precision epf
      integer isclp(*), mhalf, mfull, ibn
      integer a, b, c, i, j, isga, idmin
      double precision delx, dely, dif, xtemp, ytemp, dlst, dot
c
c     test all triangles
c
      idmin = 0
      do 300 j=1,ibn
         i = ib(j)
         if(it(i).eq.0) stop 3590
         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
      if(idmin.ne.0) stop 3620
c
c     write(*,*) '************************************************'
c     write(*,*) 'Number of orientation violations = ',idmin
c     write(*,*) '(This number should equal 0)'
c     write(*,*) '************************************************'
c     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, it, is, nv, nt,
     *                  nvmz, mhalf, mfull, isclp, epf)
c
      double precision x(*), y(*)
      integer ix(*), iy(*), ix2(*), iy2(*)
      integer icon(6,*), it(*), 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(it(i).eq.0) stop 3680
         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
      if(idmin.ne.0) stop 3790
c     write(*,*) '************************************************'
c     write(*,*) 'Number of convexity violations = ',idmin
c     write(*,*) '(This number should equal 0)'
c     write(*,*) '************************************************'
c     write(*,*) ''
c
      return
      end
*CNVTST
c
c  Subroutine to verify convexity of union of triangles
c
      subroutine cnvtst(x, y, ix, iy, ix2, iy2, icon, is, iz,
     *                  izn, nvmz, mhalf, mfull, isclp, epf)
c
      double precision x(*), y(*)
      integer ix(*), iy(*), ix2(*), iy2(*)
      integer icon(6,*), is(*), iz(*), izn, nvmz
      double precision epf
      integer isclp(*), mhalf, mfull, isga, iflag, idmin
      integer i, ivcur, ivini, isite, iscur, isadj, ivadj, ivnxt
      double precision delx, dely, dif, xtemp, ytemp, dlst, dot
c
c     test current site
c
      idmin = 0
      do 1000 i = 1, izn
         isite = iz(i)
         ivini = is(isite)
         iflag = 0
c
c        test for boundary of convex hull
c
         if(icon(4,ivini).eq.isite) then
            iscur = icon(6,ivini)
            ivadj = icon(3,ivini)
            isadj = icon(5,ivini)
            ivnxt = icon(2,ivini)
         elseif(icon(5,ivini).eq.isite) then
            iscur = icon(4,ivini)
            ivadj = icon(1,ivini)
            isadj = icon(6,ivini)
            ivnxt = icon(3,ivini)
         elseif(icon(6,ivini).eq.isite) then
            iscur = icon(5,ivini)
            ivadj = icon(2,ivini)
            isadj = icon(4,ivini)
            ivnxt = icon(1,ivini)
         else
            stop 3800
         endif
         if(ivadj.eq.0 .or. ivadj.eq.nvmz) iflag = 1
c
c        go around site
c
  100    continue
         if(ivnxt.eq.0 .or. ivnxt.eq.nvmz) then
            if(iflag.ne.1) stop 3810
            go to 500
         endif
         ivcur = iabs(ivnxt)
         if(icon(4,ivcur).eq.isite) then
            iscur = icon(6,ivcur)
            ivnxt = icon(2,ivcur)
         elseif(icon(5,ivcur).eq.isite) then
            iscur = icon(4,ivcur)
            ivnxt = icon(3,ivcur)
         elseif(icon(6,ivcur).eq.isite) then
            iscur = icon(5,ivcur)
            ivnxt = icon(1,ivcur)
         else
            stop 3820
         endif
         if(iabs(ivnxt).eq.ivini) go to 1000
         go to 100
c
c        test current vertex isite at boundary of convex hull
c
  500    continue
         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
 1000 continue
c
      if(idmin.ne.0) stop 3830
c     write(*,*) '************************************************'
c     write(*,*) 'Number of convexity violations = ',idmin
c     write(*,*) '(This number should equal 0)'
c     write(*,*) '************************************************'
c     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