module voro

  use VorCells_Polygons, only: Polygon,VoronoiCell,voronoi_cell,VorCellInterConvPol, &
       VorCellInterConvPol_Collinear,voronoi_cell_destroy,polygon_destroy

  implicit none
  
  type :: pdata_type
     integer :: npols
     real(kind=8) :: c(2),rmax,volA,xmin,xmax,ymin,ymax,drawscale
     integer, allocatable :: nvpols(:),poledges(:,:)
     real(kind=8), allocatable :: xpol(:,:),ypol(:,:)
  end type pdata_type
  
  public :: voronoi,drawvor

contains

  ! ******************************************************************
  ! ******************************************************************

  subroutine voronoi(nsites,sites,An,Ax,Ay,Aflag,nvmax,sstart,nv,vx,vy,vflag,istop)
	
    implicit none
  
    ! SCALAR ARGUMENTS
    integer, intent(in) :: An,nsites,nvmax
    integer, intent(out) :: istop,nv
    
    ! ARRAY ARGUMENTS
    integer, intent(in) :: Aflag(An+1)
    real(kind=8), intent(in) :: Ax(An+1),Ay(An+1),sites(2,nsites)
    integer, intent(out) :: sstart(nsites+1),vflag(nvmax)
    real(kind=8), intent(out) :: vx(nvmax),vy(nvmax)
    
    ! LOCAL SCALARS
    integer :: allocerr,i,ierror,ind,j,k,status,ntriag
    type(VoronoiCell) :: Vi
    type(Polygon) :: Wi,P
    
    ! LOCAL ARRAYS
    integer :: indices(nsites),sitetv(3*2*nsites),start(nsites+1),til(3,2*nsites),tnbr(3,2*nsites)
    real(kind=8) :: vorxy(2,2*nsites)
    
   istop = 0

    if ( nsites .lt. 3 ) then
      write(*,*) 'In voronoi: nsites must be at least 3!'
      stop
    end if

    ! Polygon representing the domain

    P%n = An + 1
    P%deg = .false.

    allocate(P%v(2,P%n),P%e(P%n-1),stat=allocerr)
    if ( allocerr .ne. 0 ) then
       write(*,*) 'In voronoi: Allocation error.'
       stop
    end if
    
    P%v(1,1:P%n) = Ax(1:An+1)
    P%v(2,1:P%n) = Ay(1:An+1)
    P%e(1:P%n-1) = Aflag(1:An)
          
    ! Compute Delaunay triangulation

    indices(1:nsites) = (/ (i,i=1,nsites) /)

    call dtris2(nsites,sites,indices,ntriag,til,tnbr,ierror)

    if ( ierror .ne. 0 ) then
      istop = -1
      return
    end if
    
    if ( ierror .ne. 0 .and. ierror .ne. 225 ) then
       write(*,*) 'In voronoi: dtris2 returned an error different from 0 and 225: ',ierror
       istop = -1
    end if
  
    call voronoi_cell_patch(nsites,ntriag,til,start,sitetv)
    
    ! Compute Voronoi diagram
          
    if ( ntriag .gt. 2 * nsites ) then
       write(*,*) 'In voronoi: There is something wrong. According to the dtris2, '
       write(*,*) 'documentation the number of triangles in the Delaunay trangulation '
       write(*,*) 'cannot be larger than twice the number of sites!'
       stop
    end if
  
    do k = 1,ntriag
       call triangle_circumcenter_2d(sites(1:2,til(1:3,k)),vorxy(1:2,k))
       if ( isnan(vorxy(1,k)) .or. isnan(vorxy(1,k)) ) then
          write(*,*) 'In voronoi: triangle_circumcenter_2d returned NaN.'
          write(*,*) 'Triangle vertices: ',sites(1:2,til(1:3,k))
          istop = -2
       end if
    end do

    sstart(1) = 1

    do i = 1,nsites
       call voronoi_cell(i,nsites,sites,start,sitetv,ntriag,til,tnbr,vorxy,Vi,status)
       if ( status .ne. 0 ) then
          write(*,*) 'In voronoi: Error in voronoi_cell, status = ',status
          stop
       end if

       ! Compute the intersection of the ith Voronoi cell and the domain

       if ( nsites .ge. 3 .and. ierror .eq. 0 ) then
          call VorCellInterConvPol(sites(1:2,i),P,Vi,Wi)
       else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then
          call VorCellInterConvPol_Collinear(i,nsites,sites,P,Wi)
       end if

       if ( Wi%deg ) then
          sstart(i+1) = sstart(i)

       else
          sstart(i+1) = sstart(i) + ( Wi%n - 1 )
          
          if ( sstart(i+1) .gt. nvmax + 1 ) then
             write(*,*) 'Increase nvmax.'
             stop
          end if
    
          do j = 1,Wi%n - 1
             ind = sstart(i) + j - 1
             vx(ind) = Wi%v(1,j)
             vy(ind) = Wi%v(2,j)
             vflag(ind) = -Wi%e(j)
          end do
       end if
       
       call polygon_destroy(Wi)

       call voronoi_cell_destroy(Vi)
    end do

    nv = sstart(nsites+1) - 1
    
    deallocate(P%v,P%e,stat=allocerr)
    if ( allocerr .ne. 0 ) then
       write(*,*) 'In voronoi: Deallocation error.'
       stop
    end if

  end subroutine voronoi

  ! ******************************************************************
  ! ******************************************************************

  subroutine drawvor(nsites,sites,An,Ax,Ay,sstart,nv,vx,vy)
	
    implicit none
  
    ! SCALAR ARGUMENTS
    integer, intent(in) :: An,nsites,nv
    
    ! ARRAY ARGUMENTS
    real(kind=8), intent(in) :: Ax(An+1),Ay(An+1),sites(2,nsites),vx(nv),vy(nv)
    integer, intent(in) :: sstart(nsites+1)
    
    ! LOCAL SCALARS
    integer :: i,j,jnext,k
    
    open(unit=10,file='voronoi.mp')                           

    write(10,10)

    do i = 1,nsites
       ! Draw the cell (intersected with the domain)
       do j = sstart(i),sstart(i+1) - 1
          if ( j .lt. sstart(i+1) - 1 ) then
             jnext = j + 1
          else
             jnext = sstart(i)
          end if
          write(10,20) vx(j),vy(j),vx(jnext),vy(jnext),'desire'
       end do

       ! Draw the ith site
       write(10,30) sites(1:2,i)
       write(10,31) i,sites(1:2,i)
    end do
      
    ! Draw the domain
    do k = 1,An - 1
       write(10,20) Ax(k),Ay(k),Ax(k+1),Ay(k+1),'black'
    end do

    ! Close figure
    
    write(10,40)
    
    close(10)

    ! NON-EXECUTABLE STATEMENTS
    
