      program  polygn
c
*touched 1997/01/28 time=16:24:55 for a delaunay tessel.
c  %%%%%%% date 1988.10.20: time 10:54:11 %%%%%%%%%%%%%%%
c
c  **************************************************************
c  *                                                            *
c  *         pattern  analysis  of  two-dimensional             *
c  *         system of points by  voronoi  polygons             *
c  *                                                            *
c  *          <<  main program : version 4.3 >>                 *
c  **************************************************************
c  *     <copyright: masaharu tanemura, i.s.m.  1986.12.20>     *
c  **************************************************************
c
c
c     *** comments on parameters ***
c
c        dens  :   number density(= number of particles/area of box)
c                  of a pattern to be processed.
c        eps   :   critical value which gives an allowance range of
c                  error when a measure-zero phenomenon happens
c                  (usually eps = 1.0d-12 is satisfactory).
c        np    :   number of particles of the pattern to be processed.
c                      3 < np < 1000  is assumed.
c        nef   :   efficient number of particles in the construction of
c                  voronoi polygons for respective particles.
c                              3 < nef < np
c
c                  nef [ several tens is usually effective ].
c        nrg   :   average number of particles in a circle of efficient
c                  range(abbreviated by cer).   this value is used to
c                  determine the radius of cer.
c                               nef < nrg < np.
c
c        nclx1, ncly1   :
c                  numbers of cells of the rectangular box in  x & y
c                  directions, respectively.   the box is divided into
c                  (nclx1 x ncly1) identical cells.
c
c    << inputs from file 2 >>
c        *  single data set is composed of the following list.
c                       (number of record is 1)
c                 (np)
c                 (idpat)    pattern matching identifier
c                 (bxlx, bxly)
c                            side-lengths of the rectangular box
c                 (dens)     number density
c                 (xx(i),yy(i),i=1,np)
c                       coordinates of particles confined in the box
c
c    << outputs into file 1 >>
c        *  single data set is composed of np+1 records.  the first
c           record is a simple duplication of the record of tape 1, and
c           following np records contain informations about the voronoi
c           polygons of respective particles.
c        <1>  list of the first record
c                 (np)
c                 (idpat)
c                 (bxlx,bxly)
c                 (dens)
c                 (xx(i),yy(i),i=1,np)
c        <k>  list of the k-th record  (k=2,np+1)
c                 (ln) : particle number
c                 (nedg)   :  no. of edge
c                 (spol)   :  area of the polygon
c                 (ivw(l),l=1,nedg) : set of particles contiguous to ln.
c                 (dstnc(l),l=1,nedg) : distance to ivw(l)
c                 (spj(l),l=1,nedg) : length of edges
c                 (txj(l),tyj(l),l=1,nedg) : coordinates of vertices
c                 (ts(l),l=1,nedg) : area of delaunay triangle
c                 (radt(l),l=1,nedg) : radius of circumscribing circle
c                          of each d.t.
c                 (thet(l),l=1,nedg) : interior angle of the polygon
c
c       ---------------------------------------------------------------
c
c
c
      implicit real*8 (a-h,o-p,r-z)
      integer*4  mnbr
c
      parameter (nvmx = 45)
      parameter (npnt = 12345)
      parameter (ndln = 34567)
*
      dimension radt(nvmx),ts(nvmx),spj(nvmx),thet(nvmx),dstnc(nvmx)
      dimension knbr(nvmx)
      dimension icnt(2), jcnt(3,ndln), kcnt(ndln)
      common /vtxy/ vtx(nvmx,npnt),vty(nvmx,npnt)
      common /mnbr/ mnbr(nvmx,npnt),nbr(npnt)
      common /xy/ xx(npnt),yy(npnt) /boxl/bxlx,bxly,dens,np,idpat
      common /cnst/bx2,by2, pi,pi2,pi5,eps, bxlw,bxup,bylw,byup
      common /in/in1,in2,in3,in4
      common /iperio/ iperio
      common /ilist/ ilist
      data   kcnt /ndln*1/
c
***********************************************************************
c            comments to input parameters                             *
***********************************************************************
c
c     ilist  = 1, then output of detailed information
c            = 2, then output of detailed information for test
c            = 3, then output of vertex coordinates for plot
c            = 0, then output of only the final result
c     ifile  = 1, then inf. of voronoi polygons are stored in file (1).
c                 otherwise, data are not stored.
c     iperio = 1  periodic boundary conditions are assigned.
c            = 0 ,then division of box is done under the existence of
c                 boundaries.
***********************************************************************
*            indices for the case: iperio = 0                         *
*---------------------------------------------------------------------*
*        index                  boundary edge                         *
*---------------------------------------------------------------------*
*         -1     the edge parallel with y-axis  and smallest in x     *
*         -2     the edge parallel with y-axis  and largest  in x     *
*         -3     the edge parallel with x-axis  and smallest in y     *
*         -4     the edge parallel with x-axis  and largest  in y     *
*=====================================================================*
c     incard = 1  coordinate data are read as card image (80 bytes
c                 record sizes)
c            = 0  data are read from file (2).
c     (bxlw,bxup)    boundary values for x-direction
c     (bylw,byup)    boundary values for y-direction
c
c
c     *********** input of parameters **********
c
      inum = 0
      open(3,file='delone1.conf')
c     open(3,file='tane20p.aftkobe.data')
c     open(3,file='NChina.data')
c     open(3,file='cjapan.data')
c     read(3, *)  ilist
*     read(3, *)  ilist
      ilist=1
c        ilist = 2
c     *************  1988/10/20 *****
c     read(3, *)  ifile
*     read(3, *)  ifile
      ifile=0
*       ifile = 2
c     read(3, *)  iperio
*     read(3, *)  iperio
      iperio=0
c     read(3, *)  incard
*     read(3, *)  incard
      incard=1
c     read(3, *)  bxlw,bxup
*     read(3, *)  bxlw,bxup
c     read(3, *)  bylw,byup
*     read(3, *)  bylw,byup
      bxlw=0
      bylw=0
c     read(3, *)  dens
*     read(3, *)  dens
      dens=1.d0
      read(3, *)  eps
*     read(3, *)  eps
c     read(3, *)  np
*     read(3, *)  np
      read(3, *)  nef0, nrg0
*     read(3, *)  nef0, nrg0
c     read(3, *)  nclx1, ncly1
*     read(3, *)  nclx1, ncly1
      nclx1=3
      ncly1=3
c     read(3, *)  idpat
*     read(3, *)  idpat
      idpat=1
c
c     ******************************************
c
c
c     *******  output of parameters  ***********
c
c     open(9,file='tane20p.out')
      open(9,file='delone1.out')
      write(6,60) iperio, np, dens, eps, nef0, nrg0, nclx1, ncly1
     1           ,ilist, ifile, incard, idpat
c
c
   60 format(1h0,10x,'***** input parameters *****'/
     1   1h ,10x,'iperio=',i10,5x,'np    =',i10/
     2   1h ,10x,'dens  =',f10.5,5x,'eps=',e15.5/
     3   1h ,10x,'nef0  =',i10,5x,'nrg0  =',i10/
     4   1h ,10x,'nclx1 =',i10,5x,'ncly1 =',i10/
     5   1h ,10x,'ilist =',i10,5x,'ifile =',i10/
     6   1h ,10x,'incard=',i10,5x,'idpat =',i10/)
c
c     ***********************************************************
      nrg=nrg0
c
c
c     ******** input of configuration ********
c
c     -----------------------------------
      call config (idpat, idpat1, incard)
c     -----------------------------------
c
c     ****************************************
c
      if(ilist.eq.1) then
         lpage = 1
         write(6,61) np,idpat,bxlx,bxly,dens
c        write(6,62) (l,xx(l),yy(l),l=1,np)
         write(6,63) lpage
   61 format(1h0,10x,'*** input coordinates ***'/1h ,10x,'np=',i6
     1   ,3x,'idpat=',i5,5x,'bxlx,bxly=',2f15.5,3x,'dens=',f15.5)
*  62 format(1h0,4(3x,'no.',7x,'x',11x,'y',4x)/(1h ,4(i5,2f12.5,1x)))
   62 format(1h0,3x,'no.',7x,'x',11x,'y'/(1x ,i5,2f12.5,1x))
