c23456789 123456789 123456789 123456789 123456789 123456789 123456789 12
*MAIN
c
c     main program to -
c
c     produce data files and script file for plotting in the form of a
c     movie Voronoi diagrams of dynamic data, i. e. of a set of moving
c     points in 2-d space.
c
c     It is assumed that another existing program, program dynvor, has
c     been executed thus producing the Delaunay triangulations and
c     Voronoi diagrams of the dynamic data and the corresponding output
c     files. From the output files of program dynvor the program here
c     produces data files and script file vorplot.m that when executed
c     creates from the data files a movie of the plots of the Voronoi
c     diagrams of the dynamic data in the order in which program dynvor
c     reads the input data files that define the dynamic data.
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     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     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, 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     njmax  Parameter which is the dimension of arrays 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     x1     Array that contains x-coordinates of points or vertices
c            1 through nv for previous data set of input points.
c            Each line in input file 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     y1     Array that contains y-coordinates of points or vertices
c            1 through nv for previous data set of input points.
c            (see description of array x above). nmax is the dimension
c            for this array.
c
c     x2     Array similar to array x1 (described above) for current
c            data set of input points.
c
c     y2     Array similar to array y1 (described above) for current
c            data set of input points.
c
c     vx     Array that contains x-coordinates of circumcenters
c            (Voronoi vertices) of triangles 1 through ivnxt.
c            nvmax is the dimension for this array.
c
c     vy     Array that contains y-coordinates of circumcenters
c            (Voronoi vertices) of triangles 1 through ivnxt.
c            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     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.
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.
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.
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
c***********************************************************************
c
      program main
c
      integer nmax, nvmax, nemax, nhmax, ndmax
      integer njmax, nkmax, nlmax
      parameter (nmax = 160000, nvmax = 2.1*nmax)
      parameter (nemax = 6*nmax, nhmax = 1000, ndmax = 40000)
      parameter (njmax = 3*ndmax, nkmax = ndmax, nlmax = 2*ndmax)
c
      double precision x1(nmax), y1(nmax), x2(nmax), y2(nmax)
      double precision vx(nvmax), vy(nvmax)
      integer is(nmax), icon(6,nvmax), id(nvmax)
      integer ied1(nvmax), ied2(nemax), iedn(nemax), ih(nhmax)
      integer idl(ndmax), ia(njmax), ib(njmax), ic(nkmax), iz(nlmax)
c
      integer idyn, idcr, n, nv, icnxt, ivnxt, nvmz, nd, iflag
      integer ien, iem, ihn, ian, ibn, icn, izn, icsfig
c
      character*4 az
      character*3 nu(100)
      character*10 pn(100), tn(100), vn(100)
      character*10 pe(100), ve(100), pa(100), va(100)
      character*7  fpd(100), fvd(100), fpc(100), fvc(100)
      character*5  spd(100), svd(100), spc(100), svc(100)
      character*5  spe(100), sve(100), spa(100), sva(100)
      character*6  fas
      character*2  lline
c
      integer i, j, ng, ivini, ivcur, ivadj, ivpre, iscur
      integer iecur, ienxt
      double precision xcor, ycor, xmax, xmin, ymax, ymin
      double precision vxmin, vxmax, vymin, vymax
c
      integer larg, izero
c
      izero = 0
c
      write(*,*)' '
      write(*,*)'This program is for preparing data for plotting ',
     *          'Voronoi'
      write(*,*)'diagrams of dynamic data.'
c
c     initialize for Voronoi edges data structure
c
      ien = 0
      ihn = 0
      do 4 i = 1, nvmax
         ied1(i) = 0
    4 continue
      do 6 i = 1, nemax
         ied2(i) = 0
         iedn(i) = 0
    6 continue
c
c     initialize for names of dynamic input data files and output files
c
      az ='az.m'
