program evalpoints

    use iso_c_binding, only: c_loc
    use optvormod2, only: pdata_type,compallJ

    implicit none

    ! LOCAL SCALARS 
    integer :: allocerr,nsites,n,j,flag,maxnv,p !,k 
    real(kind=8) :: c,scale,area,marg,auxper
    logical :: debug
    type(pdata_type), target :: pdata


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

   

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

    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, 1.0d0, 1.0d0, 0.0d0, 0.0d0 /)
    pdata%ypol(1,1:pdata%nvpols(1)) = (/ 0.0d0, 0.0d0, 1.0d0, 1.0d0, 0.0d0 /)

    pdata%poledges(1,1:pdata%nvpols(1)-1) = (/ 0, 0, 0, 0 /) ! 0 means frontier of A, 1 means internarl

    scale = sqrt( real(nsites) )
    area = scale ** 2
    write(*,*) 'area = ',area
    write(*,*) 'scale = ',scale
    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)) )

    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(1:2) - (/ pdata%xmin, pdata%ymin /) )
    pdata%drawscale = 1.0d0/scale

    n = 2 * nsites

    p = 1001
    ! Goal function data
    pdata%Jweight(0) = 1.0d0
    pdata%Jweight(1:3) = 0.0d0
    ! pdata%Jweight(1) = 1.0d0
    pdata%Jweight(2) = 1.0d0
    ! pdata%Jweight(3) = 1.0d0

    allocate(x(n), f(0:3,n/2), percentage(p), stat=allocerr)
    if ( allocerr .ne. 0 ) then
        write(*,*) 'Allocation error'
        stop
    end if

    ! open(unit=10,file="/Users/juansebastian/Desktop/Mestrado USP/Dissertation/optlloyd/xoptimal-bfgs-1000.txt",status="old")
    ! open(unit=10,file="/Users/juansebastian/Desktop/Mestrado USP/Dissertation/optlloyd/xoptimal-bfgsJ2-1000-c-0.10.txt",status="old")
    ! open(unit=10,file="/Users/juansebastian/Desktop/Mestrado USP/Dissertation/optlloyd/xoptimal-bfgsJ2-1000-c-0.20.txt",status="old")
    ! open(unit=10,file="/Users/juansebastian/Desktop/Mestrado USP/Dissertation/optlloyd/xoptimal-bfgsJ2-1000-c-0.30.txt",status="old")
    ! open(unit=10,file="/Users/juansebastian/Desktop/Mestrado USP/Dissertation/optlloyd/xoptimal-bfgsJ2-1000-c-0.40.txt",status="old")
    ! open(unit=10,file="/Users/juansebastian/Desktop/Mestrado USP/Dissertation/optlloyd/xoptimal-bfgsJ2-1000-c-0.50.txt",status="old")
    ! open(unit=10,file="/Users/juansebastian/Desktop/Mestrado USP/Dissertation/optlloyd/xoptimal-bfgsJ2-1000-c-0.60.txt",status="old")
    ! open(unit=10,file="/Users/juansebastian/Desktop/Mestrado USP/Dissertation/optlloyd/xoptimal-bfgsJ2-1000-c-0.70.txt",status="old")
    ! open(unit=10,file="/Users/juansebastian/Desktop/Mestrado USP/Dissertation/optlloyd/xoptimal-bfgsJ2-1000-c-0.80.txt",status="old")
    open(unit=10,file="/Users/juansebastian/Desktop/Mestrado USP/Dissertation/optlloyd/xoptimal-bfgsJ2-1000-c-0.90.txt",status="old")


    ! open(unit=10,file="/Users/juansebastian/Desktop/Mestrado USP/Dissertation/optlloyd/xoptimal-bfgsJ3-1000-c-0.50.txt",status="old")
    ! open(unit=10,file="/Users/juansebastian/Desktop/Mestrado USP/Dissertation/optlloyd/xoptimal-bfgsJ3-1000-c-0.60.txt",status="old")
    ! open(unit=10,file="/Users/juansebastian/Desktop/Mestrado USP/Dissertation/optlloyd/xoptimal-bfgsJ3-1000-c-0.70.txt",status="old")
    ! open(unit=10,file="/Users/juansebastian/Desktop/Mestrado USP/Dissertation/optlloyd/xoptimal-bfgsJ3-1000-c-0.80.txt",status="old")

    do j=1,n
        read(10,"(D24.12)") x(j)
    end do
    close(10)

    
    do j=1,p
        if (pdata%Jweight(2) .gt. 0.0d0 ) then

            pdata%J2ctol = real( j-1 )/real( p-1 )
        
            debug = .false.
            call compallJ(n,x,f,debug,flag,c_loc(pdata))
            if ( flag .ne.  0) then
                write(*,*) 'compallJ returned a non-null flag!'
            end if

            auxper = real( count( f(1,1:n/2) .gt. 1.0d-04 ) )
            percentage(j) = 1.0d0 - auxper/real(nsites)
        else if ( pdata%Jweight(3) .gt. 0.0d0 ) then

            pdata%J3ctol = real( j-1 )/real( p-1 )
        
            debug = .false.
            call compallJ(n,x,f,debug,flag,c_loc(pdata))
            if ( flag .ne.  0) then
                write(*,*) 'compallJ returned a non-null flag!'
            end if

            auxper = real( count( f(2,1:n/2) .gt. 1.0d-04 ) )
            percentage(j) = 1.0d0 - auxper/real(nsites)
        end if
    end do

    if ( pdata%Jweight(2) .gt. 0.0d0 ) then

        write(*,*) ' Write the value of c_2 used to obtain a^*: '
        read(*,*) c

        write(filename,"(A14,I0,A3,F4.2,A4)") 'percentagesJ2-',nsites,'-c-',c,'.txt'
        open(unit=20,file=filename)
        do j=1,p
            write(20,"(D24.12)") percentage(j)
        end do
        close(20)
    else if ( pdata%Jweight(3) .gt. 0.0d0 ) then

        write(*,*) ' Write the value of c_3 used to obtain a^*: '
        read(*,*) c

        write(filename,"(A14,I0,A3,F4.2,A4)") 'percentagesJ3-',nsites,'-c-',c,'.txt'
        open(unit=20,file=filename)
        do j=1,p
            write(20,"(D24.12)") percentage(j)
        end do
        close(20)
    end if

    deallocate(x,f,percentage,stat=allocerr)
    if ( allocerr .ne. 0 ) then
        write(*,*) 'Deallocation error.'
        stop
    end if
    stop
        
end program evalpoints