program voroma

  use voro, only: voronoi,drawvor
  
  implicit none

  ! PARAMETER
  integer, parameter :: nvmax = 100
  
  ! LOCAL SCALARS
  integer :: allocerr,An,i,nsites,nv,j,istop
  real(kind=8) :: seed,mindist,dist, aux1, aux2

  ! LOCAL ARRAYS
  integer :: vflag(nvmax)
  integer, allocatable :: Aflag(:),sstart(:)
  real(kind=8) :: vx(nvmax),vy(nvmax)
  real(kind=8), allocatable :: Ax(:),Ay(:),sites(:,:)
    
  ! FUNCTIONS
  real(kind=8), external :: drand
  
  ! Set domain A 

  An = 5
  
  allocate(Ax(An+1),Ay(An+1),Aflag(An+1),stat=allocerr)
  if ( allocerr .ne. 0 ) then
     write(*,*) 'Allocation error.'
     stop
  end if

  Ax(1) = 0.0d0
  Ay(1) = 0.0d0
  Aflag(1) = 4

  Ax(2) = 1.0d0
  Ay(2) = 0.0d0
  Aflag(2) = 2

  Ax(3) = 1.0d0
  Ay(3) = 1.0d0
  Aflag(3) = 3

  Ax(4) = 0.0d0
  Ay(4) = 1.0d0
  Aflag(4) = 1

  Ax(5) = Ax(1)
  Ay(5) = Ay(1)
  Aflag(5) = Aflag(1)

  ! Set sites

  write(*,*) 'Enter nsites>0: '
  read(*,*) nsites

  ! nsites = 3
  
  allocate(sites(2,nsites),sstart(nsites+1),stat=allocerr)
  if ( allocerr .ne. 0 ) then
     write(*,*) 'Allocation error.'
     stop
  end if

  ! sites(1:2,1) = (/ 0.23d0, 0.65d0 /)
  ! sites(1:2,2) = (/ 0.87d0, 0.78d0 /)
  ! sites(1:2,3) = (/ 0.72d0, 0.32d0 /)

  
  ! Sort sites within A
  
   seed = 123456.0d0
  i = 1
  do while (i .le. nsites)
    aux1 = drand(seed)
    sites(1,i) = 0.10d0 + 0.80d0 * aux1
    aux2 = drand(seed) 
    sites(2,i) = 0.10d0 + 0.80d0 * aux2
    mindist = huge(1.0d0)
    do j=1,i-1
      dist = sqrt((sites(1,i)-sites(1,j))**2 + (sites(2,i)-sites(2,j))**2)
      mindist = min(mindist,dist)
    end do
    if( mindist .gt. 1.0d-01 ) i = i + 1
  end do


  do i=1,2
    write(*,*) sites(i,:)
  end do

  call voronoi(nsites,sites,An,Ax,Ay,Aflag,nvmax,sstart,nv,vx,vy,vflag,istop)

  write(*,*)
  write(*,*) 'sstart = ',sstart(1:nsites+1)
  write(*,*) 'vx = ',vx(1:nv)
  write(*,*) 'vy = ',vy(1:nv)
  write(*,*) 'vflag = ',vflag(1:nv)

  call drawvor(nsites,sites,An,Ax,Ay,sstart,nv,vx,vy)
  
  deallocate(Ax,Ay,sites,sstart,stat=allocerr)
  if ( allocerr .ne. 0 ) then
     write(*,*) 'Deallocation error.'
     stop
  end if
  
  stop

end program voroma