c
      do 10 j = 1, 100
        do 8 i =1,3
           nu(j)(i:i) = ' '
    8   continue
   10 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 12 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'
   12 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) = ' '
            pe(j)(i:i) = ' '
            ve(j)(i:i) = ' '
            pa(j)(i:i) = ' '
            va(j)(i:i) = ' '
   20    continue
c
         pn(j)(1:2) = 'pn'
         tn(j)(1:2) = 'tn'
         vn(j)(1:2) = 'vn'
         pe(j)(1:2) = 'pe'
         ve(j)(1:2) = 've'
         pa(j)(1:2) = 'pa'
         va(j)(1:2) = 'va'
         fpd(j)(1:2) = 'pd'
         fvd(j)(1:2) = 'vd'
         fpc(j)(1:2) = 'pc'
         fvc(j)(1:2) = 'vc'
         spd(j)(1:2) = 'pd'
         svd(j)(1:2) = 'vd'
         spc(j)(1:2) = 'pc'
         svc(j)(1:2) = 'vc'
         spe(j)(1:2) = 'pe'
         sve(j)(1:2) = 've'
         spa(j)(1:2) = 'pa'
         sva(j)(1:2) = 'va'
c
         pn(j)(3:5) = nu(j)
         tn(j)(3:5) = nu(j)
         vn(j)(3:5) = nu(j)
         pe(j)(3:5) = nu(j)
         ve(j)(3:5) = nu(j)
         pa(j)(3:5) = nu(j)
         va(j)(3:5) = nu(j)
         fpd(j)(3:5) = nu(j)
         fvd(j)(3:5) = nu(j)
         fpc(j)(3:5) = nu(j)
         fvc(j)(3:5) = nu(j)
         spd(j)(3:5) = nu(j)
         svd(j)(3:5) = nu(j)
         spc(j)(3:5) = nu(j)
         svc(j)(3:5) = nu(j)
         spe(j)(3:5) = nu(j)
         sve(j)(3:5) = nu(j)
         spa(j)(3:5) = nu(j)
         sva(j)(3:5) = nu(j)
c
         if(j.lt.10) then
            pe(j)(4:5) = '.m'
            ve(j)(4:5) = '.m'
            pa(j)(4:5) = '.m'
            va(j)(4:5) = '.m'
         elseif(j.lt.100) then
            pe(j)(5:6) = '.m'
            ve(j)(5:6) = '.m'
            pa(j)(5:6) = '.m'
            va(j)(5:6) = '.m'
         else
            pe(j)(6:7) = '.m'
            ve(j)(6:7) = '.m'
            pa(j)(6:7) = '.m'
            va(j)(6:7) = '.m'
         endif
         fpd(j)(6:7) = '=['
         fvd(j)(6:7) = '=['
         fpc(j)(6:7) = '=['
         fvc(j)(6:7) = '=['
   25 continue
c
      lline = '];'
c
      fas = 'as=[  '
c
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
   30 continue
      read(11,*,end=40) xcor, ycor
      nv = nv + 1
      if(nv .gt. nmax) stop 10
      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 30
   40 continue
      close(unit=11)
      if(nv.lt.3) stop 20
      write(*,*)' '
      write(*,*)'Number of points read =',nv

      write(*,*)' '
      write(*,*)' xmin = ',xmin,' xmax = ',xmax
      write(*,*)' ymin = ',ymin,' ymax = ',ymax
      vxmin = xmin
      vxmax = xmax
      vymin = ymin
      vymax = ymax
      go to 110
c
c******************************************************************
c
c     create file for unplotting input points
c
   50 continue
      open(unit=20,file=pe(idcr))
      nv = n
      write(20,290) fpd(idcr), nv
      do 60 iscur = 1, nv
         write(20,*) iscur
   60 continue
      write(20,505) nv, lline
      close(unit=20)
