program optvorma

  use iso_c_binding, only: c_loc
  use optvor, only: pdata_type,evalfg,approxg,drawsol,compscale
  
  implicit none
	
  ! LOCAL SCALARS
  integer :: allocerr,i,inform,j,maxnv,n,nsites
  real(kind=8) :: area,f,marg,scale
  type(pdata_type), target :: pdata

  ! LOCAL ARRAYS
  real(kind=8), allocatable :: g(:),gtmp(:),x(:)
  character(len=80) :: filename

  ! Problem data
  
  pdata%Jweight(0:5) = 0.0d0
  pdata%Jweight(4) = 1.0d0

  pdata%J3ctol = 0.95d0
  pdata%J4ctol = 0.95d0
  
  ! Set sites

  nsites = 3

  ! Set disjoint convex polygons Aj whose union define A 
  
  pdata%npols = 1
  
  allocate(pdata%nvpols(pdata%npols),stat=allocerr)
  if ( allocerr .ne. 0 ) then
     write(*,*) 'Allocation error.'
     stop
  end if

  pdata%nvpols(1:pdata%npols) = 5

  maxnv = maxval( pdata%nvpols(1:pdata%npols) )

  allocate(pdata%poledges(pdata%npols,maxnv),pdata%xpol(pdata%npols,maxnv), &
       pdata%ypol(pdata%npols,maxnv),stat=allocerr)
  if ( allocerr .ne. 0 ) then
     write(*,*) 'Allocation error.'
     stop
  end if

  pdata%xpol(1,1:pdata%nvpols(1))  = (/ 0.0d0, 2.0d0, 2.0d0, 0.0d0, 0.0d0 /)
  pdata%ypol(1,1:pdata%nvpols(1))  = (/ 0.0d0, 0.0d0, 3.0d0, 3.0d0, 0.0d0 /)
  
  pdata%poledges( 1,1:pdata%nvpols(1)-1)  = (/ 0, 0, 0, 0 /) ! 0 means frontier of A, 1 means internal

  ! Scale A so cells have ideal area near unity
  
  call compscale(nsites,scale,area,c_loc(pdata))

  pdata%volA = area
  
  do j = 1,pdata%npols
     pdata%xpol(j,1:pdata%nvpols(j)) = scale * pdata%xpol(j,1:pdata%nvpols(j))
     pdata%ypol(j,1:pdata%nvpols(j)) = scale * pdata%ypol(j,1:pdata%nvpols(j))
  end do
  
  pdata%xmin = minval( pdata%xpol(1,1:pdata%nvpols(1)) )
  pdata%xmax = maxval( pdata%xpol(1,1:pdata%nvpols(1)) )
  pdata%ymin = minval( pdata%ypol(1,1:pdata%nvpols(1)) )
  pdata%ymax = maxval( pdata%ypol(1,1:pdata%nvpols(1)) )
  do j = 2,pdata%npols
     pdata%xmin = min( pdata%xmin, minval( pdata%xpol(j,1:pdata%nvpols(j)) ) )
     pdata%xmax = max( pdata%xmax, maxval( pdata%xpol(j,1:pdata%nvpols(j)) ) )
     pdata%ymin = min( pdata%ymin, minval( pdata%ypol(j,1:pdata%nvpols(j)) ) )
     pdata%ymax = max( pdata%ymax, maxval( pdata%ypol(j,1:pdata%nvpols(j)) ) )
  end do
  
  marg = 0.1d0 * max( 1.0d0, abs( pdata%xmin ), abs( pdata%xmax ), abs( pdata%ymin), abs( pdata%ymax ) )

  pdata%xmin = pdata%xmin - marg
  pdata%xmax = pdata%xmax + marg
  pdata%ymin = pdata%ymin - marg
  pdata%ymax = pdata%ymax + marg
  
  pdata%c(1:2) = 0.5d0 * (/ pdata%xmax + pdata%xmin, pdata%ymax + pdata%ymin /)
  pdata%rmax = norm2(pdata%c - (/ pdata%xmin, pdata%ymin /))
  pdata%drawscale = 1.0d0 / scale

  pdata%Jweight(0:5) = 0.0d0
  pdata%Jweight(0)  = 10.0d0
  pdata%Jweight(1)  =  1.0d0
  
  pdata%J0delta = 0.1d0
  
  n = 2 * nsites
  
  allocate(x(n),g(n),gtmp(n),stat=allocerr)
  if ( allocerr .ne. 0 ) then
     write(*,*) 'Allocation error.'
     stop
  end if

!!$  x(1:2) = scale * (/ 1.75d0, 0.05d0 /)
!!$  x(3:4) = scale * (/ 1.75d0, 1.50d0 /)
!!$  x(5:6) = scale * (/ 0.75d0, 1.50d0 /)
!!$
  x(1:2) = scale * (/ 1.0d0, 0.50d0 /)
  x(3:4) = scale * (/ 1.5d0, 1.50d0 /)
  x(5:6) = scale * (/ 0.5d0, 1.50d0 /)

  x(1:2) = scale * (/ 1.0d0, 1.0d0 /)
  x(3:4) = scale * (/ 1.5d0, 2.0d0 /)
  x(5:6) = scale * (/ 0.5d0, 2.0d0 /)

!!$  x(1:2) = scale * (/ 1.0d0, 1.0d0 /)
!!$  x(3:4) = scale * (/ 1.0d0, 1.5d0 /)
!!$  x(5:6) = scale * (/ 1.0d0, 0.5d0 /)
!!$  x(7:8) = scale * (/ 1.5d0, 1.0d0 /)
!!$  x(9:10) = scale * (/ 0.5d0, 1.0d0 /)
  
  call evalfg(n,x,f,g,inform,c_loc(pdata))

  write(*,*) 'f = ',f
  write(*,*) 'g = ',g(1:n)

  call approxg(n,x,gtmp,1.0d-06,inform,c_loc(pdata))
  write(*,*) 'approxg = ',gtmp(1:n)

  do i = 1,n
    write(*,*) i,g(i),gtmp(i),abs(g(i)-gtmp(i))
  end do

  write(*,*) 'Maximum absolute error = ',maxval( abs( g(1:n) - gtmp(1:n) ) )
  
  write(filename,"(A18,I0,A3)") 'fig-2x3-rectangle-',nsites,'.mp'
  call drawsol(n,x,filename,c_loc(pdata))
    
  deallocate(x,g,gtmp,stat=allocerr)
  if ( allocerr .ne. 0 ) then
     write(*,*) 'Deallocation error.'
     stop
  end if
  
  stop
  
end program optvorma