10  format('prologues := 3;',/, &
         'outputtemplate := "%j.mps";',/, &
         'input mpcolornames;',/, &
         'beginfig( 1 );',/, &
         'color skyblue, dollarbill, desire, jonquil, royalblue, teagreen, ', &
         'orange, plum, bluea, blueb, bluec, blued, bluee, corn, deepsaffron, ferrarired;',/, &
         'skyblue      = (135/255,206/255,235/255);',/, &
         'dollarbill   = (133/255,187/255,101/255);',/, &
         'teagreen     = (208/255,240/255,192/255);',/, &
         'desire       = (234/255, 60/255, 83/255);',/, &
         'jonquil      = (244/255,202/255, 22/255);',/, &
         'royalblue    = ( 65/255,105/255,225/255);',/, &
         'orange       = (255/255,165/255,  0/255);',/, &
         'plum         = (221/255,160/255,221/255);',/, &
         'bluea        = (  0/255,230/255,230/255);',/, &
         'blueb        = (  0/255,191/255,230/255);',/, &
         'bluec        = (  0/255,153/255,230/255);',/, &
         'blued        = (  0/255,115/255,230/255);',/, &
         'bluee        = (  0/255, 77/255,230/255);',/, &
         'corn         = (255/255,230/255,102/255);',/, &
         'deepsaffron  = (255/255,153/255, 51/255);',/, &
         'ferrarired   = (255/255, 42/255,  0/255);',/, &
         'u = 10cm;')
20  format('draw ( ',f20.10,'u,',f20.10,'u)--(',f20.10,'u,',f20.10,'u) withcolor',a12,';')
30  format('drawdot(',f20.10,'u,',f20.10,'u);')
31  format('dotlabel.bot(btex ',i2,' etex, (',f20.10,'u,',f20.10,'u));')
40  format('endfig;',/,'end;')

  end subroutine drawvor
  
end module voro