c
c     create file for unplotting Voronoi cells of input points
c
      open(unit=22,file=ve(idcr))
      nv = n
      write(22,290) fvd(idcr), nv
      do 90 iscur = 1, nv
         ivini = is(iscur)
         ivcur = ivini
         ng = 0
   70    continue
         ng = ng+1
         if(icon(4,ivcur).eq.iscur) then
            ivadj = icon(2,ivcur)
         elseif(icon(5,ivcur).eq.iscur) then
            ivadj = icon(3,ivcur)
         elseif(icon(6,ivcur).eq.iscur) then
            ivadj = icon(1,ivcur)
         else
            stop 50
         endif
         ivcur = ivadj
         if(ivcur.ne.0 .and. ivcur.ne.ivini) go to 70
         if(ivcur.eq.ivini) ng = ng+1
         write(22,*) ng
         ivcur = ivini
         write(22,*) izero
   80    continue
         if(icon(4,ivcur).eq.iscur) then
            ivadj = icon(2,ivcur)
         elseif(icon(5,ivcur).eq.iscur) then
            ivadj = icon(3,ivcur)
         elseif(icon(6,ivcur).eq.iscur) then
            ivadj = icon(1,ivcur)
         else
            stop 60
         endif
         if(ivadj.eq.0) go to 90
         ivpre = ivcur
         ivcur = ivadj
c
         iecur = ied1(ivpre)
         if(iecur.le.0) stop 70
         if(ied2(iecur).eq.ivcur) then
            ienxt = iedn(iecur)
            if(ienxt.eq.0) then
               ied1(ivpre) = 0
               ied2(iecur) = 0
            else
               if(ienxt.lt.0) stop 100
               ied1(ivpre) = ienxt
               ied2(iecur) = 0
               iedn(iecur) = 0
            endif
         else
   82       continue
            ienxt = iedn(iecur)
            if(ienxt.le.0) stop 110
            if(ied2(ienxt).eq.ivcur) go to 84
            iecur = ienxt
            go to 82
   84       continue
            iedn(iecur) = iedn(ienxt)
            iecur = ienxt
            ied2(iecur) = 0
            iedn(iecur) = 0
         endif
c
         write(22,*) iecur
         if(ivcur.ne.ivini) go to 80
   90 continue
      write(22,505) nv, lline
      close(unit=22)
      ihn = 0
c
      do 95 i = 1, n
         x1(i) = x2(i)
         y1(i) = y2(i)
   95 continue
      go to 120
c
c*****************************************************************
c
c     read information about current Delaunay triangulation
c
c     write(*,*)' '
c     write(*,*)'Reading information about current Delaunay ',
c    *          'triangulation ... '
c
  110 continue
      open(unit=16,file=tn(idcr))
      read(16,*) iflag
  120 continue
      read(16,*) n, ivnxt, nvmz, icsfig, idyn
      if(n.ne.nv) stop 150
      if(nvmz.ne.-nvmax-1) stop 160
      if(icsfig.lt.0) stop 170
      read(16,130)((icon(i,j),i=1,6),j=1,ivnxt)
      read(16,140)(is(i),i=1,n)
  130 format(6i10)
  140 format(5i10)
      close(unit=16)
c
c     read information about Voronoi vertices
c
c     write(*,*)' '
c     write(*,*)'Reading circumcenters/Voronoi vertices ... '
c
      open(unit=18,file=vn(idcr))
      read(18,*) n, ivnxt, icsfig, larg
      if(n.ne.nv) stop 180
      if(larg.ne.0) then
         write(*,*)'At least one Voronoi vertex is too large'
         write(*,*)'Index of current file = ',idcr
         write(*,*)'Program terminated'
         stop
      endif
      read(18,190) (id(i),i=1,ivnxt)
      do 180 i=1,ivnxt
         read(18,*) vx(i),vy(i)
  180 continue
  190 format(40(1x,i1))
      close(unit=18)
c
c*****************************************************************
c
c     create file for plotting input points
c
      open(unit=20,file=pa(idcr))
      nv = n
      write(20,290) fpc(idcr), nv, nv, nv
      do 250 iscur = 1, nv
         write(20,*) x1(iscur), y1(iscur), iscur
  250 continue
      write(20,515) nv, nv, nv, lline
      close(unit=20)