c  63 format(1h1,'*** detailed outputs ***',t120,'page',i6/1h0,4x,'n',
   63 format(1h0,'*** detailed outputs ***',t120,'page',i6/1h0,4x,'n',
     1 2x,'nedge',2x,'polygon area',6x,4('nbr',5x,'vertex coord.',4x))
c
         lprn = 0
      end if
c
      factr = 1.0d0
c**** transform of coordinates *******
c     factr = bxlx / bxly * 0.3
      do 123 kn=1,np
         xx(kn) = xx(kn) - bxlw
         yy(kn) = (yy(kn) - bylw) * factr
 123  continue
c******************************** 88.10.20 **
      bxup = bxup - bxlw
      byup = (byup - bylw)*factr
      bxlw = 0.0
      bylw = 0.0
      bxly = bxly * factr
      dens = dfloat(np)/bxlx/bxly
c
      densa=dens
      boxar=bxlx*bxly
      bx2=0.5d0*bxlx
      by2=0.5d0*bxly
      pi=3.14159265358979324d0
      pi2=pi+pi
      pi5=pi*0.5d0
      if(dens.le.0.0) densa=dfloat(np)/boxar
      nclx=nclx1
      ncly=ncly1
      bra=dfloat(nrg)/pi/densa
      brasq=dsqrt(bra)
c
c        write(6,1066) np,bxlx,bxly,dens,brasq
 1066 format(1h ,3x,'np=',i6,' bxly,bxly =',2f15.6,'  dens=',f10.5/
     \       1h ,3x,'brasq=',f15.6)
c        write(6,62) (l,xx(l),yy(l),l=1,np)
c
c
      if(ifile .ne. 0) then
c     **********************************************************
        if(ifile .eq. 1) then
          write(1)  np, idpat, bxlx,bxly, dens, (xx(l),yy(l),l=1,np)
        else
          write(1,6601)  np, idpat, bxlx,bxly, dens
          write(1,6602)  (xx(l),yy(l),l=1,np)
        endif
c     **********************************************************
 6601 format(2i10,3f20.10)
 6602 format((4f20.10))
      endif
c
c     ---------------------------------------
      call vorp (nef0, bra, brasq, nclx,ncly)
c     ---------------------------------------
c
      sumar=0.0d0
c
      do 900 k=1,np
      write(6,*) 'step =', k,'/',np
c
c     *****************************************************
c     *                                                   *
c     *      computation of geometrical quantities        *
c     *   concerning the voronoi polygon of atom  <k>.    *
c     *                                                   *
c     *****************************************************
c
        nedg=nbr(k)
        spol=0.0d0
        nedg1=nedg-1
        do 600 i=1,nedg
          knbr(i)=mnbr(i,k)
          ngbr1=knbr(i)
          ngbr2=mnbr(1,k)
          if(i.le.nedg1) ngbr2=mnbr(i+1,k)
          if(iperio.eq.1) go to  300
          nbd1=ngbr1-np
          nbd2=ngbr2-np
          if(nbd1.le.0) go to  300
c
c         ---------------------------------------
          call mirror (nbd1, xx(k),yy(k), xd1,yd1)
c         ---------------------------------------
c
  200     if(nbd2.le.0) go to 400
c
c         ---------------------------------------
          call mirror (nbd2, xx(k),yy(k), xd2,yd2)
c         ---------------------------------------
c
          go to 500
  300     xd1=xx(ngbr1)-xx(k)
          yd1=yy(ngbr1)-yy(k)
          if(iperio.eq.0) go to 200
  400     xd2=xx(ngbr2)-xx(k)
          yd2=yy(ngbr2)-yy(k)
          if(iperio.eq.0) go to 500
c
c         --------------------
          call period (xd1,yd1)
          call period (xd2,yd2)
c         --------------------
c
  500     continue
          rr1=xd1*xd1 + yd1*yd1
          rr2=xd2*xd2 + yd2*yd2
          dstnc(i)=dsqrt(rr1)
          rr3=(xd1-xd2)*(xd1-xd2) + (yd1-yd2)*(yd1-yd2)
          r12=dsqrt(rr1*rr2)
          cost=(rr1+rr2-rr3)*0.5d0/r12
          cos1=1.-cost*cost
          sint=dsqrt(cos1)
          thet(i)=datan2(sint,cost)
          ts(i)=r12*0.5d0*sint
          radt(i)=r12*dsqrt(rr3)*0.25d0/ts(i)
  600   continue
c
        do 700 j=1,nedg1
          t1=vtx(j,k)-vtx(j+1,k)
          t2=vty(j,k)-vty(j+1,k)
          rr=t1*t1+t2*t2
          spj(j)=dsqrt(rr)
          spol=spol+dstnc(j+1)*spj(j)*0.25d0
  700   continue
        t1=vtx(nedg,k)-vtx(1,k)
        t2=vty(nedg,k)-vty(1,k)
        rr=t1*t1+t2*t2
        spj(nedg)=dsqrt(rr)
        spol=spol+dstnc(1)*spj(nedg)*0.25
c
        sumar=sumar+spol
c
        do 750 l = 1, nbr(k)
          if(mnbr(l,k) .gt. np) then
            knbr(l) = np - mnbr(l,k)
          else
            knbr(l) = mnbr(l,k)
          end if
  750   continue
c
        if(ifile .ne. 0) then
c        ----------   output into file (1)  ---------
c     **************************************************************
c     = in order to have detailed outputs, uncomment the following =
c     =    two lines and put a comment to the next one.            =
c     write(1)  k,nedg,spol,(knbr(l),dstnc(l),spj(l),l=1,nedg)
c    1     ,(vtx(l,k),vty(l,k), ts(l), radt(l), thet(l), l=1,nedg)
        if(ifile .eq. 1) then
          write(1) k,nedg,spol,(knbr(l),vtx(l,k),vty(l,k),l=1,nedg)
        else
          write(1,6603) k,nedg,spol
          write(1,6604) (knbr(l),vtx(l,k),vty(l,k),l=1,nedg)
        endif
c
 6603 format(2i10,f20.10)
 6604 format((2(i5,2f15.7,5x)))
c     **************************************************************
c
      endif
c
      if(ilist.eq.1) then
        if(nbr(k).le.4) then
          write(6,64) k,nbr(k),spol,
     1               (knbr(l),vtx(l,k),vty(l,k),l=1,nedg)
   64     format(1h0,i5,i7,f13.5,5x,4(i5,2f10.3))
        else
          write(6,65) k,nbr(k),spol,
     1               (knbr(l),vtx(l,k),vty(l,k),l=1,nedg)
   65 format(1h0,i5,i7,f13.5,5x,4(i5,2f10.3)/(1h ,30x,4(i5,2f10.3)))
        end if
        if(spol.gt.10.) write(6,*) 'spol=',spol,k,nbr(k)
        lprn = lprn + 2 + (nbr(k)-1)/4
        if(lprn.gt.50) then
          lpage = lpage + 1
          write(6,63) lpage
          lprn = 0
        end if
      end if
c
  900 continue
c
      write(6,66)
   66 format(1h0,10x,'***** result of voronoi division *****'/
     1     1h0,10x,'idpat',4x,'np',13x,'bxlx',11x,'bxly',
     2     11x,'brasq',9x,'sum of pol.ar.',12x,'box area')
      write(6,67) idpat,np,bxlx,bxly, brasq, sumar, boxar
   67 format(1h ,10x,i5,1x,i5,2x,2f15.4,1x,f15.3,3x,1p2e20.8)
*     if(dabs(sumar/boxar-1.0d0).gt.eps) write(6,68)
*  68 format(1h ,20x,'+++++ warning +++++   division is not completed')
c
c
      if(ilist .ne. 0) then
***** output of delaunay information ************** (97/01/28)
        do 970 k = 1, np
          do 975 l = 1, nbr(k)
            icnt(1) = mnbr(l,k)
            if(l .eq. nbr(k)) then
              icnt(2) = mnbr(1,k)
            else
              icnt(2) = mnbr(l+1,k)
            end if
            call out9(np, k, icnt(1), inum, jcnt, kcnt)
  975     continue
  970   continue
*
        write(6,665) inum
  665   format(/1x,'number of delaunay triangle = ',i7)
        do 9100 l = 1, inum
          if(jcnt(2,l) .lt. 0) then
            if(kcnt(l) .ne. 1) then
              write(6, 961) l, kcnt(l)
            end if
          else if(jcnt(1,l) .lt. 0) then
            if(kcnt(l) .ne. 2) then
              write(6, 961) l, kcnt(l)
            end if
          else if(jcnt(1,l).gt.0 .and. kcnt(l).ne.3) then
            write(6, 961) l, kcnt(l)
  961       format(1x,'error in kcnt l = ',i7,2x,' kcnt(l) = ',i7)
          end if
c 99/05/24
          write(9, 960) l, (jcnt(j, l),j=1,3), kcnt(l)
  960     format(1x,i7,1x,3i7,1x,i7)
 9100   continue
***** output of delaunay information ************** (97/01/28)
      end if
*
      if(ilist .eq. 3) then
***** output for drawing voronoi tessellation ***** (97/01/28) *****
        write(6, 901)
  901   format(//)
        do 930 k = 1, np
          do 931 kv = 1, nbr(k)
            kv1 = kv + 1
            if(kv1 .gt. nbr(k)) kv1 = 1
            write(6, 902) vtx(kv,k),vty(kv,k), vtx(kv1,k),vty(kv1,k)
  902       format(2x,2f15.8,2x,2f15.8)
  931   continue
  930 continue
***** output for drawing voronoi tessellation ***** (97/01/28) *****
      end if
      stop
      end
c
      subroutine  vorp (nef0, bra,brasq, nclx,ncly)
c
c  **************************************************************
c  *                                                            *
c  *          tessellation   of  a  planar  region              *
c  *               by  voronoi  polygons(#)                     *
c  *                                                            *
c  **************************************************************
c
c      note(#)   definition of 'voronoi polygon' will be given below.
c
c     this program computes the area, the number of edges and other
c     quantities  of voronoi polygon for each particle in the plane.
c     particles are assumed to be placed on a rectangular box of
c     sidelengths <bxlx> and <bxly>.
c            the box is uniquely divided, i.e. the boundary effect is
c            considered.
c
c
c     --------------------------------------------------------------
c      following subroutines are directly called in this program :
c
c                            config
c                            rebox
c                            blking
c                            moda
c                            qsort
c                            smlteq
c                            innova
c     --------------------------------------------------------------
c
c       ---------------------------------------------------------------
c        definition of voronoi polygon of particle i  :  v(i)
c             v(i):   a set of all points  x  in the plane whose
c                  distance from x(i) is less than those from other
c                  x(j)'s (j=1,...,np; not including i).   then,
c                  v(i) is the set of all points  x  satisfying
c                         ! x - x(i) !  <  ! x - x(j) !
c                                  (j=1,...,i-1,i+1,...,np).
c
c       ---------------------------------------------------------------
c
c              *********************************************
c              this program was
c              written     by m.tanemura          1975. 6.15
c              revised version 2.0                1977.11.25
c              revised version 3.0                1982.12.09
c              revised version 3.1                1983.03.04
c              revised version 4.0                1984.04.16
c              revised version 4.1                1984.09.10
c              revised version 4.2                1986.12.20
c              *********************************************
c
c
c
c      inord(l) : inverse relation of k - iord(l)
c      next(k) : storage of sequential treatment of atoms
c      mclxy(1,k) : x-component of cell index for the atom <k>
c      mclxy(2,k) : y-component of cell index for the atom <k>
c      icomp(k) = 0, if the atom <k> is not completed by delaunay
c                 triangles (dt's).
c      icomp(k) = 1, if the atom <k> is completed by dt's; then, the
c                 atom does not appear for the after processing.
c      nbr(k) : no. of neighbor atoms at the instant
c      inbr(l,k) (l=1,...,nbr(k)) : index for the contiguity to the
c                          neighbor atoms
c              = 0, initial condition
c              = 1, singlet contiguity ( the edge (k - mnbr(l,k)) is
c                singly counted)
c              = 2, doublet contiguity ( the edge is doubly counted)
c        <note>  if inbr(k,l) = 2, then the edge is completed.
c      ndt(k) : no. of delaunay triangles (dt's) of <k> at the instant.
c      mpair(1,l,k),mpair(2,l,k) (l=1,...,ndt(k)) : a pair of atoms
c              which construct with <k> a  d.t.
c       <note>  (k, mapir(1,l,k), mpair(2,l,k) ) is a set of atoms which
c              are contiguous with one another.
c      mnbr(l,k) (l=1,...,nbr(k)) : atoms which are neighbors to <k>.
c
c
c
c
      implicit real*8 (a-h,o-p,r-z)
      integer*4  nlb,lb, mclxy,inbr,mnbr,mpair
c
      parameter (nwrk = 22345)
      parameter (nvmx = 45)
      parameter (npnt = 12345)
      dimension aa(4),bb(2)
      dimension rd(nwrk),rd1(nwrk), tx(nwrk),ty(nwrk),
     1          x0(nwrk),y0(nwrk), fi(nwrk), dis(nwrk)
      dimension iord(nwrk),jjj(nwrk),jjj1(nwrk),mord(nwrk),
     2          ivw(nwrk),iv(nwrk),itobu(nwrk),kord(nwrk)
      dimension txj(nwrk),tyj(nwrk)
      dimension next(npnt), inord(npnt)
      dimension kb(3),jb(3),  knd(3),jnd(3)
      dimension mclxy(2,npnt),inbr(nvmx,npnt), mpair(2,nvmx,npnt),
     1          mdt(nwrk), icomp(npnt), ndt(npnt)
c
c
      common /vtxy/ vtx(nvmx,npnt),vty(nvmx,npnt)
      common /mnbr/ mnbr(nvmx,npnt),nbr(npnt)
c     common /lb/ nlb(5,5), lb(18296,5,5)
      common /lb/ nlb(9,9), lb(5555,9,9)
      common /xy/ xx(npnt),yy(npnt) /boxl/bxlx,bxly,dens,np,idpat
      common /cnst/bx2,by2, pi,pi2,pi5,eps, bxlw,bxup,bylw,byup
      common /iperio/ iperio
      common /in/in1,in2,in3,in4
      common /ilist/ ilist
c
c
      nef=nef0
c
c     ----------
      call rebox
c     ----------
c
  100 continue
c
c     ------------------------------
      call blking (nclx,ncly, mclxy)
c     ------------------------------
c
      if(iperio.eq.0) go to 300
      nlmin=nlb(1,1)
      do 200 j=1,ncly
        do 200 i=1,nclx
          if(nlb(i,j).lt.nlmin)  nlmin=nlb(i,j)
  200 continue
c     if(nlmin*9.gt.nrg)  go to 300
c     if(nclx.eq.3.or.ncly.eq.3) go to 300
c     nclx=nclx-1
c     ncly=ncly-1
c     go to 100
  300 continue
c
      if(mod(nclx,2).eq.0) nxg=nclx/2
      if(mod(nclx,2).eq.1) nxg=(nclx+1)/2
      if(mod(ncly,2).eq.0) nyg=ncly/2
      if(mod(ncly,2).ne.0) nyg=(ncly+1)/2
c**************************** 1988/10/20 ********
c     nxg = 1
c     nyg = 2
c***  for clustered data (earthquake data) ******
      nlxyg=nlb(nxg,nyg)
      if(nlxyg.le.0) stop 5454
      nlxyg2 = nlxyg / 2
      if(nlxyg2 .le. 0) nlxyg2 = 1
c     lnc=lb(1,nxg,nyg)
      lnc = lb(nlxyg2,nxg,nyg)
c     write(6,1065) nxg,nyg,nlxyg, lnc
 1065 format(1h ,5x,'nxg,nyg,nlxyg, lnc=',4i6)
c**************************** 1988/10/20 ********
      if(nlxyg.eq.1) go to 410
      xg=bxlx/dfloat(nclx)*(dfloat(nxg)-0.5d0)
      yg=bxly/dfloat(ncly)*(dfloat(nyg)-0.5d0)
      dmin=1.0d60
      do 400 l=1,nlxyg
        ln=lb(l,nxg,nyg)
        dx=xx(ln)-xg
        dy=yy(ln)-yg
        dd=dx*dx+dy*dy
        if(dd .lt. dmin) then
          lnc = ln
          dmin = dd
        endif
  400 continue
  410 continue
c
      lnex=1
      if(ilist .ge. 2) then
c         write(6,1065) nxg,nyg,nlxyg, lnc
c         write(6,1064) xg,yg,dmin
 1064     format(1h ,5x,'xg,yg,dmin =',3f15.6)
      endif
      next(1)=lnc
c    ***** <lnc> is the atom which is processed at the first *****
c
      lcount=1
c
c     ***** preparation before entering to the outer loop *****
c
      do 510 l=1,np
        icomp(l)=0
        nbr(l)=0
        ndt(l)=0
        do 500 k=1,nvmx
  500   inbr(k,l)=0
  510 continue
c
      ln=lnc
      if(ilist.eq.2) write(6,66)
   66 format(1h0/)
      if(ilist.eq.2) write(6,60) icomp(ln), ln
c
c     ****** check whether the voronoi polygon for the atom <ln> ****
c     ****** is already obtained                                 ****
c
      idt=0
      index1=0
      nbr(ln)=0
      ndt(ln)=0
 1000 continue
c     ********** outer loop **********
c
c
c     ******  centre atom ******
      xc=xx(ln)
      yc=yy(ln)
      nxg=mclxy(1,ln)
      nyg=mclxy(2,ln)
      nlxyg=nlb(nxg,nyg)
c
      if(icomp(ln).eq.1) write(6,60) icomp(ln), ln
   60 format(1h0,50x,'*** icomp=', i3,2x,' ln=',i5,' ***')
      if(icomp(ln).eq.1) stop 8000
c
      jn=0
      ieff=0
c         ieff : effective no. of neighbor set
c            (set of atoms for which icomp(k) = 0 )
c         jn   : no. of neighbor set which includes those icomp(k)=1.
      jnef=0
      bra1=bra
      brasq1=brasq
      fi1=0.0d0
c ********************* 1988/10/21 *******
      if(ilist .ge. 2) then
c       write(6,1062) xc,yc,nxg,nyg,ndt(ln)
 1062   format(1h ,2x,'xc,yc,nxg,nyg,ndt(ln) =',2f12.5,3i6)
      endif
      if(ndt(ln).le.0) go to 2000
c
c
c ******************* case where <ln> has some d.t.'s *****************
c     ***** case (2) nbr(ln).gt.0 *****
      nbr1=nbr(ln)
      kdt=ndt(ln)
c
      l1=1
      if(nbr1.eq.kdt) go to 1110
      do 1100 l1=1,nbr1
        if(inbr(l1,ln).eq.1) go to 1110
 1100 continue
      l1=1
 1110 mnbr1=mnbr(l1,ln)
c
c
      ivw(1)=mnbr1
      ie = 1
      idt=0
      index1=0
 1200 do 1230 l2=1,kdt
        ivw1=0
        if(index1.eq.0) go to 1220
        do 1210 ll2=1,idt1
          if(l2.eq.mdt(ll2)) go to 1230
 1210   continue
 1220   if(ivw(ie).eq.mpair(1,l2,ln)) ivw1=mpair(2,l2,ln)
        if(ivw(ie).eq.mpair(2,l2,ln)) ivw1=mpair(1,l2,ln)
        if(ie.eq.1.and.ivw1.gt.0) go to 1240
        if(ivw1.eq.0) go to 1230
        if(ivw1.ne.ivw(ie-1)) go to 1240
 1230 continue
      if(index1.eq.1) go to 3100
      if(ilist .eq. 2) then
        write(6,6993) ln,kdt,ivw1,ie,mnbr1,idt,nbr1
 6993   format(1h ,'  ln,kdt,ivw1,ie,mnbr1,idt,nbr1=',7i8)
        write(6,6994) (inbr(ll,ln),ll=1,nbr1)
        write(6,6994) (mnbr(ll,ln),ll=1,nbr1)
 6994   format((1h ,5x,10i5))
      endif
      stop 1230
 1240 idt=idt+1
      mdt(idt)=l2
      if(ivw1.eq.ivw(1)) go to 1270
      ie=ie+1
      ivw(ie)=ivw1
      do 1250 l3=1,nbr1
        if(mnbr(l3,ln).eq.ivw1) go to 1260
 1250 continue
      stop 1250
 1260 if(inbr(l3,ln).ne.2) go to 1270
      if(idt.eq.kdt) go to 1270
      go to 1200
 1270 continue
      if(index1.eq.1) go to 3100
c
c     *****
      fi1=0.0d0
      do 1370 ii=1,ie
        jn=jn+1
        if(ii.gt.2) go to 1310
        if(ii.eq.2) fi1=fi(1)
 1310   continue
        ieff=ieff+1
        nbd=ivw(ii)-np
        if(nbd.le.0) go to 1320
c
c       ------------------------------
        call mirror (nbd, xc,yc, xd,yd)
c       ------------------------------
c
        dis(ieff)=xd*xd+yd*yd
        x0(ieff)=xd
        y0(ieff)=yd
        iord(ieff)=ivw(ii)
        inord(ivw(ii))=ieff
        if(nbd.eq.1) f0=pi-fi1
        if(nbd.eq.2) f0=-fi1
        if(nbd.eq.3) f0=pi+pi5-fi1
        if(nbd.eq.4) f0=pi5-fi1
        go to 1350
 1320   xd=xx(ivw(ii))-xc
        yd=yy(ivw(ii))-yc
        if(iperio.eq.0) go to 1330
c
        call period (xd,yd)
c
 1330   rr2=xd*xd+yd*yd
        rr=dsqrt(rr2)
        dis(ieff)=rr2
        x0(ieff)=xd
        y0(ieff)=yd
        iord(ieff)=ivw(ii)
        inord(ivw(ii))=ieff
        if(xd.lt.0.0d0) go to 1340
        if(yd.ge.0.0d0) f0=dasin(yd/rr)-fi1
        if(yd.lt.0.0d0) f0=pi2-dasin(-yd/rr)-fi1
        go to 1350
 1340   f0=pi-dasin(yd/rr)-fi1
 1350   if(ii.eq.1) go to 1360
        if(f0.lt.0.0d0) f0=f0+pi2
 1360   fi(ieff)=f0
 1370 continue
c
c
c
      ie2=ie+1
      fi(1)=0.0d0
      indx0=0
      if(fi(2).gt.fi(ie)) indx0=1
c
 2000 continue
c
c
      call moda (nxg-1, jb(1), nclx, jnd(1), iperio)
      jb(2)=nxg
      jnd(2)=0
      call moda (nxg+1, jb(3), nclx, jnd(3), iperio)
      call moda (nyg-1, kb(1), ncly, knd(1), iperio)
      kb(2)=nyg
      knd(2)=0
      call moda (nyg+1, kb(3), ncly, knd(3), iperio)
      do 2160 kko=1,3
        kbb=kb(kko)
        if(kbb.lt.0) go to 2160
        dnk=knd(kko)*bxly
        do 2150  jjo=1,3
          jbb=jb(jjo)
          if(jbb.lt.0) go to 2150
          dnj=jnd(jjo)*bxlx
          ijka=nlb(jbb,kbb)
          if(ijka.le.0)  go to 2150
          do 2140 lna=1,ijka
            lnc=lb(lna,jbb,kbb)
            if(lnc.eq.ln)  go to 2140
            if(nbr1.le.0) go to 2060
            do 2050 lk=1,nbr1
 2050       if(lnc.eq.mnbr(lk,ln)) go to 2070
 2060       if(icomp(lnc).eq.1) go to 2140
 2070       continue
            xd=xx(lnc)+dnj
            yd=yy(lnc)+dnk
            dx=xd-xc
            dy=yd-yc
            if(dabs(dx).gt.brasq1)  go to 2140
            if(dabs(dy).gt.brasq1)  go to 2140
            dist=dx*dx+dy*dy
            if(dist.gt.bra1)  go to 2140
            jn=jn+1
            if(dx.lt.0.0d0) go to 2100
            if(dy.ge.0.0d0) f0=dasin(dy/dsqrt(dist))-fi1
            if(dy.lt.0.0d0) f0=pi2-dasin(-dy/dsqrt(dist))-fi1
            go to 2110
 2100       f0=pi-dasin(dy/dsqrt(dist))-fi1
 2110       continue
            if(f0.lt.0.0d0) f0=f0+pi2
            if(lcount.le.1) go to 2130
            if(indx0.eq.1) go to 2120
            if(f0.le.fi(ie)) go to 2140
            go to 2130
 2120       if(f0.ge.fi(ie)) go to 2140
            f0=-f0
 2130       ieff=ieff+1
            iord(ieff)=lnc
            inord(lnc)=ieff
            dis(ieff)=dist
            x0(ieff)=dx
            y0(ieff)=dy
            fi(ieff)=f0
 2140     continue
 2150   continue
 2160 continue
      if(jn.ge.jnef)  go to 2200
      jnef=jnef+1
      if(jnef.gt.10)  stop 99099
      bra1=bra1*1.10d0
      brasq1=brasq1*1.10d0
      jn=0
      go to 2000
 2200 continue
c
      iepd=ieff
      if(iperio.eq.1) go to 3000
c
c     ********* inclusion of mirror images to the set <cer> ********
c
      do 2600 m=1,4
        if(ndt(ln).le.0) go to 2400
        do 2300 ii=1,ie
          if(ivw(ii).eq.np+m) go to 2600
 2300   continue
 2400   iepd=iepd+1
        iord(iepd)=np+m
        inord(np+m)=iepd
c
c       ----------------------------
        call mirror (m, xc,yc, xd,yd)
c       ----------------------------
c
        dis(iepd)=xd*xd+yd*yd
        x0(iepd)=xd
        y0(iepd)=yd
        if(m.eq.1) f0=pi
        if(m.eq.2) f0=0.0d0
        if(m.eq.3) f0=pi+pi5
        if(m.eq.4) f0=pi5
        fi(iepd)=f0
        if(lcount.le.1) go to 2600
        f0=f0-fi1
        if(f0.lt.0.0d0) f0=f0+pi2
        if(indx0.eq.1) f0=-f0
        fi(iepd)=f0
        if(ndt(ln).le.0) go to 2600
        if(indx0.eq.1) go to 2500
        if(f0.le.fi(ie)) iepd=iepd-1
        go to 2600
 2500   if(f0.ge.fi(ie)) iepd=iepd-1
 2600 continue
c
c
c      **************************************************************
c
 3000 continue
c
      if(ndt(ln).le.0) go to 4000
c
c
 3100 continue
      jtri=idt
      do 3200 l1=1,ie
        iv(l1)=inord(ivw(l1))
 3200 continue
c
      if(idt.gt.30) stop 1666
      if(index1.eq.0) idt2=1
      if(index1.eq.1) idt2=idt1+1
      do 3300 l2=idt2,idt
        mdt1=mdt(l2)
        txj(l2)=vtx(mdt1,ln)
        tyj(l2)=vty(mdt1,ln)
 3300 continue
      if(index1.eq.1) go to 6100
      index1=0
      if(idt.ne.kdt) index1=1
      idt1=idt
      k=idt
      if(k.eq.ie) go to 7100
c
c
 4000 continue
c
      do 4100 ll=1,iepd
 4100 kord(ll)=ll
      ie1=1
      if(ndt(ln).gt.0) ie1=ie2
c
c     -------------------------------
      call qsort(ie1,iepd, dis, kord)
c     -------------------------------
c
c   ***** check whether the atom <ln> has any neighbors *****
c
      if(nef.gt.iepd) nef=iepd
      if(iperio.ne.1.and.nef.lt.iepd-3) nef=iepd
c
c     ***************************************
c
      if(lcount.le.1) go to 5000
      do 4200 ll=1,nef
 4200 mord(ll)=kord(ll)
c
c     ------------------------------
      call qsort (ie1,nef, fi, mord)
c     ------------------------------
c
      mord(nef+1)=iv(1)
      nef=nef+1
      istrt=ie1
      go to 6000
c
c     ********* case (i)  first point processing **********
c
c
 5000 continue
      nr1=kord(1)
      iv(1)=nr1
      ivw(1)=iord(iv(1))
c
c
      t1=x0(nr1)
      t2=y0(nr1)
      aa(1)=-t1
      aa(2)=-t2
      bb(1)=-(t1*t1+t2*t2)*0.5d0
      in1=1
      in2=2
      in4=666
c
      fi1=fi(nr1)
      fi(nr1)=0.0d0
      rmin=1.0d50
c
c
      do 5200 j=2,nef
        nrj=kord(j)
        fi(nrj)=fi(nrj)-fi1
        if(fi(nrj).lt.0.0d0) fi(nrj)=fi(nrj)+pi2
        in3=j
        aa(3)=-x0(nrj)
        aa(4)=-y0(nrj)
        bb(2)=(x0(nrj)*aa(3)+y0(nrj)*aa(4))/2.
c
c       ---------------------------
        call smlteq(aa,bb,t1,t2,ind)
c       ---------------------------
c
        if(ind.eq.1)  go to 5100
        itobu(j)=0
        tx(j)=t1
        ty(j)=t2
        rd(j)=t1*t1+t2*t2
        if(rmin.gt.rd(j)) kr2=j
        if(rmin.gt.rd(j)) rmin=rd(j)
        go to 5200
 5100   itobu(j)=1
 5200 continue
c
c
c     ***********************************************
c     *                                             *
c     *    registration of the 1-st d.t.            *
c     *                                             *
c     ***********************************************
c
      jtri=1
c     radt(1)=dsqrt(rd(kr2))
      txj(1)=tx(kr2)+xc
      tyj(1)=ty(kr2)+yc
      nr2=kord(kr2)
      ie=2
      iv(2)=nr2
      ivw(2)=iord(iv(2))
      indx0=0
      if(fi(nr2).gt.pi) indx0=1
      mord(1)=nr1
      do 5300 ll=2,nef
        nrk=kord(ll)
        if(indx0.eq.1) fi(nrk)=-fi(nrk)
 5300 mord(ll)=nrk
c
c
c
c     --------------------------
      call qsort (2,nef,fi,mord)
c     --------------------------
c
      mord(nef+1)=nr1
      nef=nef+1
c
      if(ilist .eq. 2) then
        write(6,6006) ln,ivw(1),ivw(2),
     1    ndt(ln), ndt(ivw(1)), ndt(ivw(2)),
     2    nef, txj(1),tyj(1)
 6006   format(1h ,' *** test(19970128a) l1,l2,l3',3i10/
     1    1x,'ndt(ln), ndt(ivw1), ndt(ivw2) ',3i5/
     2    1x,'nef =',i5,' txj tyj = ',2f15.7)
      end if
c     ==============================================
      call innova (ln,ivw(1),ivw(2), inbr,mpair, ndt
     1             ,txj(1),tyj(1))
c     ==============================================
c
      k=1
      istrt=3
c
 6000 continue
c
c
c     **************** internal loop for dt's ************
c
c
 6100 k=k+1
      if(k-ie) 6200,6200,7000
 6200 nrk=iv(ie)
 6300 aa(1)=-x0(nrk)
      aa(2)=-y0(nrk)
      bb(1)=-(aa(1)*aa(1)+aa(2)*aa(2))*0.5d0
      nrj=iv(ie-1)
      tt=-x0(nrj)*aa(2)+y0(nrj)*aa(1)
      isn=0
      if(tt.ge.0.) isn=1
      jun=0
      jun1=0
      rmin=1.0d50
      rmax=-1.0d50
      do 6600 i=istrt,nef
        if(lcount.gt.1) go to 6400
        if(i.ge.2.and.itobu(i).eq.1)  go to 6600
 6400   nri=mord(i)
        if(nri.eq.nrk) go to 6600
        aa(3)=-x0(nri)
        aa(4)=-y0(nri)
        tt=aa(3)*aa(2)-aa(4)*aa(1)
        tab=dabs(tt)
        if(tab.lt.eps) go to 6600
        isn1=0
        if(tt.ge.0.) isn1=1
        if(isn.eq.isn1) go to 6600
        bb(2)=(x0(nri)*aa(3)+y0(nri)*aa(4))*0.5d0
        in3=i
        in4=7
c
c       ----------------------------
        call smlteq(aa,bb,t1,t2,ind)
c       ----------------------------
c
        if(ind.eq.1)  go to 6600
        tx(nri)=t1
        ty(nri)=t2
        tt=-t1*aa(2)+t2*aa(1)
        isnt=0
        if(tt.ge.0.) isnt=1
        if(isn.eq.isnt) go to 6500
        jun=jun+1
        rd(i)=t1*t1+t2*t2
        jjj(i)=nri
        if(rmin.gt.rd(i)) jnn=i
        if(rmin.gt.rd(i)) rmin=rd(i)
        go to 6600
 6500   jun1=jun1+1
        rd1(i)=t1*t1+t2*t2
        jjj1(i)=nri
        if(rmax.lt.rd1(i)) jnn1=i
        if(rmax.lt.rd1(i)) rmax=rd1(i)
 6600 continue
         if(jun1.eq.0) go to 6700
          if(jun.eq.0.and.jun1.eq.0)  stop 3717
c
           nrl=jjj1(jnn1)
             istrt=jnn1+1
              go to 6800
 6700         continue
c
             nrl=jjj(jnn)
         istrt=jnn+1
 6800  jtri=jtri+1
       txj(jtri)=tx(nrl)+xc
       tyj(jtri)=ty(nrl)+yc
c
c     ******** branch if new vertex is equal to iv(1) ********
c
       if(nrl.eq.iv(1))  go to 7000
c
c     ********************************************************
c
c
c
      ie=ie+1
      iv(ie)=nrl
      ivw(ie)=iord(iv(ie))
      if(ilist .eq. 2) then
c **************   1988/10/20 *****************************
        write(6,6992) ln,xx(ln),yy(ln),ie,iv(ie),ivw(ie)
 6992   format(1h ,i5,2f15.6,3i10)
c ************** earthquake data **************************
      endif
c
      if(ilist .eq. 2) write(6,6007) ln,iord(nrk),ivw(ie)
 6007 format(1h ,' *** test(19890405b) l1,l2,l3',3i10)
c     ==================================================
      call innova (ln,iord(nrk),ivw(ie), inbr,mpair, ndt
     1           ,txj(jtri),tyj(jtri))
c     ==================================================
c
      idt=idt+1
      mdt(idt)=ndt(ln)
      nrj=iv(ie-1)
      nrk=iv(ie)
      ie1=ie+1
      kdt=ndt(ln)
      if(index1.eq.1) go to 1200
      go to 6100
c
c************************************ dt's  have closed around <ln> ****
c
c
 7000 continue
c
      if(ilist .eq. 2) write(6,6008) ln, iord(iv(ie)),iord(iv(1))
 6008 format(1h ,' *** test(19890405c) l1,l2,l3',3i10)
c
c     =====================================================
      call innova (ln,iord(iv(ie)),iord(iv(1)), inbr,mpair,
     1          ndt,txj(jtri),tyj(jtri))
c     =====================================================
c
c
 7100 continue
c
c     *****************************************************
      jnef=0
      brasq1=brasq
      bra1=bra
      nedg=jtri
c
      icomp(ln)=1
c
c
      nedg=nbr(ln)
      do 7200 l=1,nedg
        mnbr(l,ln)=ivw(l)
        vtx(l,ln)=txj(l)
        vty(l,ln)=tyj(l)
 7200 continue
c
c
c
      if(ilist.eq.2) write(6,7206) ln,nedg,(ivw(l),txj(l),tyj(l),
     1        l=1,nedg)
 7206 format(1h ,5x,2i8,(1h ,4(i4,2f10.3)))
      lcount=lcount+1
      if(lcount.gt.np) go to 8000
c
c    ***** determine a new atom to be processed at the next step *****
c
      do 7500 kk=1,nedg
        ivw1=ivw(kk)
        if(icomp(ivw1).eq.1.or.ivw1.gt.np) go to 7500
        if(lnex.le.1) go to 7400
        do 7300 ll=1,lnex
          if(ivw1.eq.next(ll)) go to 7500
 7300   continue
 7400   continue
        lnex=lnex+1
        next(lnex)=ivw1
 7500 continue
c
c
      ln=next(lcount)
      nef=nef0
c ***** test ********** 1988/10/21 *******
      if(mod(lcount,1000).eq.0) write(6,1063) lcount,ln,nef
      if(ilist .eq. 2) then
        write(6,1063) lcount,ln,nef
 1063   format(1h ,3x,'lcount,ln,nef =',3i6)
      endif
c ********************* 1988/10/21 *******
c
c     ***** return to the top *****
      go to 1000
c
 8000 continue
c
c
c     ***** check whether the division is completed *****
c
      do 8100 kl=1,np
        if(icomp(kl).eq.1) go to 8100
        write(6,61) kl, nbr(kl), icomp(kl)
   61   format(1h ,10x,'ln=',i5,2x,'nbr(ln)=',i3,2x,'icomp=',i2)
 8100 continue
c
      return
      end
      subroutine  smlteq (a,b,t1,t2,index)
c
c
c     ******************************************
c     *                                        *
c     *  simultaneous equation of order one    *
c     *                                        *
c     ******************************************
c
      implicit real*8 (a-h,o-z)
      dimension a(4),b(2)
      common /in/in1,in2,in3,in4
c
c
      ad=a(1)*a(4)-a(2)*a(3)
      if(abs(ad).lt.1.0d-50)  go to 100
      index=0
      t1=b(1)*a(4)-b(2)*a(2)
      t2=b(2)*a(1)-b(1)*a(3)
      t1=t1/ad
      t2=t2/ad
      return
  100 if(in3.ne.1)  go to 200
      index=1
      write(6,10) ad,(a(l),l=1,4),b(1),b(2),in1,in2,in3,in4
   10 format(1h0,5x,'ad=',f10.3,10x'a=',4f10.3,5x'b=',2f10.3/
     1        10x,'in1='i5,5x,'in2='i5,5x'in3='i5,5x'in4='i5)
      return
  200 index=2
      t1=1.0d30
      t2=1.0d30
      return
      end
      subroutine rebox
c
c
c
c     ***********************************************
c     *                                             *
c     *     putting particles into original box     *
c     *       using periodic boundary condition     *
c     *                                             *
c     ***********************************************
c
c
      implicit real*8 (a-h,o-z)
c
      parameter (npnt = 12345)
      common /xy/ xx(npnt),yy(npnt) /boxl/bxlx,bxly,dens,np,idpat
      common /cnst/bx2,by2, pi,pi2,pi5,eps, bxlw,bxup,bylw,byup
c
      do 500 i=1,np
        if( xx(i) - bxlw .lt. 0.) go to 70
        xb = (xx(i) - bxlw)/bxlx
        ixb = xb
        xx(i) = xx(i) - bxlx*ixb
        go to 71
   70   xb = (xx(i) - bxlw)/bxlx
        ixb = -xb + 1
        xx(i) = xx(i) + bxlx*ixb
   71   if(yy(i) - bylw .lt. 0.) go to 72
        yb = (yy(i) - bylw) / bxly
        iyb = yb
        yy(i) = yy(i) - bxly*iyb
        go to 500
   72   yb = (yy(i) - bylw)/bxly
        iyb = -yb + 1
        yy(i) = yy(i) + bxly*iyb
  500 continue
      return
      end
*
      subroutine  config (idp, idp1, incard)
c
c
c     *************************************
c     *                                   *
c     *     input of coordinate data      *
c     *                                   *
c     *************************************
c
c     if(incard.eq.1) then read data as card image from file (2)
c     if(incard.eq.0) then read data as binary image from file (2).
c
      implicit real*8 (a-h,o,p,r-z), real*4 (q)
c
      parameter (npnt = 12345)
      common /xy/ xx(npnt),yy(npnt) /boxl/bxlx,bxly,dens,np,idpat
      common /cnst/bx2,by2, pi,pi2,pi5,eps, bxlw,bxup,bylw,byup
      dimension d(7,npnt)
      real*4 ux,uy
c
      if(incard.eq.1) go to 400
c
c   ****** input binary data from file(2) ********
      read(2,end=300)  np,idp1, bxlx,bxly,dens
     1           ,(xx(l),yy(l),l=1,np)
c
c   **********************************************
      if(idp.ne.idp1)  go to 300
c
c
      idpat = idp1
      return
  300 write(6,60)  np, idp1
   60 format(1h0,20x,'\\\\\  data  end  \\\\\, np=',i6,2x,'idpat1=',i5)
      stop 2222
c
  400 continue
c
c    ********* input card image data from file (3) ********
c
*     read(3, *)  bxlx,bxly
c     read(3, *)  xmin,ymin,amag0
      read(3, *)  xmin,ymin
      read(3, *)  bxlx,bxly
      bxup=bxlx
      byup=bxly
*     read(3, *)  np,idp1
c     read(3, *)  np,idp1
      read(3, *)  np
      idp1=1
      close(3)
c
      open(3,file='etas2aniso.out3')
      read(3,*)
c
      dens=np/bxlx/bxly
      if(idp.ne.idp1) go to 300
c     read(3, *)  (xx(l),yy(l),l=1,np); changed format (1990/05/08)
c     read(3, *)  (xx(l),yy(l),ndmmy,l=1,np)
c     read(5, *)  (xx(l),yy(l),dmmy,dmmy,l=1,np)
c     read(3, *)  (xx(l),yy(l),l=1,np)
c     read(3, 1)  (idum,xx(l),yy(l),(d(i,l),i=1,5), l=1,np)
      read(3, *)  (idum,xx(l),yy(l),(d(i,l),i=1,5), l=1,np)
c     read(3, *)  (idum,xx(l),yy(l),(d(i,l),i=1,6), l=1,np)
c
    1 format(i5,4f12.5,3f9.4)
c     call xydchk(xx,yy,np)
c     epp =5.0d-1
      epp=0.00166666
c     epp=0.1d-5
      epp=0.1d-2
*     write(9,1) np,bxlx,bxly,xmin,ymin
      ix=1874
      iy=265
      iz=367
      kk=0
c
      xeps=0.1d-10*bxlx
      yeps=0.1d-10*bxly
c
      do 710 k=1,np
c       if(d(1,k).lt.amag0) go to 710
c       kk=kk+1
c       xx(kk) = xx(k) - xmin
c       yy(kk) = yy(k) - ymin
c
        xxk=xx(k)-xmin
        yyk=yy(k)-ymin
        if(xxk.gt.xeps .and. xxk.lt.bxlx-xeps .and.
     &     yyk.gt.yeps .and. yyk.lt.bxly-yeps) then
          kk=kk+1
          xx(kk)=xx(k)-xmin
          yy(kk)=yy(k)-ymin
c
          do 730 i=1,5
            d(i,kk)=d(i,k)
  730       continue  
c           d(3,kk)=1.0
c           d(4,kk)=1.0
c           d(5,kk)=0.0
  777     continue
c
c         call unifor(ux)
c         call unifor(uy)
          ux=random(ix,iy,iz)
          uy=random(ix,iy,iz)
c
          xxkk = xx(kk) - epp*(ux - 0.5d0)
          yykk = yy(kk) - epp*(uy - 0.5d0)
          if(xxkk.lt.0.0 .or. xxkk.gt.bxlx) go to 777
          if(yykk.lt.0.0 .or. yykk.gt.bxly) go to 777
          xx(kk) = xxkk
          yy(kk) = yykk
c         write(9,1) kk,xx(kk),yy(kk),(d(i,k),i=1,5)
  776     format(i5,2f20.10)
        else
        write(6,*) k,xx(k),yy(k)
        endif
  710 continue
      np=kk
      write(6,*) 'np',np
c     stop
      write(9,1) np,bxlx,bxly,xmin,ymin
      do 720 k=1,np
        write(9,111) k,xx(k),yy(k),(d(i,k),i=1,5)
  720 continue
  111 format(i5,2f12.8,2f12.5,3f9.4)
c
c    ***********************************
c
      return
      end
c++++++++++++++++++++++++++++++++++++++++++++++++++++++
      function random(ix,iy,iz)
c ix,iy,iz are seeds for generating uniform random numbers
c
c::::     Wichmann+Hill Appl. Statist (31) 188-? (1982)
c
      implicit real*8 (a-h,o-z)
c     call unifor(random)
c     return
      ix=171*mod(ix,177)-2*(ix/177)
      iy=172*mod(iy,176)-35*(iy/176)
      iz=170*mod(iz,178)-63*(iz/178)
      if (ix.lt.0) ix=ix+30269
      if (iy.lt.0) iy=iy+30307
      if (iz.lt.0) iz=iz+30323
      random=mod(float(ix)/30269.0+float(iy)/30307.0+
     1            float(iz)/30323.0,1.0)
      return
      end
c++++++++++++++++++++++++++++++++++++++++++++++++++++++
      subroutine  moda (k,j,ncl,indx, ip)
c
c     **********************************
c     *                                *
c     *      find neighbor cells       *
c     *                                *
c     **********************************
c
c
      if(k.lt.1)  go to 10
      if(k.gt.ncl)  go to 20
      j=k
      indx=0
      return
   10 j=ncl+k
c     if(ip.ne.1) j=-j
      indx=-1
      if(ip.ne.1) indx=0
      return
   20 j=k-ncl
c     if(ip.ne.1) j=-j
      indx=1
      if(ip.ne.1) indx=0
      return
      end
      subroutine  blking (nclx,ncly, mclxy)
c
c
c     ************************************************
c     *                                              *
c     *   division of the original box into small    *
c     *        sub-cells.                            *
c     *                                              *
c     *     no. of cells =  nclx  *  ncly            *
c     *     no. of particles in the cell (i,j) :     *
c     *               nlb(i,j)                       *
c     *                ( nlb  <= 1000 is assumed )  *
c     *     member of atoms in the cell (i,j)  :     *
c     *              (lb(k, i,j),k=1,nlb(i,j) )      *
c     *                                              *
c     ************************************************
c
c
      implicit real*8 (a-h,o-z)
      integer*4  nlb, lb, mclxy
      parameter (npnt = 12345)
      dimension mclxy(2,npnt)
      common /xy/ xx(npnt),yy(npnt)
c     common /lb/ nlb(5,5),lb(18296,5,5) /boxl/bxlx,bxly,dens,np,idpat
      common /lb/ nlb(9,9),lb(5555,9,9) /boxl/bxlx,bxly,dens,np,idpat
      common /cnst/bx2,by2, pi,pi2,pi5,eps, bxlw,bxup,bylw,byup
c
c
      do 10 k=1,ncly
      do 10 j=1,nclx
   10 nlb(j,k)=0
c
c
      pclx=nclx
      pcly=ncly
      ulx=pclx/bxlx
      uly=pcly/bxly
      do 1000 n=1,np
        fx = (xx(n) - bxlw) * ulx
        fy = (yy(n) - bylw) * uly
        j=fx+1
c       if(j.lt.1.or.j.gt.nclx)  stop 10101
        k=fy+1
c       if(k.lt.1.or.k.gt.ncly)  stop 10101
        nlb(j,k)=nlb(j,k)+1
        mclxy(1,n)=j
        mclxy(2,n)=k
        im=nlb(j,k)
        if(im.gt.5555) then
          write(6,60) np,idpat,bxlx,bxly, j,k, im
   60 format(1h0,10x,'*** sub blking ***   np=',i6,3x,'idpat=',i5,
     1       5x,'bxlx,bxly=',2f10.5,/
     2      1h ,10x,'##### stop 1012 #####  nlb(',2i3,' )=',i6/)
          stop 1012
        end if
        lb(im,j,k)=n
 1000 continue
      write(6,66) ((i1,i2,nlb(i1,i2),i1=1,nclx),i2=1,ncly)
   66 format(1h /(1h ,3x,'(',2i5,')',i8))
      return
      end
      subroutine innova (l1,l2,l3, inbr, mpair, ndt,t1,t2)
c
c   **** innovation of the array for the atoms <l1>, <l2> and <l3> ****
c          arrays to be registered :
c             nbr
c            mnbr
c            inbr
c           mpair
c             ndt
c          -----------------------
c
      real*8 t1,t2,bxlx,bxly,dens
      integer*4 mnbr,inbr,mpair
      parameter (nvmx = 45)
      parameter (npnt = 12345)
      common /boxl/bxlx,bxly,dens,np
      common /mnbr/ mnbr(nvmx,npnt),nbr(npnt)
      dimension inbr(nvmx,npnt), mpair(2,nvmx,npnt), ndt(npnt)
c
c
      if(l1.le.0.or.l1.gt.np) go to 10
      call innov1 (l1,l2,l3, inbr, mpair, ndt,t1,t2)
   10 continue
      if(l2.le.0.or.l2.gt.np)go to 20
      call innov1 (l2,l3,l1, inbr, mpair, ndt,t1,t2)
   20 continue
      if(l3.le.0.or.l3.gt.np)return
      call innov1 (l3,l1,l2, inbr, mpair, ndt,t1,t2)
      return
      end
*
      subroutine innov1 (l1,l2,l3, inbr,mpair,ndt,t1,t2)
c
c      ******************************************************
c      *                                                    *
c      *  innovation of arrays for the atom <l1>            *
c      *        <l2> and <l3> are contiguous to <l1>.       *
c      *                                                    *
c      ******************************************************
c
      implicit real*8 (a-h,o-z)
      integer*4 mnbr,inbr,mpair
      parameter (nvmx = 45)
      parameter (npnt = 12345)
*
      dimension inbr(nvmx,npnt), mpair(2,nvmx,npnt), ndt(npnt)
      common /xy/ xx(npnt),yy(npnt) /boxl/bxlx,bxly,dens,np,idpat
      common /cnst/bx2,by2
      common /iperio/ iperio
      common  /vtxy/ vtx(nvmx,npnt),vty(nvmx,npnt)
      common /mnbr/ mnbr(nvmx,npnt),nbr(npnt)
c
      nbr1=nbr(l1)
      if(nbr1.le.0) go to 5
      do 100 k1=1,nbr1
        if(mnbr(k1,l1).ne.l2) go to 100
        inbr(k1,l1)=inbr(k1,l1)+1
        go to 110
  100 continue
      go to 5
  110 if(nbr1.le.0) go to 6
      do 120 k1=1,nbr1
        if(mnbr(k1,l1).ne.l3) go to 120
        inbr(k1,l1)=inbr(k1,l1)+1
        go to 130
  120 continue
      go to 6
    5 nbr(l1)=nbr(l1)+1
      mnbr(nbr(l1),l1)=l2
      nbr2=nbr(l1)
      inbr(nbr2,l1)=inbr(nbr2,l1)+1
      go to 110
    6 nbr(l1)=nbr(l1)+1
      mnbr(nbr(l1),l1)=l3
      nbr2=nbr(l1)
      inbr(nbr2,l1)=inbr(nbr2,l1)+1
  130 continue
      nbr1=nbr(l1)
      nbr2=nbr(l2)
c
      if(nbr1.gt.50) then
        write(6,60) l1,l2,l3,nbr1,nbr2
   60   format(1h ,' *** error in innov1 ***'/
     1         1h ,' l1,l2,l3,nbr1,nbr2 =',5i8)
        stop 1555
      end if
      ndt(l1)=ndt(l1)+1
      mpair(1,ndt(l1),l1)=l2
      mpair(2,ndt(l1),l1)=l3
      kdt=ndt(l1)
      if(iperio.eq.1) go to 140
      vtx(kdt,l1)=t1
      vty(kdt,l1)=t2
      return
  140 continue
      tx=t1
      ty=t2
      dx=xx(l1)-tx
      dy=yy(l1)-ty
      if(dabs(dx).gt.bx2) tx=tx+dsign(bxlx,dx)
      if(dabs(dy).gt.by2) ty=ty+dsign(bxly,dy)
      vtx(kdt,l1)=tx
      vty(kdt,l1)=ty
      return
      end
*
      subroutine qsort (ii,jj,dm,iord)
c
c     ***********************************************
c     *                                             *
c     *     quick sort                              *
c     *                                             *
c     *       sorting of array dm(k) , k=ii,jj      *
c     *                                             *
c     *       result :  iord(k),k=ii,jj             *
c     *           array <iord(k)> is arranged       *
c     *           in ascending order of values of   *
c     *           <dm>.                             *
c     *                                             *
c     ***********************************************
c
c
      implicit real*8 (a-h,o-z)
      dimension dm(1),iord(1),iu(16),il(16)
c
      m=1
      i=ii
      j=jj
    5 if(i.ge.j) go to 70
   10 k=i
      ij=(j+i)/2
      a=dm(iord(ij))
      it=iord(ij)
      if(dm(iord(i)).le.a) go to 20
      iord(ij)=iord(i)
      a=dm(iord(ij))
      iord(i)=it
      it=iord(ij)
   20 l=j
      if(dm(iord(j)).ge.a) go to 40
      iord(ij)=iord(j)
      iord(j)=it
      a=dm(iord(ij))
      it=iord(ij)
      if(dm(iord(i)).le.a) go to 40
      iord(ij)=iord(i)
      iord(i)=it
      a=dm(iord(ij))
      it=iord(ij)
      go to 40
   30 iord(l)=iord(k)
      iord(k)=itt
   40 l=l-1
      if(dm(iord(l)).gt.a) go to 40
      itt=iord(l)
   50 k=k+1
      if(dm(iord(k)).lt.a) go to 50
      if(k.le.l) go to 30
      if(l-i.le.j-k) go to 60
      il(m)=i
      iu(m)=l
      i=k
      m=m+1
      go to 80
   60 il(m)=k
      iu(m)=j
      j=l
      m=m+1
      go to 80
   70 m=m-1
      if(m.eq.0) return
      i=il(m)
      j=iu(m)
   80 if(j-i.ge.ii) go to 10
      if(i.eq.ii) go to 5
      i=i-1
   90 i=i+1
      if(i.eq.j) go to 70
      a=dm(iord(i+1))
      it=iord(i+1)
      if(dm(iord(i)).le.a) go to 90
      k=i
  100 iord(k+1)=iord(k)
      k=k-1
      if(a.lt.dm(iord(k))) go to 100
      iord(k+1)=it
      go to 90
      end
       subroutine period (xd,yd)
c
c     *******************************************
c     *                                         *
c     *      periodic boundary condition        *
c     *     xd, yd : distance on the torus      *
c     *         smaller distance is selected    *
c     *                                         *
c     *******************************************
c
      implicit real*8 (a-h,o-z)
      common /cnst/ bx2,by2     /boxl/ bxlx, bxly   /iperio/ iperio
c
      if(iperio.ne.1) return
      if(dabs(xd).gt.bx2) xd=xd-dsign(bxlx,xd)
      if(dabs(yd).gt.by2) yd=yd-dsign(bxly,yd)
      return
      end
*
      subroutine mirror (im, xc,yc, xd,yd)
c
c     *******************************************
c     *                                         *
c     *   mirror image construction             *
c     *      centre atom  (xc, yc)              *
c     *      image atom   (xd, yd)              *
c     *          : with respect to the          *
c     *            boundary <im>                *
c     *                                         *
c     *******************************************
c
c
      implicit real*8 (a-h,o-z)
      common /cnst/bx2,by2, pi,pi2,pi5,eps, bxlw,bxup,bylw,byup
c
      xd=0.0d0
      yd=0.0d0
      if(im.eq.1) xd=2.0d0*(bxlw-xc)
      if(im.eq.2) xd=2.0d0*(bxup-xc)
      if(im.eq.3) yd=2.0d0*(bylw-yc)
      if(im.eq.4) yd=2.0d0*(byup-yc)
      return
      end
*
      subroutine out9(np, ii,icnt,inum,jcnt,kcnt)
      implicit double precision (a-h, o-z)
*
      parameter (ndln = 34567)
      dimension icnt(2),jcnt(3,ndln),ip(3)
      dimension kcnt(ndln)
*
      isw=0
      ip(1)=ii
      do 10 i=1,2
        if(icnt(i) .gt. np) then
          ip(i+1) = np - icnt(i)
        else
          ip(i+1) = icnt(i)
        end if
*     if(icnt(i).lt.0) isw=1
   10 continue
*
      call orderk(3, ip)
*
      if(isw.eq.1.or.inum.eq.0) go to 40
      do 50 i=1,inum
        do 60 j=1,3
          if(ip(j).ne.jcnt(j,i)) go to 50
   60   continue
        kcnt(i) = kcnt(i) + 1
        go to 80
   50 continue
   40 continue
      inum=inum+1
      do 70 i=1,3
        jcnt(i,inum)=ip(i)
   70 continue
c 99/05/24
c     ts=1.
c     write(9,1) inum,ii,(icnt(i),i=1,2),ts
    1 format(i7,3x,3i5,1pd15.7)
   80 return
      end
*
      subroutine orderk(k,ind)
      implicit logical (a-z)
      integer k,ind(k)
c
c     purpose: order k elements of array ind in nondecreasing order.
c        it is assume that k is small, say <= 15, so that insertion sort
c        is used. if k is larger, a faster sort such as heapsort should
c        be used.
c
c     input parameters:
c        k - size of array ind
c        ind(1:k) - k integers
c
c     output parameters:
c        ind(1:k) - sorted in nondecreasing order
c
      integer i,j,s,t
c
      do 20 i = 2,k
         t = ind(i)
         j = i
   10    continue
         s = ind(j-1)
         if (t .lt. s) then
            ind(j) = s
            j = j - 1
            if (j .gt. 1) go to 10
         endif
         ind(j) = t
   20 continue
      end
      subroutine xydchk(xx,yy,np)
      implicit real*8 (a-h,o-p,r-z)
      parameter (npnt = 12345)
      dimension xx(npnt),yy(npnt)
      ii=0
      do 100 i=1,np
      xx1=xx(i)
      yy1=yy(i)
      if(i.eq.1) go to 300
      do 200 j=1,ii
      if(abs(xx1-xx(j)).gt.1.d-5) go to 300
      if(abs(yy1-yy(j)).gt.1.d-5) go to 300
      write(6,*) 'data  chk',i,j,xx1,yy1
      go to 100
  200 continue
  300 ii=ii+1
      xx(ii)=xx1
      yy(ii)=yy1
  100 continue
      np=ii
      return
      end