c
c     create file for plotting Voronoi cells of input points
c
      open(unit=22,file=va(idcr))
  290 format(a7,i7,i7,i7)
      nv = n
      write(22,290) fvc(idcr), nv, nv, nv
      do 500 iscur = 1, nv
         ivini = is(iscur)
         ivcur = ivini
         ng = 0
  300    continue
         ng = ng+1
         if(icon(4,ivcur).eq.iscur) then
            ivadj = icon(2,ivcur)
         elseif(icon(5,ivcur).eq.iscur) then
            ivadj = icon(3,ivcur)
         elseif(icon(6,ivcur).eq.iscur) then
            ivadj = icon(1,ivcur)
         else
            stop 200
         endif
         ivcur = ivadj
         if(ivcur.ne.0 .and. ivcur.ne.ivini) go to 300
         if(ivcur.eq.ivini) ng = ng+1
         write(22,*) ng, ng, ng
         ivcur = ivini
         write(22,*) vx(ivcur), vy(ivcur), izero
  400    continue
         if(icon(4,ivcur).eq.iscur) then
            ivadj = icon(2,ivcur)
         elseif(icon(5,ivcur).eq.iscur) then
            ivadj = icon(3,ivcur)
         elseif(icon(6,ivcur).eq.iscur) then
            ivadj = icon(1,ivcur)
         else
            stop 220
         endif
         if(ivadj.eq.0) go to 500
         ivpre = ivcur
         ivcur = ivadj
c
         ien = ien+1
         if(ien.gt.nemax) stop 240
         if(ied1(ivpre).eq.0) then
            ied1(ivpre) = ien
            ied2(ien) = ivcur
            iedn(ien) = 0
         else
            iecur = ied1(ivpre)
  410       continue
            if(ied2(iecur).le.0) stop 260
            ienxt = iedn(iecur)
            if(ienxt.lt.0) stop 280
            if(ienxt.eq.0) go to 420
            iecur = ienxt
            go to 410
  420       continue
            iedn(iecur) = ien
            ied2(ien) = ivcur
            iedn(ien) = 0
         endif
c
         write(22,*) vx(ivcur), vy(ivcur), ien
         if(ivcur.ne.ivini) go to 400
  500 continue
      write(22,515) nv, nv, nv, lline
      close(unit=22)
      ihn = 0
c
  505 format(i7,a2)
  510 format(i7,i7,a2)
  515 format(i7,i7,i7,a2)
c
c*****************************************************************
c
c     read next points file
c
  595 continue
      if(idcr.ge.idyn) go to 1160
      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
  600 continue
      read(11,*,end=610) xcor, ycor
      nv = nv + 1
      if(nv .gt. nmax) stop 300
      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 600
  610 continue
      close(unit=11)
      if(n.ne.nv) stop 320

      write(*,*)' '
      write(*,*)' xmin = ',xmin,' xmax = ',xmax
      write(*,*)' ymin = ',ymin,' ymax = ',ymax
      if(xmin.lt.vxmin) vxmin = xmin
      if(xmax.gt.vxmax) vxmax = xmax
      if(ymin.lt.vymin) vymin = ymin
      if(ymax.gt.vymax) vymax = ymax
c
c**********************************************************************
c
c     read information about current Delaunay triangulation
c
c     write(*,*)' '
c     write(*,*)'Reading partial information about current Delaunay ',
c    *          'triangulation ... '
c
      open(unit=16,file=tn(idcr))
      read(16,*) iflag
      if(iflag.eq.0) go to 50
      read(16,*) n, ivnxt, nvmz, icsfig
      if(n.ne.nv) stop 340
      read(16,*) nd, izn, ian, ibn, icn, icnxt
      if(icnxt.gt.ivnxt) stop 360
      read(16,140)(idl(i),i=1,nd)
      read(16,140)(iz(i),i=1,izn)
c
c*********************************************************************
c
c     create file for unplotting points that have moved using
c     arrays idl
c
      open(unit=20,file=pe(idcr))
      write(20,290) fpd(idcr), nd
      do 660 i = 1, nd
         iscur = idl(i)
         write(20,*) iscur
  660 continue
      write(20,505) nd, lline
      close(unit=20)
c
c     create file for unplotting Voronoi cells of points affected
c     by moving points using arrays iz, is, icon
c
      open(unit=22,file=ve(idcr))
      write(22,290) fvd(idcr), izn
      do 690 i = 1, izn
         iscur = iz(i)
         ivini = is(iscur)
         ivcur = ivini
         ng = 0
  670    continue
         ng = ng+1
         if(icon(4,ivcur).eq.iscur) then
            ivadj = icon(2,ivcur)
         elseif(icon(5,ivcur).eq.iscur) then
            ivadj = icon(3,ivcur)
         elseif(icon(6,ivcur).eq.iscur) then
            ivadj = icon(1,ivcur)
         else
            stop 380
         endif
         ivcur = ivadj
         if(ivcur.ne.0 .and. ivcur.ne.ivini) go to 670
         if(ivcur.eq.ivini) ng = ng+1
         write(22,*) ng
         ivcur = ivini
         write(22,*) izero
  680    continue
         if(icon(4,ivcur).eq.iscur) then
            ivadj = icon(2,ivcur)
         elseif(icon(5,ivcur).eq.iscur) then
            ivadj = icon(3,ivcur)
         elseif(icon(6,ivcur).eq.iscur) then
            ivadj = icon(1,ivcur)
         else
            stop 400
         endif
         if(ivadj.eq.0) go to 690
         ivpre = ivcur
         ivcur = ivadj
c
         iecur = ied1(ivpre)
         if(iecur.le.0) stop  410
         if(ied2(iecur).eq.ivcur) then
            ienxt = iedn(iecur)
            if(ienxt.eq.0) then
               ied1(ivpre) = 0
               ied2(iecur) = 0
            else
               if(ienxt.lt.0) stop 420
               ied1(ivpre) = ienxt
               ied2(iecur) = 0
               iedn(iecur) = 0
            endif
         else
  682       continue
            ienxt = iedn(iecur)
            if(ienxt.le.0) stop 424
            if(ied2(ienxt).eq.ivcur) go to 684
            iecur = ienxt
            go to 682
  684       continue
            iedn(iecur) = iedn(ienxt)
            iecur = ienxt
            ied2(iecur) = 0
            iedn(iecur) = 0
         endif
         ihn = ihn+1
         if(ihn.gt.nhmax) stop 426
         ih(ihn) = iecur
c
         write(22,*) iecur
         if(ivcur.ne.ivini) go to 680
  690 continue
      write(22,505) izn, lline
      close(unit=22)
c
c
c*********************************************************************
c     
      read(16,140)(ia(i),i=1,ian)
      read(16,140)(ib(i),i=1,ibn)
      if(icn.ne.0) read(16,140)(ic(i),i=1,icn)
      read(16,130)((icon(i,ib(j)),i=1,6),j=1,ibn)
      if(icn.ne.0) read(16,130)((icon(i,ic(j)),i=1,6),j=1,icn)
      read(16,140)(is(iz(i)),i=1,izn)
      close(unit=16)
c
c     read information about Voronoi vertices
c
c     write(*,*)' '
c     write(*,*)'Reading circumcenters/Voronoi vertices ... '
c
      open(unit=18,file=vn(idcr))
      read(18,*) n, ivnxt, icsfig, larg
      if(n.ne.nv) stop 430
      read(18,1050) (id(ib(i)),i=1,ibn)
      do 1040 i=1,ibn
         j = ib(i)
         read(18,*) vx(j),vy(j)
 1040 continue
 1050 format(40(1x,i1))
      close(unit=18)
      write(*,*)' vxmin = ',vxmin,' vxmax = ',vxmax
      write(*,*)' vymin = ',vymin,' vymax = ',vymax
c
      if(ia(1).eq.izero .or. id(1).eq.izero) izero = 0
c
c*********************************************************************
c
c     create file for plotting points that have moved using
c     arrays idl, x2, y2
c
      open(unit=20,file=pa(idcr))
      write(20,290) fpc(idcr), nd, nd, nd
      do 1060 i = 1, nd
         iscur = idl(i)
         write(20,*) x2(iscur), y2(iscur), iscur
 1060 continue
      write(20,515) nd, nd, nd, lline
      close(unit=20)
c
c
c     create file for plotting Voronoi cells of points affected
c     by moving points using arrays iz, is, icon
c
      open(unit=22,file=va(idcr))
      write(22,290) fvc(idcr), izn, izn, izn
      do 1100 i = 1, izn
         iscur = iz(i)
         ivini = is(iscur)
         ivcur = ivini
         ng = 0
 1070    continue
         ng = ng+1
         if(icon(4,ivcur).eq.iscur) then
            ivadj = icon(2,ivcur)
         elseif(icon(5,ivcur).eq.iscur) then
            ivadj = icon(3,ivcur)
         elseif(icon(6,ivcur).eq.iscur) then
            ivadj = icon(1,ivcur)
         else
            stop 440
         endif
         ivcur = ivadj
         if(ivcur.ne.0 .and. ivcur.ne.ivini) go to 1070
         if(ivcur.eq.ivini) ng = ng+1
         write(22,*) ng, ng, ng
         ivcur = ivini
         write(22,*) vx(ivcur), vy(ivcur), izero
 1080    continue
         if(icon(4,ivcur).eq.iscur) then
            ivadj = icon(2,ivcur)
         elseif(icon(5,ivcur).eq.iscur) then
            ivadj = icon(3,ivcur)
         elseif(icon(6,ivcur).eq.iscur) then
            ivadj = icon(1,ivcur)
         else
            stop 460
         endif
         if(ivadj.eq.0) go to 1100
         ivpre = ivcur
         ivcur = ivadj
c
         if(ihn.eq.0) then
            ien = ien+1
            if(ien.gt.nemax) stop 480
            iem = ien
         else
            iem = ih(ihn)
            ihn = ihn-1
         endif
         if(ied1(ivpre).eq.0) then
            ied1(ivpre) = iem
            ied2(iem) = ivcur
            iedn(iem) = 0
         else
            iecur = ied1(ivpre)
 1090       continue
            if(ied2(iecur).le.0) stop 500
            ienxt = iedn(iecur)
            if(ienxt.lt.0) stop 520
            if(ienxt.eq.0) go to 1095
            iecur = ienxt
            go to 1090
 1095       continue
            iedn(iecur) = iem
            ied2(iem) = ivcur
            iedn(iem) = 0
         endif
c
         write(22,*) vx(ivcur), vy(ivcur), iem
         if(ivcur.ne.ivini) go to 1080
 1100 continue
      write(22,515) izn, izn, izn, lline
      close(unit=22)
c
c
c*********************************************************************
c
c     save coordinates of points that have moved in arrays x1 and y1
c
      do 1155 j = 1, nd
         i = idl(j)
         x1(i) = x2(i)
         y1(i) = y2(i)
 1155 continue
      go to 595
c
 1160 continue
c
c     create axis file
c
      open(unit=19,file=az)
      vxmin = vxmin-10.0
      vxmax = vxmax+10.0
      vymin = vymin-10.0
      vymax = vymax+10.0
      write(*,*)' vxmin = ',vxmin,' vxmax = ',vxmax
      write(*,*)' vymin = ',vymin,' vymax = ',vymax
      write(19,290) fas, nv, nv
      write(19,*) vxmin, vxmax
      write(19,*) vymin, vymax
      write(19,510) nv, nv, lline
      close(unit=19)
c
c     create script file for doing actual plots
c
      open(unit=25,file='vorplot.m')
      write(25,*)'clear all;'
      write(25,*)'hold on'
      write(25,*)'az;'
c
c     axis data
c
      write(25,*)'axis([as(2,1) as(2,2) as(3,1) as(3,2)]);'
      write(25,*)'pause;'
c
c     first file data
c
      idcr = 1
      write(25,*)'pa1;'
      write(25,*)'va1;'
      write(25,*)'n=pc1(1,1);'
      write(25,*)'for i=1:n'
      write(25,*)'l=pc1(i+1,3);'
      write(25,*)"hp(l)=plot(pc1(i+1,1),pc1(i+1,2),'.');"
      write(25,*)'end'
      write(25,*)'k=1;'
      write(25,*)'m=vc1(1,1);'
      write(25,*)'for i=1:m'
      write(25,*)'j=k+2;'
      write(25,*)'k=j-1+vc1(j-1,1);'
      write(25,*)'for n=j:k-1'
      write(25,*)'l2=vc1(n+1,3);'
      write(25,*)'hv(l2)=plot(vc1(n:n+1,1),vc1(n:n+1,2));'
      write(25,*)'end'
      write(25,*)'end'
      write(25,*)'iv=1'
      write(25,*)'mv(1)=getframe;'
c
c     next file data
c
 1210 format(a5,';')
 1220 format('n=',a5,'(1);')
 1230 format('l=',a5,'(i+1);')
 1240 format('m=',a5,'(1);')
 1250 format('k=j-1+',a5,'(j-1);')
 1270 format('l2=',a5,'(n+1);')
 1280 format('n=',a5,'(1,1);')
 1290 format('l=',a5,'(i+1,3);')
 1300 format("hp(l)=plot(",a5,"(i+1,1),",a5,"(i+1,2),'.');")
 1310 format('m=',a5,'(1,1);')
 1320 format('k=j-1+',a5,'(j-1,1);')
 1340 format('l2=',a5,'(n+1,3);')
 1350 format('hv(l2)=plot(',a5,'(n:n+1,1),',a5,'(n:n+1,2));')
 1400 continue
      if(idcr.ge.idyn) go to 1500
      idcr = idcr+1
c     write(25,*)'pause(1);'
      write(25,1210)spe(idcr)
      write(25,1210)sve(idcr)
      write(25,1220)spd(idcr)
      write(25,*)'for i=1:n'
      write(25,1230)spd(idcr)
      write(25,*)'delete(hp(l));'
      write(25,*)'end'
      write(25,*)'k=1;'
      write(25,1240)svd(idcr)
      write(25,*)'for i=1:m'
      write(25,*)'j=k+2;'
      write(25,1250)svd(idcr)
      write(25,*)'for n=j:k-1'
      write(25,1270)svd(idcr)
      write(25,*)'delete(hv(l2));'
      write(25,*)'end'
      write(25,*)'end'
      write(25,1210)spa(idcr)
      write(25,1210)sva(idcr)
      write(25,1280)spc(idcr)
      write(25,*)'for i=1:n'
      write(25,1290)spc(idcr)
      write(25,1300)spc(idcr),spc(idcr)
      write(25,*)'end'
      write(25,*)'k=1;'
      write(25,1310)svc(idcr)
      write(25,*)'for i=1:m'
      write(25,*)'j=k+2;'
      write(25,1320)svc(idcr)
      write(25,*)'for n=j:k-1'
      write(25,1340)svc(idcr)
      write(25,1350)svc(idcr),svc(idcr)
      write(25,*)'end'
      write(25,*)'end'
      write(25,*)'iv=iv+1'
      write(25,*)'mv(iv)=getframe;'
      go to 1400
c
 1500 continue
      write(25,*)'hold off'
      write(25,*)'pause;'
      write(25,*)'cla'
      write(25,*)'movie(mv)'
      close(unit=25)
c
      stop
      end