module optvormod2

    use iso_c_binding, only: c_ptr,c_f_pointer
    use VorCells_Polygons, only: Polygon,VoronoiCell,voronoi_cell,voronoi_cell_old,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,J2ctol,J3ctol,J4c(2),J4r
       integer, allocatable :: nvpols(:),poledges(:,:)
       real(kind=8), allocatable :: xpol(:,:),ypol(:,:)
       real(kind=8) :: Jweight(0:4)
       ! 0: CVT
       ! 1: cells with identical volume given by vol(A)/m
       ! 2: each cell with equally sized edges
       ! 3: each cell with equally sized angles
       ! 4: cells with volume proportional to density function f evaluated at their sites

    end type pdata_type
    
    public :: evalfg,approxg,drawsol,drawsolJ1,drawsolJ2
  
  contains
  
   ! ******************************************************************
   ! ******************************************************************

  logical function InA(p,pdataptr)
      
  implicit none

  ! SCALAR ARGUMENTS
  type(c_ptr), optional, intent(in) :: pdataptr
  
  ! ARRAY ARGUMENTS
  real(kind=8), intent(in) :: p(2)

  ! LOCAL SCALARS
  integer :: allocerr,j,maxnv
  type(pdata_type), pointer :: pdata

  ! LOCAL ARRAYS
  real(kind=8), allocatable :: vertices(:,:)
  
  call c_f_pointer(pdataptr,pdata)
    
  maxnv = maxval( pdata%nvpols(1:pdata%npols) )
    
  allocate(vertices(2,maxnv),stat=allocerr)
  if ( allocerr .ne. 0 ) then
     write(*,*) 'InA: Allocation error.'
     stop
  end if

  j = 1
  InA = .false.
  do while ( j .le. pdata%npols .and. .not. InA )
     vertices(1,1:pdata%nvpols(j)-1) = pdata%xpol(j,1:pdata%nvpols(j)-1)
     vertices(2,1:pdata%nvpols(j)-1) = pdata%ypol(j,1:pdata%nvpols(j)-1)
     call polygon_contains_point_2d_convex(pdata%nvpols(j)-1,vertices,p,inA)
     j = j + 1
  end do
      
  deallocate(vertices,stat=allocerr)
  if ( allocerr .ne. 0 ) then
     write(*,*) 'InA: Deallocation error.'
     stop
  end if
  
end function InA
   
! ******************************************************************
! ******************************************************************

subroutine approxg(n,x,g,h,flag,pdataptr)
      
      implicit none
    
      ! SCALAR ARGUMENTS
      integer, intent(in) :: n
      integer, intent(out) :: flag
      real(kind=8), intent(in) :: h
      type(c_ptr), optional, intent(in) :: pdataptr
        
      ! ARRAY ARGUMENTS
      real(kind=8), intent(inout) :: x(n)
      real(kind=8), intent(out) :: g(n)
  
      ! LOCAL SCALARS
      integer :: i
      real(kind=8) :: fplus,fminus,tmp
  
      ! LOCAL ARRAYS
      real(kind=8) :: gtmp(n)
  
      do i = 1,n
         tmp = x(i)
         x(i) = tmp + h
         call evalfg(n,x,fplus,gtmp,flag,pdataptr)
         x(i) = tmp - h
         call evalfg(n,x,fminus,gtmp,flag,pdataptr)
         g(i) = ( fplus - fminus ) / ( 2.0d0 * h )
         x(i) = tmp
      end do
        
   end subroutine approxg
      
   ! ******************************************************************
   ! ******************************************************************
  
   subroutine evalfg(n,x,f,g,flag,pdataptr)
      
      implicit none
    
      ! SCALAR ARGUMENTS
      integer, intent(in) :: n
      integer, intent(out) :: flag
      real(kind=8), intent(out) :: f
      type(c_ptr), optional, intent(in) :: pdataptr
        
      ! ARRAY ARGUMENTS
      real(kind=8), intent(in) :: x(n)
      real(kind=8), intent(out) :: g(n)
  
      ! PARAMETERS
      logical, parameter :: compgrad = .true., debug = .false.
      integer, parameter :: gtmpnnzmax = 100000
      
      ! LOCAL SCALARS
      integer :: allocerr,gtmpnnz,i,ierror,k,maxnv,nsites,ntriag,status
      real(kind=8) :: ftmp
      type(pdata_type), pointer :: pdata
      type(VoronoiCell) :: V
      type(Polygon) :: P,W
      
      ! LOCAL ARRAYS
      integer :: gtmpind(gtmpnnzmax),indices(n/2),sitetv(3*2*(n/2)),start((n/2)+1),til(3,2*(n/2)),tnbr(3,2*(n/2))
      real(kind=8) :: fden(n/2),gden(2,n/2),gtmpval(gtmpnnzmax),sites(2,n/2),vorxy(2,2*(n/2)),valf(n/2)
          
      flag = 0
      
      call c_f_pointer(pdataptr,pdata)
  
      nsites = n / 2
      indices(1:nsites) = (/ (i,i=1,nsites) /)
      sites(1:2,1:nsites) = reshape( x(1:n), (/ 2, nsites /) )
          
      if ( nsites .lt. 3 ) then
         write(*,*) 'In evalfg: nsites must be at least 3!'
         flag = - 1
         return
      end if
  
      ! Compute Delaunay triangulation
  
      call dtris2(nsites,sites,indices,ntriag,til,tnbr,ierror)
      
      if ( ierror .ne. 0 .and. ierror .ne. 225 ) then
         write(*,*) 'In evalfg: dtris2 returned an error different from 0 and 225: ',ierror
         flag = - 1
         return
      end if
  
      call voronoi_cell_patch(nsites,ntriag,til,start,sitetv)
      
      ! Compute vertices of the Voronoi diagram
            
      if ( ntriag .gt. 2 * nsites ) then
         write(*,*) 'In evalfg: 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!'
         flag = - 1
         return
      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 evalfg: triangle_circumcenter_2d returned NaN.'
            flag = - 1
            return
         end if
      end do

      ! Compute density function at sites if required
      
      if ( pdata%Jweight(4) .gt. 0.0d0 ) then
         do i = 1,nsites
            call fgcompJ4(sites(1:2,i),fden(i),gden(1:2,i),compgrad,pdataptr)
         end do
      end if

      f = 0.0d0
  
      if ( compgrad ) g(1:n) = 0.0d0
      
      maxnv = maxval( pdata%nvpols(1:pdata%npols) )

      allocate(P%v(2,maxnv),P%e(maxnv-1),stat=allocerr)

      if ( allocerr .ne. 0 ) then
         write(*,*) 'In evalfg: Allocation error.'
         flag = - 1
         return
      end if

      P%n = pdata%nvpols(1)
      P%deg = .false.
         
      P%v(1,1:P%n) = pdata%xpol(1,1:P%n)
      P%v(2,1:P%n) = pdata%ypol(1,1:P%n)
      P%e(1:P%n-1) = pdata%poledges(1,1:P%n-1)
  
      do i = 1,nsites
        if ( debug ) then
            write(*,*)
            write(*,*) 'SITE i = ',i
            write(*,*) 'Generator: ',sites(1:2,i)
            write(*,*)
        end if
  
        call voronoi_cell(i,nsites,sites,start,sitetv,ntriag,til,tnbr,vorxy,V,status)
        if ( status .ne. 0 ) then
            write(*,*) 'In evalfg: Error in voronoi_cell, status = ',status
            flag = - 1
            return
        end if


        if ( nsites .ge. 3 .and. ierror .eq. 0 ) then
            call VorCellInterConvPol(sites(1:2,i),P,V,W)
        else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then 
            call VorCellInterConvPol_Collinear(i,nsites,sites,P,W)
        end if
        
        if ( pdata%Jweight(4) .gt. 0.0d0) then
            call compJ4(nsites,sites,i,W,pdata%volA,fden,gden,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug)
            if ( flag .ne. 0) return

            f = f + pdata%Jweight(4) * ftmp ** 2
            if ( compgrad ) g(gtmpind(1:gtmpnnz)) = g(gtmpind(1:gtmpnnz)) + pdata%Jweight(4) * 2.0d0 * ftmp * gtmpval(1:gtmpnnz)
        end if
            
        if ( W%deg ) then
            if ( debug ) then
               write(*,*)
               write(*,*) 'POLYGON Wi with i = ',i, ' IS EMPTY.'
            end if
        
            f = f + 1.0d0

        else
            if ( debug ) then
               write(*,*)
               write(*,*) 'POLYGON Wi with i = ',i
               write(*,*) 'Number of vertices: ',W%n - 1
               write(*,*)
            end if
        
            call compG(nsites,sites,i,W,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug)
            ! call compGalternative(nsites,sites,i,W,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug)
            if ( flag .ne. 0 ) return

            f = f + pdata%Jweight(0)*ftmp   
            
            if ( compgrad ) g(gtmpind(1:gtmpnnz)) = g(gtmpind(1:gtmpnnz)) + pdata%Jweight(0)*gtmpval(1:gtmpnnz)

            if ( pdata%Jweight(1) .gt. 0.0d0) then
               call compJ1(nsites,sites,i,W,pdata%volA,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug)
               if ( flag .ne. 0) return

               f = f + pdata%Jweight(1) * ftmp ** 2
               if ( compgrad ) g(gtmpind(1:gtmpnnz)) = g(gtmpind(1:gtmpnnz)) + pdata%Jweight(1) * 2.0d0 * ftmp * gtmpval(1:gtmpnnz)
            end if
            
            if ( pdata%Jweight(2) .gt. 0.0d0 ) then
               call compJ2(nsites,sites,i,W,pdata%J2ctol,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug)
               if ( flag .ne. 0 ) return

               f = f + pdata%Jweight(2) * ftmp 
               if ( compgrad ) g(gtmpind(1:gtmpnnz)) = g(gtmpind(1:gtmpnnz)) + pdata%Jweight(2) * gtmpval(1:gtmpnnz)
            end if

            if ( pdata%Jweight(3) .gt. 0.0d0 ) then
               call compJ3(nsites,sites,i,W,pdata%J3ctol,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug)
               if ( flag .ne. 0 ) return

               f = f + pdata%Jweight(3) * ftmp 
               if ( compgrad ) g(gtmpind(1:gtmpnnz)) = g(gtmpind(1:gtmpnnz)) + pdata%Jweight(3) * gtmpval(1:gtmpnnz)
            end if

         end if

        call voronoi_cell_destroy(V)
        
        call polygon_destroy(W)

      end do

      f = f / nsites
      
      if ( compgrad ) g(1:n) = g(1:n) / nsites
      
      deallocate(P%v,P%e,stat=allocerr)
      if ( allocerr .ne. 0 ) then
         write(*,*) 'In evalfg: Deallocation error.'
         flag = - 1
         return
      end if
  
    end subroutine evalfg

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

   !! This function works to calculate the values of G, J1, J2, J3

    subroutine evalJg(n,x,G1,j2,j3,j1,j4,ji,g,flag,pdataptr)
      
      implicit none
    
      ! SCALAR ARGUMENTS
      integer, intent(in) :: n
      integer, intent(out) :: flag
      real(kind=8), intent(out) :: j2,G1,j3,j1,j4
      type(c_ptr), optional, intent(in) :: pdataptr
        
      ! ARRAY ARGUMENTS
      real(kind=8), intent(in) :: x(n)
      real(kind=8), intent(out) :: g(n),ji(n/2)
  
      ! PARAMETERS
      logical, parameter :: compgrad = .false., debug = .false.
      integer, parameter :: gtmpnnzmax = 100000
      
      ! LOCAL SCALARS
      integer :: allocerr,gtmpnnz,i,ierror,k,maxnv,nsites,ntriag,status
      real(kind=8) :: ftmp
      type(pdata_type), pointer :: pdata
      type(VoronoiCell) :: V
      type(Polygon) :: P,W
      
      ! LOCAL ARRAYS
      integer :: gtmpind(gtmpnnzmax),indices(n/2),sitetv(3*2*(n/2)),start((n/2)+1),til(3,2*(n/2)),tnbr(3,2*(n/2))
      real(kind=8) :: fden(n/2),gden(2,n/2),gtmpval(gtmpnnzmax),sites(2,n/2),vorxy(2,2*(n/2)),valf(n/2)
          
      flag = 0
      
      call c_f_pointer(pdataptr,pdata)
  
      nsites = n / 2
      indices(1:nsites) = (/ (i,i=1,nsites) /)
      sites(1:2,1:nsites) = reshape( x(1:n), (/ 2, nsites /) )
          
      if ( nsites .lt. 3 ) then
         write(*,*) 'In evalfg: nsites must be at least 3!'
         flag = - 1
         return
      end if
  
      ! Compute Delaunay triangulation
  
      call dtris2(nsites,sites,indices,ntriag,til,tnbr,ierror)
      
      if ( ierror .ne. 0 .and. ierror .ne. 225 ) then
         write(*,*) 'In evalfg: dtris2 returned an error different from 0 and 225: ',ierror
         flag = - 1
         return
      end if
  
      call voronoi_cell_patch(nsites,ntriag,til,start,sitetv)
      
      ! Compute vertices of the Voronoi diagram
            
      if ( ntriag .gt. 2 * nsites ) then
         write(*,*) 'In evalfg: 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!'
         flag = - 1
         return
      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 evalfg: triangle_circumcenter_2d returned NaN.'
            flag = - 1
            return
         end if
      end do

      ! Compute density function at sites if required
   
      do i = 1,nsites
         call fgcompJ4(sites(1:2,i),fden(i),gden(1:2,i),compgrad,pdataptr)
      end do

      
      G1 = 0.0d0
      j1 = 0.0d0
      j2 = 0.0d0
      j3 = 0.0d0
      j4 = 0.0d0

      if ( compgrad ) g(1:n) = 0.0d0
      
      maxnv = maxval( pdata%nvpols(1:pdata%npols) )

      allocate(P%v(2,maxnv),P%e(maxnv-1),stat=allocerr)

      if ( allocerr .ne. 0 ) then
         write(*,*) 'In evalfg: Allocation error.'
         flag = - 1
         return
      end if

      P%n = pdata%nvpols(1)
      P%deg = .false.
         
      P%v(1,1:P%n) = pdata%xpol(1,1:P%n)
      P%v(2,1:P%n) = pdata%ypol(1,1:P%n)
      P%e(1:P%n-1) = pdata%poledges(1,1:P%n-1)
  
      do i = 1,nsites
        if ( debug ) then
            write(*,*)
            write(*,*) 'SITE i = ',i
            write(*,*) 'Generator: ',sites(1:2,i)
            write(*,*)
        end if
  
        call voronoi_cell(i,nsites,sites,start,sitetv,ntriag,til,tnbr,vorxy,V,status)
        if ( status .ne. 0 ) then
            write(*,*) 'In evalfg: Error in voronoi_cell, status = ',status
            flag = - 1
            return
        end if


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

        call compJ4(nsites,sites,i,W,pdata%volA,fden,gden,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug)
        if ( flag .ne. 0 ) return
        ! write(*,*) 'ftmpevalJg = ', ftmp
        j4 = j4 + ftmp ** 2
        if ( compgrad ) g(gtmpind(1:gtmpnnz)) = g(gtmpind(1:gtmpnnz)) + pdata%Jweight(4) * 2.0d0 * ftmp * gtmpval(1:gtmpnnz)
            
        if ( W%deg ) then

            if ( debug ) then
               write(*,*)
               write(*,*) 'POLYGON Wi with i = ',i, ' IS EMPTY.'
            end if

            ! write(*,*) 'Degenerate polygon, do not anything.'
            G1 = G1 + 1.0d0
            j1 = j1 + 1.0d0
            j2 = j2 + 1.0d0
            j3 = j3 + 1.0d0

        else
            if ( debug ) then
               write(*,*)
               write(*,*) 'POLYGON Wi with i = ',i
               write(*,*) 'Number of vertices: ',W%n - 1
               write(*,*)
            end if
        
            call compG(nsites,sites,i,W,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug)
            if ( flag .ne. 0 ) return

            G1 =  G1 + ftmp 
         
            if ( compgrad ) g(gtmpind(1:gtmpnnz)) = g(gtmpind(1:gtmpnnz)) + gtmpval(1:gtmpnnz)

            call compJ1(nsites,sites,i,W,pdata%volA,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug)
            if ( flag .ne. 0 ) return

            j1 = j1 + ftmp ** 2
            ji(i) =  ftmp
            if ( compgrad ) g(gtmpind(1:gtmpnnz)) = g(gtmpind(1:gtmpnnz)) + pdata%Jweight(1) * 2.0d0 * ftmp * gtmpval(1:gtmpnnz)
            

            call compJ2(nsites,sites,i,W,pdata%J2ctol,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug)
            if ( flag .ne. 0 ) return

            j2 = j2 + ftmp
            if ( compgrad ) g(gtmpind(1:gtmpnnz)) = g(gtmpind(1:gtmpnnz)) + pdata%Jweight(1) * gtmpval(1:gtmpnnz)

            call compJ3(nsites,sites,i,W,pdata%J3ctol,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug)
            if ( flag .ne. 0 ) return

            j3 = j3 + ftmp
            if ( compgrad ) g(gtmpind(1:gtmpnnz)) = g(gtmpind(1:gtmpnnz)) + pdata%Jweight(2) * gtmpval(1:gtmpnnz)

         end if

        call voronoi_cell_destroy(V)
        
        call polygon_destroy(W)

      end do

      G1 = G1 / nsites
      j2 = j2 / nsites
      j3 = j3 / nsites
      j1 = j1 / nsites
      j4 = j4 / nsites
      if ( compgrad ) g(1:n) = g(1:n) / nsites
      
      deallocate(P%v,P%e,stat=allocerr)
      if ( allocerr .ne. 0 ) then
         write(*,*) 'In evalfg: Deallocation error.'
         flag = - 1
         return
      end if
  
    end subroutine evalJg

   !! This function works to calculate the volume of each cell, and the J^1_i(a) recalling that
   !! the function J^1(a)=1/nsites sum^{nsites}_{i=1}[J^1_i(a)]^2. 

    subroutine evalvolJi1(n,x,vol,ji1,flag,pdataptr)
      
      implicit none
    
      ! SCALAR ARGUMENTS
      integer, intent(in) :: n
      integer, intent(out) :: flag
      type(c_ptr), optional, intent(in) :: pdataptr
        
      ! ARRAY ARGUMENTS
      real(kind=8), intent(in) :: x(n)
      real(kind=8), intent(out) :: ji1(n/2), vol(n/2)
  
      ! PARAMETERS
      logical, parameter :: compgrad = .false.,  debug = .false.
      integer, parameter :: gtmpnnzmax = 100000
      
      ! LOCAL SCALARS
      integer :: allocerr,i,ierror,k,maxnv,nsites,ntriag,status,ed,ednext,gtmpnnz
      real(kind=8) :: perim,ftmp,volW
      type(pdata_type), pointer :: pdata
      type(VoronoiCell) :: V
      type(Polygon) :: P,W
      
      ! LOCAL ARRAYS
      integer :: gtmpind(gtmpnnzmax),indices(n/2),sitetv(3*2*(n/2)),start((n/2)+1),til(3,2*(n/2)),tnbr(3,2*(n/2))
      real(kind=8) :: sites(2,n/2),vorxy(2,2*(n/2)),vsizeed(10),taued(2),gtmpval(gtmpnnzmax)
          
      flag = 0
      
      call c_f_pointer(pdataptr,pdata)
  
      nsites = n / 2
      indices(1:nsites) = (/ (i,i=1,nsites) /)
      sites(1:2,1:nsites) = reshape( x(1:n), (/ 2, nsites /) )
          
      if ( nsites .lt. 3 ) then
         write(*,*) 'In evalfg: nsites must be at least 3!'
         flag = - 1
         return
      end if
  
      ! Compute Delaunay triangulation
  
      call dtris2(nsites,sites,indices,ntriag,til,tnbr,ierror)
      
      if ( ierror .ne. 0 .and. ierror .ne. 225 ) then
         write(*,*) 'In evalfg: dtris2 returned an error different from 0 and 225: ',ierror
         flag = - 1
         return
      end if
  
      call voronoi_cell_patch(nsites,ntriag,til,start,sitetv)
      
      ! Compute vertices of the Voronoi diagram
            
      if ( ntriag .gt. 2 * nsites ) then
         write(*,*) 'In evalfg: 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!'
         flag = - 1
         return
      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 evalfg: triangle_circumcenter_2d returned NaN.'
            flag = - 1
            return
         end if
      end do

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

      allocate(P%v(2,maxnv),P%e(maxnv-1),stat=allocerr)

      if ( allocerr .ne. 0 ) then
         write(*,*) 'In evalfg: Allocation error.'
         flag = - 1
         return
      end if

      P%n = pdata%nvpols(1)
      P%deg = .false.
         
      P%v(1,1:P%n) = pdata%xpol(1,1:P%n)
      P%v(2,1:P%n) = pdata%ypol(1,1:P%n)
      P%e(1:P%n-1) = pdata%poledges(1,1:P%n-1)

      
  
      do i = 1,nsites
        if ( debug ) then
            write(*,*)
            write(*,*) 'SITE i = ',i
            write(*,*) 'Generator: ',sites(1:2,i)
            write(*,*)
        end if
  
        call voronoi_cell(i,nsites,sites,start,sitetv,ntriag,til,tnbr,vorxy,V,status)
        if ( status .ne. 0 ) then
            write(*,*) 'In evalfg: Error in voronoi_cell, status = ',status
            flag = - 1
            return
        end if


        if ( nsites .ge. 3 .and. ierror .eq. 0 ) then
            call VorCellInterConvPol(sites(1:2,i),P,V,W)
        else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then 
            call VorCellInterConvPol_Collinear(i,nsites,sites,P,W)
        end if
        
            
        if ( W%deg ) then
            if ( debug ) then
               write(*,*)
               write(*,*) 'POLYGON Wi with i = ',i, ' IS EMPTY.'
            end if
            write(*,*) 'Degenerate polygon, do not anything.'
        else
            if ( debug ) then
               write(*,*)
               write(*,*) 'POLYGON Wi with i = ',i
               write(*,*) 'Number of vertices: ',W%n - 1
               write(*,*)
            end if
            call polygon_area_2d(W%n-1,W%v(1:2,1:W%n-1),volW)
            vol(i) = volW

            call compJ1(nsites,sites,i,W,pdata%volA,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug)
            if ( flag .ne. 0 ) return

            ji1(i) = ftmp
        end if

        call voronoi_cell_destroy(V)
        
        call polygon_destroy(W)

      end do
      
      deallocate(P%v,P%e,stat=allocerr)
      if ( allocerr .ne. 0 ) then
         write(*,*) 'In evalfg: Deallocation error.'
         flag = - 1
         return
      end if
  
   end subroutine evalvolJi1

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

   !! This function works to calculate the average of each cell, the size of the edges of each cell and the J^2_i(a) recalling that
   !! the function J^2(a)=1/nsites sum^{nsites}_{i=1}J^2_i(a). 

   subroutine evalsizeedJi2(n,x,avg,sizes,ji2,flag,pdataptr)
      
      implicit none
    
      ! SCALAR ARGUMENTS
      integer, intent(in) :: n
      integer, intent(out) :: flag
      type(c_ptr), optional, intent(in) :: pdataptr
        
      ! ARRAY ARGUMENTS
      real(kind=8), intent(in) :: x(n)
      real(kind=8), intent(out) :: ji2(n/2), avg(n/2), sizes(n/2,10)
  
      ! PARAMETERS
      logical, parameter :: compgrad = .false.,  debug = .false.
      integer, parameter :: gtmpnnzmax = 100000
      
      ! LOCAL SCALARS
      integer :: allocerr,i,ierror,k,maxnv,nsites,ntriag,status,ed,ednext,gtmpnnz
      real(kind=8) :: perim,ftmp
      type(pdata_type), pointer :: pdata
      type(VoronoiCell) :: V
      type(Polygon) :: P,W
      
      ! LOCAL ARRAYS
      integer :: gtmpind(gtmpnnzmax),indices(n/2),sitetv(3*2*(n/2)),start((n/2)+1),til(3,2*(n/2)),tnbr(3,2*(n/2))
      real(kind=8) :: sites(2,n/2),vorxy(2,2*(n/2)),vsizeed(10),taued(2),gtmpval(gtmpnnzmax)
          
      flag = 0
      
      call c_f_pointer(pdataptr,pdata)
  
      nsites = n / 2
      indices(1:nsites) = (/ (i,i=1,nsites) /)
      sites(1:2,1:nsites) = reshape( x(1:n), (/ 2, nsites /) )
          
      if ( nsites .lt. 3 ) then
         write(*,*) 'In evalfg: nsites must be at least 3!'
         flag = - 1
         return
      end if
  
      ! Compute Delaunay triangulation
  
      call dtris2(nsites,sites,indices,ntriag,til,tnbr,ierror)
      
      if ( ierror .ne. 0 .and. ierror .ne. 225 ) then
         write(*,*) 'In evalfg: dtris2 returned an error different from 0 and 225: ',ierror
         flag = - 1
         return
      end if
  
      call voronoi_cell_patch(nsites,ntriag,til,start,sitetv)
      
      ! Compute vertices of the Voronoi diagram
            
      if ( ntriag .gt. 2 * nsites ) then
         write(*,*) 'In evalfg: 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!'
         flag = - 1
         return
      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 evalfg: triangle_circumcenter_2d returned NaN.'
            flag = - 1
            return
         end if
      end do

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

      allocate(P%v(2,maxnv),P%e(maxnv-1),stat=allocerr)

      if ( allocerr .ne. 0 ) then
         write(*,*) 'In evalfg: Allocation error.'
         flag = - 1
         return
      end if

      P%n = pdata%nvpols(1)
      P%deg = .false.
         
      P%v(1,1:P%n) = pdata%xpol(1,1:P%n)
      P%v(2,1:P%n) = pdata%ypol(1,1:P%n)
      P%e(1:P%n-1) = pdata%poledges(1,1:P%n-1)

      
  
      do i = 1,nsites
        if ( debug ) then
            write(*,*)
            write(*,*) 'SITE i = ',i
            write(*,*) 'Generator: ',sites(1:2,i)
            write(*,*)
        end if
  
        call voronoi_cell(i,nsites,sites,start,sitetv,ntriag,til,tnbr,vorxy,V,status)
        if ( status .ne. 0 ) then
            write(*,*) 'In evalfg: Error in voronoi_cell, status = ',status
            flag = - 1
            return
        end if


        if ( nsites .ge. 3 .and. ierror .eq. 0 ) then
            call VorCellInterConvPol(sites(1:2,i),P,V,W)
        else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then 
            call VorCellInterConvPol_Collinear(i,nsites,sites,P,W)
        end if
        
            
        if ( W%deg ) then
            if ( debug ) then
               write(*,*)
               write(*,*) 'POLYGON Wi with i = ',i, ' IS EMPTY.'
            end if

            write(*,*) 'Degenerate polygon, do not anything.'

        else
            if ( debug ) then
               write(*,*)
               write(*,*) 'POLYGON Wi with i = ',i
               write(*,*) 'Number of vertices: ',W%n - 1
               write(*,*)
            end if
            vsizeed(1:10) = 0.0d0
            perim = 0.0d0
            do ed = 1,W%n - 1
               if ( ed .eq. W%n-1) then
                  ednext = 1
               else 
                  ednext = ed + 1
               end if 
               taued(1:2) = W%v(1:2,ednext) - W%v(1:2,ed)
               vsizeed(ed) = norm2( taued(1:2) )
               perim = perim + vsizeed(ed)
            end do
            sizes(i,1:10)=vsizeed(1:10)
            avg(i) = perim / (W%n-1)

            call compJ2(nsites,sites,i,W,pdata%J2ctol,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug)
            if ( flag .ne. 0 ) return

            ji2(i) = ftmp
        end if

        call voronoi_cell_destroy(V)
        
        call polygon_destroy(W)

      end do
      
      deallocate(P%v,P%e,stat=allocerr)
      if ( allocerr .ne. 0 ) then
         write(*,*) 'In evalfg: Deallocation error.'
         flag = - 1
         return
      end if
  
    end subroutine evalsizeedJi2
  
   ! *******************************************************************
   ! *******************************************************************


   !! This subroutine give us the minimun value of |E| and minimum value of \theta (angle between two edges) in the whole diagram
   !! Also gives a vector with the average size of the cell's of the edges of each Voronoi cell.

    subroutine eval(n,x,e,a,avgsizeed,flag,pdataptr)
      
      implicit none
    
      ! SCALAR ARGUMENTS
      integer, intent(in) :: n
      integer, intent(out) :: flag
      real(kind=8), intent(out) :: e,a
      type(c_ptr), optional, intent(in) :: pdataptr
        
      ! ARRAY ARGUMENTS
      real(kind=8), intent(in) :: x(n)
      real(kind=8), intent(out) :: avgsizeed(n/2)

      ! PARAMETERS
      logical, parameter ::  debug = .false.
      
      ! LOCAL SCALARS
      integer :: allocerr,i,ierror,k,maxnv,nsites,ntriag,status,ed,ednext,edprev
      real(kind=8) :: perim
      type(pdata_type), pointer :: pdata
      type(VoronoiCell) :: V
      type(Polygon) :: P,W
      
      ! LOCAL ARRAYS
      integer :: indices(n/2),sitetv(3*2*(n/2)),start((n/2)+1),til(3,2*(n/2)),tnbr(3,2*(n/2))
      real(kind=8) :: fden(n/2),gden(2,n/2),sites(2,n/2),vorxy(2,2*(n/2)),valf(n/2),taued(2),tauedprev(2), &
                     vminang(n/2),vminsizeed(n/2),vsizeed(10),vang(10)
          
      flag = 0
      
      call c_f_pointer(pdataptr,pdata)
  
      nsites = n / 2
      indices(1:nsites) = (/ (i,i=1,nsites) /)
      sites(1:2,1:nsites) = reshape( x(1:n), (/ 2, nsites /) )
          
      if ( nsites .lt. 3 ) then
         write(*,*) 'In evalfg: nsites must be at least 3!'
         flag = - 1
         return
      end if
  
      ! Compute Delaunay triangulation
  
      call dtris2(nsites,sites,indices,ntriag,til,tnbr,ierror)
      
      if ( ierror .ne. 0 .and. ierror .ne. 225 ) then
         write(*,*) 'In evalfg: dtris2 returned an error different from 0 and 225: ',ierror
         flag = - 1
         return
      end if
  
      call voronoi_cell_patch(nsites,ntriag,til,start,sitetv)
      
      ! Compute vertices of the Voronoi diagram
            
      if ( ntriag .gt. 2 * nsites ) then
         write(*,*) 'In evalfg: 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!'
         flag = - 1
         return
      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 evalfg: triangle_circumcenter_2d returned NaN.'
            flag = - 1
            return
         end if
      end do

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

      allocate(P%v(2,maxnv),P%e(maxnv-1),stat=allocerr)

      if ( allocerr .ne. 0 ) then
         write(*,*) 'In evalfg: Allocation error.'
         flag = - 1
         return
      end if

      P%n = pdata%nvpols(1)
      P%deg = .false.
         
      P%v(1,1:P%n) = pdata%xpol(1,1:P%n)
      P%v(2,1:P%n) = pdata%ypol(1,1:P%n)
      P%e(1:P%n-1) = pdata%poledges(1,1:P%n-1)
  
      do i = 1,nsites
        if ( debug ) then
            write(*,*)
            write(*,*) 'SITE i = ',i
            write(*,*) 'Generator: ',sites(1:2,i)
            write(*,*)
        end if
  
        call voronoi_cell(i,nsites,sites,start,sitetv,ntriag,til,tnbr,vorxy,V,status)
        if ( status .ne. 0 ) then
            write(*,*) 'In evalfg: Error in voronoi_cell, status = ',status
            flag = - 1
            return
        end if


        if ( nsites .ge. 3 .and. ierror .eq. 0 ) then
            call VorCellInterConvPol(sites(1:2,i),P,V,W)
        else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then 
            call VorCellInterConvPol_Collinear(i,nsites,sites,P,W)
        end if
        
            
        if ( W%deg ) then
            if ( debug ) then
               write(*,*)
               write(*,*) 'POLYGON Wi with i = ',i, ' IS EMPTY.'
            end if
            write(*,*) 'Degenerate polygon, do not anything.'

        else
            if ( debug ) then
               write(*,*)
               write(*,*) 'POLYGON Wi with i = ',i
               write(*,*) 'Number of vertices: ',W%n - 1
               write(*,*)
            end if
            perim = 0.0d0
            do ed=1,W%n-1
               if ( ed .eq. W%n-1 ) then
                  edprev = ed - 1
                  ednext = 1
               else if ( ed .eq. 1) then
                  edprev = W%n-1
                  ednext = ed + 1
               else
                  edprev = ed - 1
                  ednext = ed + 1
               end if
               taued(1:2) = W%v(1:2,ednext) - W%v(1:2,ed)
               tauedprev(1:2)= W%v(1:2,ed) - W%v(1:2,edprev)
               vsizeed(ed) = norm2( taued(1:2) ) 
               vang(ed) = acos( max( -1.0d0, min( dot_product( -tauedprev(1:2)/norm2( tauedprev(1:2) ), taued(1:2)/vsizeed(ed) ), 1.0d0) ) )
               perim = perim + vsizeed(ed)
            end do
            vminang(i) = minval( vang(1:W%n-1) )
            vminsizeed(i) = minval( vsizeed(1:W%n-1) )
            avgsizeed(i) = perim/(W%n-1)

            if ( flag .ne. 0 ) return

         end if

        call voronoi_cell_destroy(V)
        
        call polygon_destroy(W)

      end do
      e = minval(vminsizeed(1:nsites))
      a = minval(vminang(1:nsites))
      
      deallocate(P%v,P%e,stat=allocerr)
      if ( allocerr .ne. 0 ) then
         write(*,*) 'In evalfg: Deallocation error.'
         flag = - 1
         return
      end if
  
    end subroutine eval

   ! ******************************************************************
   ! ******************************************************************
   
   subroutine compG(nsites,sites,i,W,f,gnnzmax,gnnz,gind,gval,flag,compgrad,debug)

      implicit none

      ! SCALAR ARGUMENTS
      logical, intent(in) :: compgrad,debug
      integer, intent(in) :: gnnzmax,i,nsites
      integer, intent(out) :: gnnz
      integer, intent(inout) :: flag
      real(kind=8), intent(out) :: f
      type(polygon), intent(in) :: W

      ! ARRAY ARGUMENTS
      integer, intent(out) :: gind(gnnzmax)
      real(kind=8), intent(in) :: sites(2,nsites)
      real(kind=8), intent(out) :: gval(gnnzmax)

      ! LOCAL SCALARS
      integer :: ed,ednext,nx,ny,j,k,r
      real(kind=8) :: vdet,volW,tmp,tinteg,tintegex

      ! LOCAL ARRAYS
      real(kind=8) :: diff1(2),diff2(2),cent(2),u(2),diff(nsites),v1(2),v2(2),v3(2),v4(3)
      integer :: minindex(1)
      ! LOCAL PARAMETERS
      real(kind=8), parameter :: h = 1.0d-4

      if ( compgrad ) then        
         
         call centroid(W,cent,volW)

         if ( debug ) then
            write(*,*)
            write(*,*) 'volWi= ',volW
            write(*,*) 'centi= ',cent(1:2)
         end if

         gnnz = 0

         if ( gnnz + 2 .gt. gnnzmax ) then
               write(*,*) 'In compG: Increase gnnzmax to at least ',gnnz+2
               flag = - 1
               return
         end if
         
         gind(gnnz+1:gnnz+2) = (/ 2*i-1, 2*i /)
         gval(gnnz+1:gnnz+2) = 0.0
         gnnz = gnnz + 2

         gval(1:2) = gval(1:2) + 2.0d0 * volW * ( sites(1:2,i) - cent(1:2)  )

      end if
      f = 0.0d0
      do ed = 1,W%n-1

         if ( ed .eq. W%n-1 ) then
            ednext = 1
         else 
            ednext = ed + 1
         end if

         if ( debug ) then
               write(*,*)
               write(*,*) '===================='
               write(*,*) 'Processing edge ed = ',ed
               write(*,*) 'Its inner vertex, named ',ed,' as well, corresponds to the point ',W%v(1:2,ed)
               write(*,*) 'Its outer vertex, named ',ednext,'as well, corresponds to the point ',W%v(1:2,ednext)
               write(*,*) '===================='
               write(*,*)
         end if



         if ( debug ) then
               write(*,*) 'Next: ',ednext
               write(*,*) 'It corresponds to the point',W%v(1:2,ednext)
         end if

         diff1(1:2) = W%v(1:2,ed) - sites(1:2,i)
         diff2(1:2) = W%v(1:2,ednext) - sites(1:2,i)
         vdet = ( diff1(1) * diff2(2) ) - ( diff1(2) * diff2(1) )  
         ! vdet = ( W%v(1,ed) - sites(1,i) ) * ( W%v(2,ednext) - sites(2,i) ) - &
         !        ( W%v(2,ed) - sites(2,i) ) * ( W%v(1,ednext) - sites(1,i) )
         !call det(W%v(1:2,ed), sites(1:2,i), W%v(1:2,ednext),vdet)

         if ( debug ) then
               write(*,*)
               write(*,*)'v_E - a_i = ',diff1(1:2)
               write(*,*)'w_E - a_i = ',diff2(1:2)
               write(*,*)'det = ',vdet
               write(*,*)'||v_E-a_i||^2= ', sum( diff1(1:2)**2 )
               write(*,*)'(v_E-a_i)^T(w_E-a_i)', dot_product( diff1(1:2),diff2(1:2) )
               write(*,*)'||w_E-a_i||^2= ',sum( diff2(1:2)**2 ) 
               write(*,*)
         end if 

         tmp = ( abs( vdet ) ) * ( sum( diff1(1:2)**2 ) + &
         dot_product( diff1(1:2),diff2(1:2) ) + sum( diff2(1:2)**2 ) )

         f = f  + tmp

         if ( debug ) then
            write(*,*)
            write(*,*)'tmp = ', tmp
            write(*,*)'f = ', f
            write(*,*)
         end if
      
      end do
   
      f = f /12.0d0
      if ( debug ) then
         write(*,*)'f= ', f
      end if

   end subroutine compG

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

   subroutine normalvec(v,normv)

      implicit none

      ! ARRAY ARGUMENTS
      real(kind=8), intent(in) :: v(2)
      real(kind=8), intent(out) :: normv(2)

      normv(1:2) = (/ -v(2), v(1) /)

   end subroutine normalvec

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

   subroutine compGalternative(nsites,sites,i,W,f,gnnzmax,gnnz,gind,gval,flag,compgrad,debug)
    
      implicit none

      ! SCALAR ARGUMENTS
      logical, intent(in) :: compgrad,debug
      integer, intent(in) :: gnnzmax,i,nsites
      integer, intent(out) :: gnnz
      integer, intent(inout) :: flag
      real(kind=8), intent(out) :: f
      type(polygon), intent(in) :: W

      ! ARRAY ARGUMENTS
      integer, intent(out) :: gind(gnnzmax)
      real(kind=8), intent(in) :: sites(2,nsites)
      real(kind=8), intent(out) :: gval(gnnzmax)

      ! LOCAL SCALARS
      integer :: ed, edprev, ednext
      real(kind=8) :: phia,phib,phic,xi,vdet,tmp,kap,absvdet,kappe,sig

      ! LOCAL ARRAYS
      real(kind=8) :: gradphi(2),M(2,2),wrot(2),vrot(2),p(2),diff1(2),diff2(2),ai(2),q(2)
      

      xi = 1.0d0/12.0d0
      f = 0.0d0

      if ( debug ) then
         write(*,*)
         write(*,*) 'xi = ', xi
         write(*,*)
      end if

      if ( compgrad ) then
         gnnz = 0
      
         if ( gnnz + 2 .gt. gnnzmax ) then
            write(*,*) 'In compGalternative: Increase gnnzmax to at least ',gnnz+2
            flag = - 1
            return
         end if
         
         gind(gnnz+1:gnnz+2) = (/ 2*i-1, 2*i /)
         gval(gnnz+1:gnnz+2) = 0.0d0
         gnnz = gnnz + 2
      end if
         
      do ed = 1,W%n-1

         if ( debug ) then
            write(*,*)
            write(*,*) '===================='
            write(*,*) 'Processing edge ed = ',ed
            write(*,*) 'Its inner vertex, named ',ed,' as well, corresponds to the point ',W%v(1:2,ed)
            write(*,*) '===================='
            write(*,*)
         end if

         if ( ed .eq. 1 ) then
            edprev = W%n - 1
            ednext = ed + 1
         else if ( ed .eq. W%n-1 ) then
            edprev = ed - 1
            ednext = 1
         else 
            edprev = ed - 1
            ednext = ed + 1
         end if

         ! Calculating the objective function

         diff1(1:2) = W%v(1:2,ed) - sites(1:2,i)
         diff2(1:2) = W%v(1:2,ednext) - sites(1:2,i)
         vdet = ( diff1(1) * diff2(2) ) - ( diff1(2) * diff2(1) )  

         if ( debug ) then
               write(*,*)
               write(*,*)'v_E - a_i = ',diff1(1:2)
               write(*,*)'w_E - a_i = ',diff2(1:2)
               write(*,*)'det = ',vdet
               write(*,*)'||v_E-a_i||^2= ', sum( diff1(1:2)**2 )
               write(*,*)'(v_E-a_i)^T(w_E-a_i)', dot_product( diff1(1:2),diff2(1:2) )
               write(*,*)'||w_E-a_i||^2= ',sum( diff2(1:2)**2 ) 
               write(*,*)
         end if 

         kap = sum( diff1(1:2)**2 ) + dot_product( diff1(1:2),diff2(1:2) ) + sum( diff2(1:2)**2 )
         absvdet = abs(vdet) 
         tmp = absvdet * kap

         f = f + tmp
         ! f = f  + absvdet

         if ( debug ) then
            write(*,*)
            write(*,*)'tmp = ', tmp
            write(*,*)'f = ', f
            write(*,*)
         end if
         
         if ( compgrad ) then

            ! Processing edge (v_E,w_E)

            if ( debug ) then
               write(*,*) 'Previous: ', edprev, 'next: ',ednext
               write(*,*) 'They corresponds to the points ',W%v(1:2,edprev),' and ', W%v(1:2,ednext),' respectively.'
            end if

            vrot(1:2) = (/ -W%v(2,ed), W%v(1,ed) /)
            wrot(1:2) = (/ W%v(2,ednext), -W%v(1,ednext) /)
            ai(1:2) = (/ -sites(2,i), sites(1,i)/)
            kappe = sign(1.0d0,vdet) * kap
            p(1:2) = kappe*(wrot(1:2)+ai(1:2))+absvdet*(2.0d0*W%v(1:2,ed) + W%v(1:2,ednext)-3.0d0*sites(1:2,i))
            q(1:2) = kappe*(vrot(1:2)-ai(1:2))+absvdet*(W%v(1:2,ed)+2.0d0*W%v(1:2,ednext)-3.0d0*sites(1:2,i))

            sig = sign(1.0d0,vdet)


            if ( debug ) then
               write(*,*)
               write(*,*) 'vrot= ',vrot
               write(*,*) 'wrot= ',wrot
               write(*,*) 'ai*= ', ai
               write(*,*) 'k(E)= ',kappe
               write(*,*)
            end if
            
            gval(1:2) = gval(1:2) - xi*(3.0d0*absvdet * (W%v(1:2,ed) + W%v(1:2,ednext)-2.0d0*sites(1:2,i)) + &
                                    kappe*(wrot(1:2)+vrot(1:2)))
            ! if ( i .eq. 1 ) then
            !    write(*,*)'first term just depends of \delta ai.'
            !    write(*,*)'summed(1:2)= ', - xi*(absvdet * (9.0d0*W%v(1:2,ed) + 5.0d0*W%v(1:2,ednext)-14.0d0*sites(1:2,i)) + &
            !    kappe*(wrot(1:2)+vrot(1:2)))
            !    write(*,*)'gi(1:2)= ',gval(1:2)
            !    write(*,*)
            ! end if

            ! gval(1:2) = gval(1:2) - sig * xi * ( wrot(1:2) + vrot(1:2))
         
            ! Starting point v_E

            if ( debug ) then
               write(*,*)
               write(*,*) 'Processing v_E'
               write(*,*)
            end if

            if ( W%e(ed) .lt. 0 .and. W%e(edprev) .lt. 0) then
               if ( debug ) then
                  write(*,*) 'We are in the case |P(v_E)|=3. ', &
                  'Neighbour cells are: ',abs(W%e(edprev)),' and ',abs(W%e(ed))
               end if

               call matrixint(sites(1:2,abs(W%e(edprev))),sites(1:2,abs(W%e(ed))),sites(1:2,i),W%v(1:2,ed),M,flag)
               if ( flag .ne. 0 ) return

               gval(1:2) = gval(1:2) + xi * matmul( transpose( M(1:2,1:2) ), p(1:2) )
               ! gval(1:2) = gval(1:2) + sig * xi * matmul( transpose( M(1:2,1:2) ), wrot(1:2) + ai(1:2) )
               ! if ( i .eq. 1 ) then
               !    write(*,*)'Case |P(v_E)|=3.'
               !    write(*,*)'summed(1:2)= ', xi * matmul( transpose( M(1:2,1:2) ), p(1:2) ) 
               !    write(*,*)'gi(1:2)= ',gval(1:2)
               !    write(*,*)
               ! end if

               call matrixint(sites(1:2,i),sites(1:2,abs(W%e(edprev))),sites(1:2,abs(W%e(ed))),W%v(1:2,ed),M,flag)
               if ( flag .ne. 0 ) return
            
               if ( gnnz + 2 .gt. gnnzmax) then
                  write(*,*) 'In compGalternative: Increase gnnzmax to at least ',gnnz + 2
                  flag = -1
                  return
               end if

               gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ed))-1, 2*abs(W%e(ed)) /)
               gval(gnnz+1:gnnz+2) = xi * matmul( transpose( M(1:2,1:2) ), p(1:2) )
               !gval(gnnz+1:gnnz+2) = sig * xi * matmul( transpose( M(1:2,1:2) ), wrot(1:2) + ai(1:2) )
               gnnz = gnnz + 2

               call matrixint(sites(1:2,abs(W%e(ed))),sites(1:2,i),sites(1:2,abs(W%e(edprev))),W%v(1:2,ed),M,flag)
               if ( flag .ne. 0 ) return

               if ( gnnz + 2 .gt. gnnzmax) then
                  write(*,*) 'In compGalternative: Increase gnnzmax to at least ',gnnz + 2
                  flag = -1
                  return
               end if

               gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(edprev))-1, 2*abs(W%e(edprev)) /)
               gval(gnnz+1:gnnz+2) = xi * matmul( transpose( M(1:2,1:2) ), p(1:2) )
               !gval(gnnz+1:gnnz+2) = sig * xi * matmul( transpose( M(1:2,1:2) ), wrot(1:2) + ai(1:2) )
               gnnz = gnnz + 2

            else if ( W%e(ed) .lt. 0 ) then
               if ( debug ) then
                  write(*,*) 'We are in the case |P(v_E)|=2 (type 1). Neighbour cell is: ',abs(W%e(ed))
               end if

               call line_exp2imp_2d(W%v(1:2,ed),W%v(1:2,edprev),phia,phib,phic)

               gradphi(1:2) = (/ phia, phib /)

               if ( gradphi(1) * sites(1,i) + gradphi(2) * sites(2,i) + phic .gt. 0.0d0 ) then
                  gradphi(1:2) = -gradphi(1:2)
               end if

               call matrixbd(sites(1:2,abs(W%e(ed))),sites(1:2,i),W%v(1:2,ed),gradphi,M,flag)
               if ( flag .ne. 0 ) return

               gval(1:2) = gval(1:2) + xi * matmul( transpose( M(1:2,1:2) ), p(1:2) )
               !gval(1:2) = gval(1:2) + sig * xi * matmul( transpose( M(1:2,1:2) ), wrot(1:2) + ai(1:2) )
               ! if ( i .eq. 1 ) then
               !    write(*,*)'Case |P(v_E)|=2 (type 1).'
               !    write(*,*)'summed(1:2)= ', xi * matmul( transpose( M(1:2,1:2) ), p(1:2) ) 
               !    write(*,*)'gi(1:2)= ',gval(1:2)
               !    write(*,*)
               ! end if

               call matrixbd(sites(1:2,i),sites(1:2,abs(W%e(ed))),W%v(1:2,ed),gradphi,M,flag)
               if ( flag .ne. 0 ) return

               if ( gnnz + 2 .gt. gnnzmax ) then
                  write(*,*) 'In compGalternative: Increase gnnzmax to at least ', gnnz+2
                  flag = -1
                  return
               end if

               gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ed))-1, 2*abs(W%e(ed)) /)
               gval(gnnz+1:gnnz+2) = xi * matmul( transpose( M(1:2,1:2) ), p(1:2) )
               !gval(gnnz+1:gnnz+2) = sig * xi * matmul( transpose( M(1:2,1:2) ), wrot(1:2) + ai(1:2) )
               gnnz = gnnz + 2

            else if ( W%e(edprev) .lt. 0 ) then
               if ( debug ) then
                  write(*,*) 'We are in the case |P(v_E)|=2 (type 2). Neighbour cell is: ', abs(W%e(edprev))
               end if

               call line_exp2imp_2d(W%v(1:2,ed),W%v(1:2,ednext),phia,phib,phic)

               gradphi(1:2) = (/ phia, phib /)

               if ( gradphi(1) * sites(1,i) + gradphi(2) * sites(2,i) + phic .gt. 0.0d0 ) then
                  gradphi(1:2) = -gradphi(1:2)
               end if

               call matrixbd(sites(1:2,abs(W%e(edprev))),sites(1:2,i),W%v(1:2,ed),gradphi,M,flag)
               if ( flag .ne. 0 ) return

               gval(1:2) = gval(1:2) + xi * matmul( transpose( M(1:2,1:2) ), p(1:2) )

               !gval(1:2) = gval(1:2) + sig * xi * matmul( transpose( M(1:2,1:2) ), wrot(1:2) + ai(1:2) )

               ! if ( i .eq. 1 ) then
               !    write(*,*)'Case |P(v_E)|=2 (type 2).'
               !    write(*,*)'summed(1:2)= ', xi * matmul( transpose( M(1:2,1:2) ), p(1:2) ) 
               !    write(*,*)'gi(1:2)= ',gval(1:2)
               !    write(*,*)
               ! end if

               call matrixbd(sites(1:2,i),sites(1:2,abs(W%e(edprev))),W%v(1:2,ed),gradphi,M,flag)
               if ( flag .ne. 0 ) return

               if ( gnnz + 2 .gt. gnnzmax ) then
                  write(*,*) 'In compGalternative: Increase gnnzmax to at least ', gnnz+2
                  flag = -1
                  return
               end if

               gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(edprev))-1, 2*abs(W%e(edprev)) /)
               gval(gnnz+1:gnnz+2) = xi * matmul( transpose( M(1:2,1:2) ), p(1:2) )
               !gval(gnnz+1:gnnz+2) = sig * xi * matmul( transpose( M(1:2,1:2) ), wrot(1:2) + ai(1:2) )
               gnnz = gnnz + 2
            end if

            ! Ending point w_E

            if ( debug ) then
               write(*,*)
               write(*,*) 'Processing w_E'
               write(*,*)
            end if

            if (W%e(ed) .lt. 0 .and. W%e(ednext) .lt. 0) then
               if ( debug ) then
                  write(*,*) 'We are in the case |P(w_E)|=3. ',&
                  'Neighbour cells are: ',abs(W%e(ed)),' and ',abs(W%e(ednext))
               end if

               call matrixint(sites(1:2,abs(W%e(ed))),sites(1:2,abs(W%e(ednext))),sites(1:2,i),W%v(1:2,ednext),M,flag)
               if ( flag .ne. 0 ) return

               gval(1:2) = gval(1:2) + xi * matmul( transpose( M(1:2,1:2) ), q(1:2) )
               !gval(1:2) = gval(1:2) + sig * xi * matmul( transpose( M(1:2,1:2) ), vrot(1:2) - ai(1:2) )

               ! if ( i .eq. 1 ) then
               !    write(*,*)'Case |P(w_E)|=3.'
               !    write(*,*)'summed(1:2)= ', xi * matmul( transpose( M(1:2,1:2) ), p(1:2) ) 
               !    write(*,*)'gi(1:2)= ',gval(1:2)
               !    write(*,*)
               ! end if

               call matrixint(sites(1:2,i),sites(1:2,abs(W%e(ed))),sites(1:2,abs(W%e(ednext))),W%v(1:2,ednext),M,flag)
               if ( flag .ne. 0) return

               if ( gnnz + 2 .gt. gnnzmax ) then
                  write(*,*) 'In compGalternative: Increase gnnzmax to at least ', gnnz+2
                  flag = -1
                  return
               end if

               gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ednext))-1, 2*abs(W%e(ednext)) /)
               gval(gnnz+1:gnnz+2) = xi * matmul( transpose( M(1:2,1:2) ), q(1:2) )
               !gval(gnnz+1:gnnz+2) = sig * xi * matmul( transpose( M(1:2,1:2) ), vrot(1:2) - ai(1:2) )
               gnnz = gnnz + 2

               call matrixint(sites(1:2,abs(W%e(ednext))),sites(1:2,i),sites(1:2,abs(W%e(ed))),W%v(1:2,ednext),M,flag)
               if ( flag .ne. 0 ) return

               if ( gnnz + 2 .gt. gnnzmax ) then
                  write(*,*) 'In compGalternative: Increase gnnzmax to at least ', gnnz+2
                  flag = -1
                  return
               end if

               gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ed))-1, 2*abs(W%e(ed)) /)
               gval(gnnz+1:gnnz+2) = xi * matmul( transpose( M(1:2,1:2) ), q(1:2) )
               !gval(gnnz+1:gnnz+2) = sig * xi * matmul( transpose( M(1:2,1:2) ), vrot(1:2) - ai(1:2) )
               gnnz = gnnz + 2

            else if ( W%e(ed) .lt. 0 ) then
               if ( debug ) then
                  write(*,*) 'We are in the case |P(w_E)|=2 (type 1). Neighbour cell is: ',abs(W%e(ed))
               end if

               call line_exp2imp_2d(W%v(1:2,ednext),W%v(1:2,ednext+1),phia,phib,phic)

               gradphi(1:2) = (/ phia, phib /)

               if ( gradphi(1) * sites(1,i) + gradphi(2) * sites(2,i) + phic .gt. 0.0d0 ) then
                  gradphi(1:2) = -gradphi(1:2)
               end if
               
               call matrixbd(sites(1:2,abs(W%e(ed))),sites(1:2,i),W%v(1:2,ednext),gradphi,M,flag)
               if ( flag .ne. 0 ) return

               gval(1:2) = gval(1:2) + xi * matmul( transpose( M(1:2,1:2) ), q(1:2) )
               !gval(1:2) = gval(1:2) + sig * xi * matmul( transpose( M(1:2,1:2) ), vrot(1:2) - ai(1:2) )

               ! if ( i .eq. 1 ) then
               !    write(*,*)'Case |P(w_E)|=2 (type 1).'
               !    write(*,*)'summed(1:2)= ', xi * matmul( transpose( M(1:2,1:2) ), p(1:2) ) 
               !    write(*,*)'gi(1:2)= ',gval(1:2)
               !    write(*,*)
               ! end if

               call matrixbd(sites(1:2,i),sites(1:2,abs(W%e(ed))),W%v(1:2,ednext),gradphi,M,flag)
               if ( flag .ne. 0 ) return

               if ( gnnz + 2 .gt. gnnzmax ) then
                  write(*,*) 'In compGalternative: Increase gnnzmax to at least ', gnnz+2
                  flag = -1
                  return
               end if

               gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ed))-1, 2*abs(W%e(ed)) /)
               gval(gnnz+1:gnnz+2) = xi * matmul( transpose( M(1:2,1:2) ), q(1:2) )
               !gval(gnnz+1:gnnz+2) = sig * xi * matmul( transpose( M(1:2,1:2) ), vrot(1:2) - ai(1:2) )
               gnnz = gnnz + 2

            else if ( W%e(ednext) .lt. 0 ) then
               if ( debug ) then
                  write(*,*) 'We are in the case |P(w_E)|=2 (type 2). Neighbour cell is: ',abs(W%e(ednext))
               end if

               call line_exp2imp_2d(W%v(1:2,ednext),W%v(1:2,ed),phia,phib,phic)

               gradphi(1:2) = (/ phia, phib /)

               if ( gradphi(1) * sites(1,i) + gradphi(2) * sites(2,i) + phic .gt. 0.0d0 ) then
                  gradphi(1:2) = -gradphi(1:2)
               end if
               
               call matrixbd(sites(1:2,abs(W%e(ednext))),sites(1:2,i),W%v(1:2,ednext),gradphi,M,flag)
               if ( flag .ne. 0 ) return

               gval(1:2) = gval(1:2) + xi * matmul( transpose( M(1:2,1:2) ), q(1:2) )
               !gval(1:2) = gval(1:2) + sig * xi * matmul( transpose( M(1:2,1:2) ), vrot(1:2) - ai(1:2) )

               ! if ( i .eq. 1 ) then
               !    write(*,*)'Case |P(w_E)|=2 (type 2).'
               !    write(*,*)'summed(1:2)= ', xi * matmul( transpose( M(1:2,1:2) ), p(1:2) ) 
               !    write(*,*)'gi(1:2)= ',gval(1:2)
               !    write(*,*)
               ! end if

               call matrixbd(sites(1:2,i),sites(1:2,abs(W%e(ednext))),W%v(1:2,ednext),gradphi,M,flag)
               if ( flag .ne. 0 ) return

               if ( gnnz + 2 .gt. gnnzmax ) then
                  write(*,*) 'In compGalternative: Increase gnnzmax to at least ', gnnz+2
                  flag = -1
                  return
               end if

               gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ednext))-1, 2*abs(W%e(ednext)) /)
               gval(gnnz+1:gnnz+2) = xi * matmul( transpose( M(1:2,1:2) ), q(1:2) )
               !gval(gnnz+1:gnnz+2) = sig * xi * matmul( transpose( M(1:2,1:2) ), vrot(1:2) - ai(1:2) )
               gnnz = gnnz + 2
            end if
         end if

      end do

      f = xi * f

      if ( debug ) then
         write(*,*) 'f = ',f
      end if
   end subroutine compGalternative

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

   subroutine compallJ(n,x,f,debug,flag,pdataptr)
      
      implicit none
    
      ! SCALAR ARGUMENTS
      logical, intent(in) :: debug
      integer, intent(in) :: n
      integer, intent(out) :: flag
      real(kind=8), intent(out) :: f(0:3,n/2)
      type(c_ptr), optional, intent(in) :: pdataptr
        
      ! ARRAY ARGUMENTS
      real(kind=8), intent(in) :: x(n)
  
      ! PARAMETERS
      logical, parameter :: compgrad = .false.
      integer, parameter :: gtmpnnzmax = 1
      
      ! LOCAL SCALARS
      integer :: allocerr,gtmpnnz,i,ierror,k,maxnv,nsites,ntriag,status
      real(kind=8) :: ftmp
      type(pdata_type), pointer :: pdata
      type(VoronoiCell) :: V
      type(Polygon) :: P,W
      
      ! LOCAL ARRAYS
      integer :: gtmpind(gtmpnnzmax),indices(n/2),sitetv(3*2*(n/2)),start((n/2)+1),til(3,2*(n/2)),tnbr(3,2*(n/2))
      real(kind=8) :: fden(n/2),gden(2,n/2),gtmpval(gtmpnnzmax),sites(2,n/2),vorxy(2,2*(n/2))
   
          
      flag = 0
      
      call c_f_pointer(pdataptr,pdata)
  
      nsites = n / 2
      indices(1:nsites) = (/ (i,i=1,nsites) /)
      sites(1:2,1:nsites) = reshape( x(1:n), (/ 2, nsites /) )
          
      if ( nsites .lt. 3 ) then
         write(*,*) 'In evalfg: nsites must be at least 3!'
         flag = - 1
         return
      end if
  
      ! Compute Delaunay triangulation
  
      call dtris2(nsites,sites,indices,ntriag,til,tnbr,ierror)
      
      if ( ierror .ne. 0 .and. ierror .ne. 225 ) then
         write(*,*) 'In evalfg: dtris2 returned an error different from 0 and 225: ',ierror
         flag = - 1
         return
      end if
    
      call voronoi_cell_patch(nsites,ntriag,til,start,sitetv)
      
      ! Compute vertices of the Voronoi diagram
            
      if ( ntriag .gt. 2 * nsites ) then
         write(*,*) 'In evalfg: 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!'
         flag = - 1
         return
      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 evalfg: triangle_circumcenter_2d returned NaN.'
            flag = - 1
            return
         end if
      end do

      ! Compute density function at sites if required

      do i = 1, nsites
         call fgcompJ4(sites(1:2,i),fden(i),gden(1:2,i),compgrad,pdataptr)
      end do
      
      maxnv = maxval( pdata%nvpols(1:pdata%npols) )
        
      allocate(P%v(2,maxnv),P%e(maxnv-1),stat=allocerr)

      if ( allocerr .ne. 0 ) then
         write(*,*) 'In evalfg: Allocation error.'
         flag = - 1
         return
      end if

      P%n = pdata%nvpols(1)
      P%deg = .false.
         
      P%v(1,1:P%n) = pdata%xpol(1,1:P%n)
      P%v(2,1:P%n) = pdata%ypol(1,1:P%n)
      P%e(1:P%n-1) = pdata%poledges(1,1:P%n-1)
  
      do i = 1,nsites
         if ( debug ) then
            write(*,*)
            write(*,*) 'SITE i = ',i
            write(*,*) 'Generator: ',sites(1:2,i)
            write(*,*)
         end if
  
         call voronoi_cell(i,nsites,sites,start,sitetv,ntriag,til,tnbr,vorxy,V,status)
         if ( status .ne. 0 ) then
            write(*,*) 'In evalfg: Error in voronoi_cell, status = ',status
            flag = - 1
            return
         end if
  
         
         if ( nsites .ge. 3 .and. ierror .eq. 0 ) then
            call VorCellInterConvPol(sites(1:2,i),P,V,W)
         else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then
            call VorCellInterConvPol_Collinear(i,nsites,sites,P,W)
         end if

         call compJ4(nsites,sites,i,W,pdata%volA,fden,gden,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug)
         if ( flag .ne. 0 ) return
         
         f(3,i) = ftmp

         if ( W%deg ) then
            if ( debug ) then
               write(*,*)
               write(*,*) 'POLYGON Wi with i = ',i, ' IS EMPTY.'
            end if
        
            f(0,i) = 1.0d0
            f(1,i) = 1.0d0
            f(2,i) = 1.0d0
            ! f(3,i) = 1.0d0

        else
            if ( debug ) then
               write(*,*)
               write(*,*) 'POLYGON Wi with i = ',i
               write(*,*) 'Number of vertices: ',W%n - 1
               write(*,*)
            end if
  
            call compJ1(nsites,sites,i,W,pdata%volA,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug)
            if ( flag .ne. 0 ) return
  
            f(0,i) = ftmp
  
            call compJ2(nsites,sites,i,W,pdata%J2ctol,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug)
            if ( flag .ne. 0 ) return

            f(1,i) = ftmp
                        
            call compJ3(nsites,sites,i,W,pdata%J3ctol,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug)
            if ( flag .ne. 0 ) return
  
            f(2,i) = ftmp

         end if

         call voronoi_cell_destroy(V)
  
         call polygon_destroy(W)
      end do
      
      deallocate(P%v,P%e,stat=allocerr)
      if ( allocerr .ne. 0 ) then
         write(*,*) 'In evalfg: Deallocation error.'
         flag = - 1
         return
      end if
    end subroutine compallJ

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

   subroutine compJ1(nsites,sites,i,W,volA,f,gnnzmax,gnnz,gind,gval,flag,compgrad,debug)
      
      implicit none
    
      ! SCALAR ARGUMENTS
      logical, intent(in) :: compgrad,debug
      integer, intent(in) :: gnnzmax,i,nsites
      integer, intent(out) :: gnnz
      integer, intent(inout) :: flag
      real(kind=8), intent(in) :: volA
      real(kind=8), intent(out) :: f
      type(Polygon), intent(in) :: W
        
      ! ARRAY ARGUMENTS
      integer, intent(out) :: gind(gnnzmax)
      real(kind=8), intent(in) :: sites(2,nsites)
      real(kind=8), intent(out) :: gval(gnnzmax)
  
      ! LOCAL SCALARS
      integer :: ed,ednext,j
      real(kind=8) :: mu,vol,xi
          
      ! LOCAL ARRAYS
      real(kind=8) :: p(2)
  
            
      call polygon_area_2d(W%n-1,W%v(1:2,1:W%n-1),vol)

      xi = nsites / volA
  
      if ( debug ) then
         write(*,*)
         write(*,*) 'xi = ',xi
         write(*,*)
      end if
         
      f = xi * vol - 1.0d0
  
      if ( .not. compgrad ) then
         return
      end if
      
      if ( 2 .gt. gnnzmax ) then
         write(*,*) 'In compJ1: Increase gnnzmax to at least ',2
         flag = - 1
         return
      end if
      
      gind(1:2) = (/ 2*i-1, 2*i /)
      gval(1:2) = 0.0d0
      gnnz = 2
      
      if ( .not. W%deg ) then
         if ( debug ) then
            write(*,*)
            write(*,*) 'POLYGON Wij with i = ',i,' and j = ',j
            write(*,*) 'Number of vertices: ',W%n - 1
            write(*,*)
         end if
            
         do ed = 1,W%n - 1
            if ( debug ) then
               write(*,*)
               write(*,*) '===================='
               write(*,*) 'Processing edge ed = ',ed
               write(*,*) 'Its inner vertex, named ',ed,' as well, corresponds to the point ',W%v(1:2,ed)
               write(*,*) 'Its outer vertex, named ',ednext,' as well, corresponds to the point ',W%v(1:2,ednext)
               write(*,*) '===================='
               write(*,*)
            end if
            
            if ( W%e(ed) .lt. 0 ) then
               
               if ( ed .eq. W%n - 1 ) then
                  ednext = 1
               else
                  ednext = ed + 1
               end if
               
               if ( debug ) then
                  write(*,*) 'Next: ',ednext
                  write(*,*) 'It corresponds to the point ',W%v(1:2,ednext)
               end if
               
               mu = norm2( W%v(1:2,ed) - W%v(1:2,ednext) ) / norm2( sites(1:2,i) - sites(1:2,abs(W%e(ed))) )
                  
               if ( debug ) then
                  write(*,*) 'mu = ',mu
               end if

               p(1:2) = 0.5d0 * ( W%v(1:2,ed) + W%v(1:2,ednext) )
               
               gval(1:2) = gval(1:2) + xi * mu * ( p(1:2) - sites(1:2,i) )
                  
               if ( gnnz + 2 .gt. gnnzmax ) then
                  write(*,*) 'In compJ1: Increase gnnzmax to at least ',gnnz+2
                  flag = - 1
                  return
               end if
         
               gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ed))-1, 2*abs(W%e(ed)) /)
               gval(gnnz+1:gnnz+2) = - xi * mu * ( p(1:2) - sites(1:2,abs(W%e(ed))) )
               gnnz = gnnz + 2

            else
               if ( debug ) then
                  write(*,*) 'This is a fixed edge.'
                  write(*,*) 'So, its vertices v_E and w_E are ignored.'
               end if
            end if
         end do

      else
         if ( debug ) then
            write(*,*)
            write(*,*) 'POLYGON Wi with i = ',i,' IS EMPTY.'
            write(*,*)
         end if
      end if
      
   end subroutine compJ1
  
    ! ******************************************************************
    ! ******************************************************************

   subroutine compJ2(nsites,sites,i,W,ctol,f,gnnzmax,gnnz,gind,gval,flag,compgrad,debug)
      
      implicit none
    
      ! SCALAR ARGUMENTS
      logical, intent(in) :: compgrad,debug
      integer, intent(in) :: gnnzmax,i,nsites
      integer, intent(out) :: gnnz
      integer, intent(inout) :: flag
      real(kind=8), intent(in) :: ctol
      real(kind=8), intent(out) :: f
      type(Polygon), intent(in) :: W
        
      ! ARRAY ARGUMENTS
      integer, intent(out) :: gind(gnnzmax)
      real(kind=8), intent(in) :: sites(2,nsites)
      real(kind=8), intent(out) :: gval(gnnzmax)
  
      ! LOCAL SCALARS
      integer :: ed,edprev,ednext
      real(kind=8) :: glcnst,mu,perim,phia,phib,phic,sizeed,tmp,xi
          
      ! LOCAL ARRAYS
      real(kind=8) :: gradphi(2),M(2,2),tang(2)
  
      !write(*,*) 'Site i = ',i
      perim = 0.0d0
      do ed = 1,W%n - 1
         if ( ed .lt. W%n - 1 ) then
            ednext = ed + 1
         else
            ednext = 1
         end if
  
         sizeed = norm2( W%v(1:2,ednext) - W%v(1:2,ed) )
         perim = perim + sizeed
  
        !write(*,*) 'Edge ',ed,' size = ',sizeed 
      end do
     !write(*,*) 'Average size = ',perim / ( W%n - 1 )
  
      if ( debug ) then
         write(*,*)
         write(*,*) 'Perimeter = ',perim
         write(*,*)
      end if
         
      xi = 2.0d0 / perim
            
      if ( debug ) then
         write(*,*)
         write(*,*) 'xi = ',xi
         write(*,*)
      end if
  
      f = 0.0d0
  
      if ( compgrad ) then
         glcnst = 0.0d0
         do ed = 1,W%n - 1
            if ( ed .lt. W%n - 1 ) then
               ednext = ed + 1
            else
               ednext = 1
            end if
            
            sizeed = norm2( W%v(1:2,ednext) - W%v(1:2,ed) )
            glcnst = glcnst + min( 0.0d0, ( W%n - 1 ) * sizeed / perim - ctol ) * sizeed / perim
         end do
         
         if ( debug ) then
            write(*,*)
            write(*,*) 'Global constant = ',glcnst
            write(*,*)
         end if
         
         gnnz = 0
         
         if ( gnnz + 2 .gt. gnnzmax ) then
            write(*,*) 'In compJ3: Increase gnnzmax to at least ',gnnz+2
            flag = - 1
            return
         end if
         
         gind(gnnz+1:gnnz+2) = (/ 2*i-1, 2*i /)
         gval(gnnz+1:gnnz+2) = 0.0d0
         gnnz = gnnz + 2
      end if
      
      do ed = 1,W%n - 1
         if ( debug ) then
            write(*,*)
            write(*,*) '===================='
            write(*,*) 'Processing edge ed = ',ed
            write(*,*) 'Its inner vertex, named ',ed,' as well, corresponds to the point ',W%v(1:2,ed)
            write(*,*) '===================='
            write(*,*)
         end if
         
         ! Processing edge (v_E,w_E)
                     
         if ( ed .eq. 1 ) then
            edprev = W%n - 1
            ednext = ed + 1
         else if ( ed .eq. W%n - 1 ) then
            edprev = ed - 1
            ednext = 1
         else
            edprev = ed - 1
            ednext = ed + 1
         end if
                    
         if ( debug ) then
            write(*,*) 'Previous: ',edprev,' next: ',ednext
            write(*,*) 'They correspond to the points ',W%v(1:2,edprev),' and ',W%v(1:2,ednext),' respectively.'
         end if
                  
         sizeed = norm2( W%v(1:2,ednext) - W%v(1:2,ed) )
  
         if ( debug ) then
            write(*,*) 'Edge size = ',sizeed
         end if
  
         tmp = min( 0.0d0, ( W%n - 1 ) * sizeed / perim - ctol )
  
         f = f + tmp ** 2 / ( W%n - 1 )
  
         if ( .not. compgrad ) then
            cycle
         end if
         
         mu = glcnst - tmp
               
         if ( debug ) then
            write(*,*) 'mu associated wit the edge = ',mu
         end if
                  
         if ( mu .eq. 0.0d0 ) then
            cycle
         end if
         
         tang(1:2) = ( W%v(1:2,ednext) - W%v(1:2,ed) ) / sizeed
  
         if ( debug ) then
            write(*,*) 'Tangent vector: ',tang(1:2)
         end if
  
         ! Starting point v_E
  
         if ( debug ) then
            write(*,*)
            write(*,*) 'Processing v_E'
            write(*,*)
         end if
  
         if ( W%e(ed) .lt. 0 .and. W%e(edprev) .lt. 0 ) then
            if ( debug ) then
               write(*,*) 'We are in the case |P(v_E)|=3. ', &
                    'Neighbour cells are: ',abs(W%e(edprev)),' and ',abs(W%e(ed))
            end if
                  
            call matrixint(sites(1:2,abs(W%e(edprev))),sites(1:2,abs(W%e(ed))),sites(1:2,i),W%v(1:2,ed),M,flag)
            if ( flag .ne. 0 ) return
                        
            gval(1:2) = gval(1:2) + xi * mu * matmul( transpose( M(1:2,1:2) ), tang(1:2) )
                     
            call matrixint(sites(1:2,i),sites(1:2,abs(W%e(edprev))),sites(1:2,abs(W%e(ed))),W%v(1:2,ed),M,flag)
            if ( flag .ne. 0 ) return
  
            if ( gnnz + 2 .gt. gnnzmax ) then
               write(*,*) 'In compJ3: Increase gnnzmax to at least ',gnnz+2
               flag = - 1
               return
            end if
          
            gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ed))-1, 2*abs(W%e(ed)) /)
            gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), tang(1:2) )
            gnnz = gnnz + 2
               
            call matrixint(sites(1:2,abs(W%e(ed))),sites(1:2,i),sites(1:2,abs(W%e(edprev))),W%v(1:2,ed),M,flag)
            if ( flag .ne. 0 ) return
  
            if ( gnnz + 2 .gt. gnnzmax ) then
               write(*,*) 'In compJ3: Increase gnnzmax to at least ',gnnz+2
               flag = - 1
               return
            end if
          
            gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(edprev))-1, 2*abs(W%e(edprev)) /)
            gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), tang(1:2) )
            gnnz = gnnz + 2
               
         else if ( W%e(ed) .lt. 0 ) then
            if ( debug ) then
               write(*,*) 'We are in the case |P(v_E)|=2 (type 1). Neighbour cell is: ',abs(W%e(ed))
            end if
                  
            call line_exp2imp_2d(W%v(1:2,ed),W%v(1:2,edprev),phia,phib,phic)
  
            gradphi(1:2) = (/ phia, phib /)
            
            if ( gradphi(1) * sites(1,i) + gradphi(2) * sites(2,i) + phic .gt. 0.0d0 ) then
               gradphi(1:2) = - gradphi(1:2)
            end if
      
            call matrixbd(sites(1:2,abs(W%e(ed))),sites(1:2,i),W%v(1:2,ed),gradphi,M,flag)
            if ( flag .ne. 0 ) return
  
            gval(1:2) = gval(1:2) + xi * mu * matmul( transpose( M(1:2,1:2) ), tang(1:2) )
                     
            call matrixbd(sites(1:2,i),sites(1:2,abs(W%e(ed))),W%v(1:2,ed),gradphi,M,flag)
            if ( flag .ne. 0 ) return
               
            if ( gnnz + 2 .gt. gnnzmax ) then
               write(*,*) 'In compJ3: Increase gnnzmax to at least ',gnnz+2
               flag = - 1
               return
            end if
          
            gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ed))-1, 2*abs(W%e(ed)) /)
            gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), tang(1:2) )
            gnnz = gnnz + 2
               
         else if ( W%e(edprev) .lt. 0 ) then
            if ( debug ) then
               write(*,*) 'We are in the case |P(v_E)|=2 (type 2). Neighbour cell is: ',abs(W%e(edprev))
            end if
                  
            call line_exp2imp_2d(W%v(1:2,ed),W%v(1:2,ednext),phia,phib,phic)
  
            gradphi(1:2) = (/ phia, phib /)
                  
            if ( gradphi(1) * sites(1,i) + gradphi(2) * sites(2,i) + phic .gt. 0.0d0 ) then
               gradphi(1:2) = - gradphi(1:2)
            end if
      
            call matrixbd(sites(1:2,abs(W%e(edprev))),sites(1:2,i),W%v(1:2,ed),gradphi,M,flag)
            if ( flag .ne. 0 ) return
  
            gval(1:2) = gval(1:2) + xi * mu * matmul( transpose( M(1:2,1:2) ), tang(1:2) )
                     
            call matrixbd(sites(1:2,i),sites(1:2,abs(W%e(edprev))),W%v(1:2,ed),gradphi,M,flag)
            if ( flag .ne. 0 ) return
                        
            if ( gnnz + 2 .gt. gnnzmax ) then
               write(*,*) 'In compJ3: Increase gnnzmax to at least ',gnnz+2
               flag = - 1
               return
            end if
          
            gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(edprev))-1, 2*abs(W%e(edprev)) /)
            gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), tang(1:2) )
            gnnz = gnnz + 2
         end if
                     
         ! Ending point w_E
                    
         if ( debug ) then
            write(*,*)
            write(*,*) 'Processing w_E'
            write(*,*)
         end if
            
         if ( W%e(ed) .lt. 0 .and. W%e(ednext) .lt. 0 ) then
            if ( debug ) then
               write(*,*) 'We are in the case |P(w_E)|=3. ', &
                    'Neighbour cells are: ',abs(W%e(ed)),' and ',abs(W%e(ednext))
            end if
                  
            call matrixint(sites(1:2,abs(W%e(ed))),sites(1:2,abs(W%e(ednext))),sites(1:2,i),W%v(1:2,ednext),M,flag)
            if ( flag .ne. 0 ) return
               
            gval(1:2) = gval(1:2) - xi * mu * matmul( transpose( M(1:2,1:2) ), tang(1:2) )
  
            call matrixint(sites(1:2,i),sites(1:2,abs(W%e(ed))),sites(1:2,abs(W%e(ednext))),W%v(1:2,ednext),M,flag)
            if ( flag .ne. 0 ) return
                   
            if ( gnnz + 2 .gt. gnnzmax ) then
               write(*,*) 'In compJ3: Increase gnnzmax to at least ',gnnz+2
               flag = - 1
               return
            end if
          
            gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ednext))-1, 2*abs(W%e(ednext)) /)
            gval(gnnz+1:gnnz+2) = - xi * mu * matmul( transpose( M(1:2,1:2) ), tang(1:2) )
            gnnz = gnnz + 2
            
            call matrixint(sites(1:2,abs(W%e(ednext))),sites(1:2,i),sites(1:2,abs(W%e(ed))),W%v(1:2,ednext),M,flag)
            if ( flag .ne. 0 ) return
                   
            if ( gnnz + 2 .gt. gnnzmax ) then
               write(*,*) 'In compJ3: Increase gnnzmax to at least ',gnnz+2
               flag = - 1
               return
            end if
          
            gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ed))-1, 2*abs(W%e(ed)) /)
            gval(gnnz+1:gnnz+2) = - xi * mu * matmul( transpose( M(1:2,1:2) ), tang(1:2) )
            gnnz = gnnz + 2
            
         else if ( W%e(ed) .lt. 0 ) then 
            if ( debug ) then
               write(*,*) 'We are in the case |P(w_E)|=2 (type 1). Neighbour cell is: ',abs(W%e(ed))
            end if
                  
            call line_exp2imp_2d(W%v(1:2,ednext),W%v(1:2,ednext+1),phia,phib,phic)
  
            gradphi(1:2) = (/ phia, phib /)
            
            if ( gradphi(1) * sites(1,i) + gradphi(2) * sites(2,i) + phic .gt. 0.0d0 ) then
               gradphi(1:2) = - gradphi(1:2)
            end if
      
            call matrixbd(sites(1:2,abs(W%e(ed))),sites(1:2,i),W%v(1:2,ednext),gradphi,M,flag)
            if ( flag .ne. 0 ) return
  
            gval(1:2) = gval(1:2) - xi * mu * matmul( transpose( M(1:2,1:2) ), tang(1:2) )
                     
            call matrixbd(sites(1:2,i),sites(1:2,abs(W%e(ed))),W%v(1:2,ednext),gradphi,M,flag)
            if ( flag .ne. 0 ) return
                     
            if ( gnnz + 2 .gt. gnnzmax ) then
               write(*,*) 'In compJ3: Increase gnnzmax to at least ',gnnz+2
               flag = - 1
               return
            end if
          
            gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ed))-1, 2*abs(W%e(ed)) /)
            gval(gnnz+1:gnnz+2) = - xi * mu * matmul( transpose( M(1:2,1:2) ), tang(1:2) )
            gnnz = gnnz + 2
            
         else if ( W%e(ednext) .lt. 0 ) then 
            if ( debug ) then
               write(*,*) 'We are in the case |P(w_E)|=2 (type 2). Neighbour cell is: ',abs(W%e(ednext))
            end if
            
            call line_exp2imp_2d(W%v(1:2,ednext),W%v(1:2,ed),phia,phib,phic)
               
            gradphi(1:2) = (/ phia, phib /)
                  
            if ( gradphi(1) * sites(1,i) + gradphi(2) * sites(2,i) + phic .gt. 0.0d0 ) then
               gradphi(1:2) = - gradphi(1:2)
            end if
               
            call matrixbd(sites(1:2,abs(W%e(ednext))),sites(1:2,i),W%v(1:2,ednext),gradphi,M,flag)
            if ( flag .ne. 0 ) return
  
            gval(1:2) = gval(1:2) - xi * mu * matmul( transpose( M(1:2,1:2) ), tang(1:2) )
            
            call matrixbd(sites(1:2,i),sites(1:2,abs(W%e(ednext))),W%v(1:2,ednext),gradphi,M,flag)
            if ( flag .ne. 0 ) return
                     
            if ( gnnz + 2 .gt. gnnzmax ) then
               write(*,*) 'In compJ3: Increase gnnzmax to at least ',gnnz+2
               flag = - 1
               return
            end if
          
            gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ednext))-1, 2*abs(W%e(ednext)) /)
            gval(gnnz+1:gnnz+2) = - xi * mu * matmul( transpose( M(1:2,1:2) ), tang(1:2) )
            gnnz = gnnz + 2
         end if
      end do
         
   end subroutine compJ2

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

   subroutine compJ4(nsites,sites,i,W,volA,fden,gden,f,gnnzmax,gnnz,gind,gval,flag,compgrad,debug)

      implicit none
      
      ! SCALAR ARGUMENTS
      logical, intent(in) :: compgrad,debug
      integer, intent(in) :: gnnzmax,i,nsites
      integer, intent(out) :: gnnz
      integer, intent(inout) :: flag
      real(kind=8), intent(in) :: volA
      real(kind=8), intent(out) :: f
      type(Polygon), intent(in) :: W

      ! ARRAY ARGUMENTS
      integer, intent(out) :: gind(gnnzmax)
      real(kind=8), intent(in) :: fden(nsites),gden(2,nsites),sites(2,nsites)
      real(kind=8), intent(out) :: gval(gnnzmax)

      ! LOCAL SCALARS
      integer :: ed,ednext,j
      real(kind=8) :: mu,vol,xi
          
      ! LOCAL ARRAYS
      real(kind=8) :: p(2)

      call polygon_area_2d(W%n-1,W%v(1:2,1:W%n-1),vol)

      xi = nsites / volA
  
      if ( debug ) then
         write(*,*)
         write(*,*) 'xi = ',xi
         write(*,*)
      end if

      f = xi * vol - fden(i)
  
      if ( .not. compgrad ) then
         return
      end if
      
      if ( 2 .gt. gnnzmax ) then
         write(*,*) 'In compJ4: Increase gnnzmax to at least ',2
         flag = - 1
         return
      end if

      gind(1:2) = (/ 2*i-1, 2*i /)
      gval(1:2) = -gden(1:2,i)
      gnnz = 2

      if ( .not. W%deg ) then
         if ( debug ) then
            write(*,*)
            write(*,*) 'POLYGON Wij with i = ',i,' and j = ',j
            write(*,*) 'Number of vertices: ',W%n - 1
            write(*,*)
         end if
            
         do ed = 1,W%n - 1
            if ( debug ) then
               write(*,*)
               write(*,*) '===================='
               write(*,*) 'Processing edge ed = ',ed
               write(*,*) 'Its inner vertex, named ',ed,' as well, corresponds to the point ',W%v(1:2,ed)
               write(*,*) 'Its outer vertex, named ',ednext,' as well, corresponds to the point ',W%v(1:2,ednext)
               write(*,*) '===================='
               write(*,*)
            end if
            
            if ( W%e(ed) .lt. 0 ) then
               
               if ( ed .eq. W%n - 1 ) then
                  ednext = 1
               else
                  ednext = ed + 1
               end if
               
               if ( debug ) then
                  write(*,*) 'Next: ',ednext
                  write(*,*) 'It corresponds to the point ',W%v(1:2,ednext)
               end if
               
               mu = norm2( W%v(1:2,ed) - W%v(1:2,ednext) ) / norm2( sites(1:2,i) - sites(1:2,abs(W%e(ed))) )
                  
               if ( debug ) then
                  write(*,*) 'mu = ',mu
               end if

               p(1:2) = 0.5d0 * ( W%v(1:2,ed) + W%v(1:2,ednext) )
               
               gval(1:2) = gval(1:2) + xi * mu * ( p(1:2) - sites(1:2,i) )
                  
               if ( gnnz + 2 .gt. gnnzmax ) then
                  write(*,*) 'In compJ1: Increase gnnzmax to at least ',gnnz+2
                  flag = - 1
                  return
               end if
         
               gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ed))-1, 2*abs(W%e(ed)) /)
               gval(gnnz+1:gnnz+2) = - xi * mu * ( p(1:2) - sites(1:2,abs(W%e(ed))) )
               gnnz = gnnz + 2

            else
               if ( debug ) then
                  write(*,*) 'This is a fixed edge.'
                  write(*,*) 'So, its vertices v_E and w_E are ignored.'
               end if
            end if
         end do

      else
         if ( debug ) then
            write(*,*)
            write(*,*) 'POLYGON Wi with i = ',i,' IS EMPTY.'
            write(*,*)
         end if
      end if

   end subroutine compJ4

  

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

   subroutine compJ3(nsites,sites,i,W,ctol,f,gnnzmax,gnnz,gind,gval,flag,compgrad,debug)
      
      implicit none
    
      ! SCALAR ARGUMENTS
      logical, intent(in) :: compgrad,debug
      integer, intent(in) :: gnnzmax,i,nsites
      integer, intent(out) :: gnnz
      integer, intent(inout) :: flag
      real(kind=8), intent(in) :: ctol
      real(kind=8), intent(out) :: f
      type(Polygon), intent(in) :: W
        
      ! ARRAY ARGUMENTS
      integer, intent(out) :: gind(gnnzmax)
      real(kind=8), intent(in) :: sites(2,nsites)
      real(kind=8), intent(out) :: gval(gnnzmax)
  
      ! LOCAL SCALARS
      integer :: ed,edprev,edpp,ednext,nmov
      real(kind=8) :: ang,avgang,glcnst,mu,phia,phib,phic,sizeed,sizeedprev,tmp,xi
          
      ! LOCAL ARRAYS
      real(kind=8) :: gradphi(2),M(2,2),normaled(2),normaledprev(2),taued(2),tauedprev(2),zeta(2)
  
      !write(*,*) 'Site i = ',i
      nmov = 0
      avgang = 0.0d0
      do ed = 1,W%n - 1
         if ( ed .eq. 1) then
            edprev = W%n - 1
            ednext = ed + 1
         else if ( ed .eq. W%n - 1 ) then
            edprev = ed - 1
            ednext = 1
         else
            edprev = ed - 1
            ednext = ed + 1
         end if
  
         ! If none is negative, the angle is fixed
         if ( W%e(ed) .lt. 0 .or. W%e(edprev) .lt. 0 ) then
            tauedprev(1:2) = W%v(1:2,ed) - W%v(1:2,edprev)
            taued(1:2)     = W%v(1:2,ednext) - W%v(1:2,ed)
  
            sizeedprev = norm2( tauedprev(1:2) )
            sizeed     = norm2( taued(1:2) )
         
            ang = acos( max( -1.0d0, min( dot_product( - tauedprev(1:2) / sizeedprev, taued(1:2) / sizeed ), 1.0d0 ) ) )
            !write(*,*) 'vertex = ',ed,' angle = ',ang
  
            if ( debug ) then
               write(*,*) 'vertex = ',ed,' angle = ',ang
            end if
         
            nmov = nmov + 1
            avgang = avgang + ang
         end if
      end do
      avgang = avgang / nmov
      !write(*,*) 'Average angle = ',avgang
  
      if ( debug ) then
         write(*,*)
         write(*,*) 'The polygon has nmov = ',nmov,' non-fixed angles.'
         write(*,*) 'Average angle = ',avgang
         write(*,*)
      end if
  
      xi = ( 2.0d0 / ( avgang * nmov ) )
             
      if ( debug ) then
         write(*,*)
         write(*,*) 'xi = ',xi
         write(*,*)
      end if
  
      f = 0.0d0
  
      if ( compgrad ) then
         glcnst = 0.0d0
         do ed = 1,W%n - 1
            if ( ed .eq. 1) then
               edprev = W%n - 1
               ednext = ed + 1
            else if ( ed .eq. W%n - 1 ) then
               edprev = ed - 1
               ednext = 1
            else
               edprev = ed - 1
               ednext = ed + 1
            end if
  
            ! If none is negative, the angle is fixed
            if ( W%e(ed) .lt. 0 .or. W%e(edprev) .lt. 0 ) then
               tauedprev(1:2) = W%v(1:2,ed) - W%v(1:2,edprev)
               taued(1:2)     = W%v(1:2,ednext) - W%v(1:2,ed)
               
               sizeedprev = norm2( tauedprev(1:2) )
               sizeed     = norm2( taued(1:2) )
               
               ang = acos( max( -1.0d0, min( dot_product( - tauedprev(1:2) / sizeedprev, taued(1:2) / sizeed ), 1.0d0 ) ) )
               
               glcnst = glcnst + ( ang / avgang ) * min( 0.0d0, ang / avgang - ctol ) / nmov
            end if
         end do
      
         gnnz = 0
      
         if ( gnnz + 2 .gt. gnnzmax ) then
            write(*,*) 'In compJ4: Increase gnnzmax to at least ',gnnz+2
            flag = - 1
            return
         end if
      
         gind(gnnz+1:gnnz+2) = (/ 2*i-1, 2*i /)
         gval(gnnz+1:gnnz+2) = 0.0d0
         gnnz = gnnz + 2
      end if
      
      do ed = 1,W%n - 1
         if ( debug ) then
            write(*,*)
            write(*,*) '===================='
            write(*,*) 'Processing edge ed = ',ed
            write(*,*) 'Its inner vertex, named ',ed,' as well, corresponds to the point ',W%v(1:2,ed)
            write(*,*) '===================='
            write(*,*)
         end if
         
         ! Processing edge (v_E,w_E)
                     
         if ( ed .eq. 1) then
            edprev = W%n - 1
            ednext = ed + 1
         else if ( ed .eq. W%n - 1 ) then
            edprev = ed - 1
            ednext = 1
         else
            edprev = ed - 1
            ednext = ed + 1
         end if
                    
         if ( debug ) then
            write(*,*) 'Previous: ',edprev,' next: ',ednext
            write(*,*) 'They correspond to the points ',W%v(1:2,edprev),' and ',W%v(1:2,ednext),' respectively.'
         end if
  
         ! If none is negative, the angle is fixed
         if ( W%e(ed) .lt. 0 .or. W%e(edprev) .lt. 0 ) then
            
            tauedprev(1:2) = W%v(1:2,ed) - W%v(1:2,edprev)
            taued(1:2)     = W%v(1:2,ednext) - W%v(1:2,ed)
  
            sizeedprev = norm2( tauedprev(1:2) )
            sizeed     = norm2( taued(1:2) )
         
            tauedprev(1:2) = tauedprev(1:2) / sizeedprev
            taued(1:2)     = taued(1:2) / sizeed
  
            ang = acos( max( -1.0d0, min( dot_product( - tauedprev(1:2), taued(1:2) ), 1.0d0 ) ) )
  
            tmp = min( 0.0d0, ang / avgang - ctol )
  
            f = f + tmp ** 2 / nmov
  
            if ( .not. compgrad ) then
               cycle
            end if
            
            mu = glcnst - tmp
  
            if ( debug ) then
               write(*,*) 'mu = ',mu
            end if
  
            if ( mu .eq. 0.0d0 ) then
               cycle
            end if
         
            call line_exp_normal_2d(W%v(1:2,ednext),W%v(1:2,ed),normaled)
            call line_exp_normal_2d(W%v(1:2,ed),W%v(1:2,edprev),normaledprev)
  
            if ( debug ) then
               write(*,*) 'Normal to ed     = ',normaled(1:2)
               write(*,*) 'Normal to edprev = ',normaledprev(1:2)
            end if
  
            ! Starting point v_E
  
            if ( debug ) then
               write(*,*)
               write(*,*) 'Processing v_E = w_{\hat E}'
               write(*,*)
            end if
  
          ! zeta(1:2) = ( dot_product( normaled(1:2), tauedprev(1:2) ) / sizeed ) * normaled(1:2) - &
          !             ( dot_product( normaledprev(1:2), taued(1:2) ) / sizeedprev ) * normaledprev(1:2)
            zeta(1:2) = normaled(1:2) / sizeed + normaledprev(1:2) / sizeedprev
  
            if ( debug ) then
               write(*,*) 'zeta = ',zeta(1:2)
            end if
                  
            if ( W%e(ed) .lt. 0 .and. W%e(edprev) .lt. 0 ) then
               if ( debug ) then
                  write(*,*) 'We are in the case |P(v_E)|=3. ', &
                       'Neighbour cells are: ',abs(W%e(edprev)),' and ',abs(W%e(ed))
               end if
                  
               call matrixint(sites(1:2,abs(W%e(edprev))),sites(1:2,abs(W%e(ed))),sites(1:2,i),W%v(1:2,ed),M,flag)
               if ( flag .ne. 0 ) return
                        
               gval(1:2) = gval(1:2) + xi * mu * matmul( transpose( M(1:2,1:2) ), zeta(1:2) )
                     
               call matrixint(sites(1:2,i),sites(1:2,abs(W%e(edprev))),sites(1:2,abs(W%e(ed))),W%v(1:2,ed),M,flag)
               if ( flag .ne. 0 ) return
  
               if ( gnnz + 2 .gt. gnnzmax ) then
                  write(*,*) 'In compJ4: Increase gnnzmax to at least ',gnnz+2
                  flag = - 1
                  return
               end if
          
               gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ed))-1, 2*abs(W%e(ed)) /)
               gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), zeta(1:2) )
               gnnz = gnnz + 2
  
               call matrixint(sites(1:2,abs(W%e(ed))),sites(1:2,i),sites(1:2,abs(W%e(edprev))),W%v(1:2,ed),M,flag)
               if ( flag .ne. 0 ) return
  
               if ( gnnz + 2 .gt. gnnzmax ) then
                  write(*,*) 'In compJ4: Increase gnnzmax to at least ',gnnz+2
                  flag = - 1
                  return
               end if
          
               gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(edprev))-1, 2*abs(W%e(edprev)) /)
               gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), zeta(1:2) )
               gnnz = gnnz + 2
  
            else if ( W%e(ed) .lt. 0 ) then
               if ( debug ) then
                  write(*,*) 'We are in the case |P(v_E)|=2 (type 1). Neighbour cell is: ',abs(W%e(ed))
               end if
               
               call line_exp2imp_2d(W%v(1:2,ed),W%v(1:2,edprev),phia,phib,phic)
  
               gradphi(1:2) = (/ phia, phib /)
            
               if ( gradphi(1) * sites(1,i) + gradphi(2) * sites(2,i) + phic .gt. 0.0d0 ) then
                  gradphi(1:2) = - gradphi(1:2)
               end if
      
               call matrixbd(sites(1:2,abs(W%e(ed))),sites(1:2,i),W%v(1:2,ed),gradphi,M,flag)
               if ( flag .ne. 0 ) return
  
               gval(1:2) = gval(1:2) + xi * mu * matmul( transpose( M(1:2,1:2) ),zeta(1:2) )
                     
               call matrixbd(sites(1:2,i),sites(1:2,abs(W%e(ed))),W%v(1:2,ed),gradphi,M,flag)
               if ( flag .ne. 0 ) return
               
               if ( gnnz + 2 .gt. gnnzmax ) then
                  write(*,*) 'In compJ4: Increase gnnzmax to at least ',gnnz+2
                  flag = - 1
                  return
               end if
          
               gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ed))-1, 2*abs(W%e(ed)) /)
               gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), zeta(1:2) )
               gnnz = gnnz + 2
               
            else if ( W%e(edprev) .lt. 0 ) then
               if ( debug ) then
                  write(*,*) 'We are in the case |P(v_E)|=2 (type 2). Neighbour cell is: ',abs(W%e(edprev))
               end if
                  
               call line_exp2imp_2d(W%v(1:2,ed),W%v(1:2,ednext),phia,phib,phic)
  
               gradphi(1:2) = (/ phia, phib /)
                  
               if ( gradphi(1) * sites(1,i) + gradphi(2) * sites(2,i) + phic .gt. 0.0d0 ) then
                  gradphi(1:2) = - gradphi(1:2)
               end if
      
               call matrixbd(sites(1:2,abs(W%e(edprev))),sites(1:2,i),W%v(1:2,ed),gradphi,M,flag)
               if ( flag .ne. 0 ) return
  
               gval(1:2) = gval(1:2) + xi * mu * matmul( transpose( M(1:2,1:2) ),zeta(1:2) )
                     
               call matrixbd(sites(1:2,i),sites(1:2,abs(W%e(edprev))),W%v(1:2,ed),gradphi,M,flag)
               if ( flag .ne. 0 ) return
                        
               if ( gnnz + 2 .gt. gnnzmax ) then
                  write(*,*) 'In compJ4: Increase gnnzmax to at least ',gnnz+2
                  flag = - 1
                  return
               end if
          
               gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(edprev))-1, 2*abs(W%e(edprev)) /)
               gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), zeta(1:2) )
               gnnz = gnnz + 2
            end if
                     
            ! Ending point w_E
                    
            if ( debug ) then
               write(*,*)
               write(*,*) 'Processing w_E'
               write(*,*)
            end if
                  
          ! zeta(1:2) = ( dot_product( normaled(1:2), tauedprev(1:2) ) / sizeed ) * normaled(1:2)
            zeta(1:2) = normaled(1:2) / sizeed
  
            if ( debug ) then
               write(*,*) 'zeta = ',zeta(1:2)
            end if
                  
            if ( W%e(ed) .lt. 0 .and. W%e(ednext) .lt. 0 ) then
               if ( debug ) then
                  write(*,*) 'We are in the case |P(w_E)|=3. ', &
                       'Neighbour cells are: ',abs(W%e(ed)),' and ',abs(W%e(ednext))
               end if
                  
               call matrixint(sites(1:2,abs(W%e(ed))),sites(1:2,abs(W%e(ednext))),sites(1:2,i),W%v(1:2,ednext),M,flag)
               if ( flag .ne. 0 ) return
                     
               gval(1:2) = gval(1:2) - xi * mu * matmul( transpose( M(1:2,1:2) ),zeta(1:2) )
               
               call matrixint(sites(1:2,i),sites(1:2,abs(W%e(ed))),sites(1:2,abs(W%e(ednext))),W%v(1:2,ednext),M,flag)
               if ( flag .ne. 0 ) return
                   
               if ( gnnz + 2 .gt. gnnzmax ) then
                  write(*,*) 'In compJ4: Increase gnnzmax to at least ',gnnz+2
                  flag = - 1
                  return
               end if
          
               gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ednext))-1, 2*abs(W%e(ednext)) /)
               gval(gnnz+1:gnnz+2) = - xi * mu * matmul( transpose( M(1:2,1:2) ), zeta(1:2) )
               gnnz = gnnz + 2
            
               call matrixint(sites(1:2,abs(W%e(ednext))),sites(1:2,i),sites(1:2,abs(W%e(ed))),W%v(1:2,ednext),M,flag)
               if ( flag .ne. 0 ) return
                   
               if ( gnnz + 2 .gt. gnnzmax ) then
                  write(*,*) 'In compJ4: Increase gnnzmax to at least ',gnnz+2
                  flag = - 1
                  return
               end if
          
               gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ed))-1, 2*abs(W%e(ed)) /)
               gval(gnnz+1:gnnz+2) = - xi * mu * matmul( transpose( M(1:2,1:2) ), zeta(1:2) )
               gnnz = gnnz + 2
            
            else if ( W%e(ed) .lt. 0 ) then 
               if ( debug ) then
                  write(*,*) 'We are in the case |P(w_E)|=2 (type 1). Neighbour cell is: ',abs(W%e(ed))
               end if
                  
               call line_exp2imp_2d(W%v(1:2,ednext),W%v(1:2,ednext+1),phia,phib,phic)
  
               gradphi(1:2) = (/ phia, phib /)
            
               if ( gradphi(1) * sites(1,i) + gradphi(2) * sites(2,i) + phic .gt. 0.0d0 ) then
                  gradphi(1:2) = - gradphi(1:2)
               end if
      
               call matrixbd(sites(1:2,abs(W%e(ed))),sites(1:2,i),W%v(1:2,ednext),gradphi,M,flag)
               if ( flag .ne. 0 ) return
  
               gval(1:2) = gval(1:2) - xi * mu * matmul( transpose( M(1:2,1:2) ),zeta(1:2) )
                     
               call matrixbd(sites(1:2,i),sites(1:2,abs(W%e(ed))),W%v(1:2,ednext),gradphi,M,flag)
               if ( flag .ne. 0 ) return
                     
               if ( gnnz + 2 .gt. gnnzmax ) then
                  write(*,*) 'In compJ4: Increase gnnzmax to at least ',gnnz+2
                  flag = - 1
                  return
               end if
          
               gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ed))-1, 2*abs(W%e(ed)) /)
               gval(gnnz+1:gnnz+2) = - xi * mu * matmul( transpose( M(1:2,1:2) ), zeta(1:2) )
               gnnz = gnnz + 2
               
            else if ( W%e(ednext) .lt. 0 ) then 
               if ( debug ) then
                  write(*,*) 'We are in the case |P(w_E)|=2 (type 2). Neighbour cell is: ',abs(W%e(ednext))
               end if
               
               call line_exp2imp_2d(W%v(1:2,ednext),W%v(1:2,ed),phia,phib,phic)
               
               gradphi(1:2) = (/ phia, phib /)
                  
               if ( gradphi(1) * sites(1,i) + gradphi(2) * sites(2,i) + phic .gt. 0.0d0 ) then
                  gradphi(1:2) = - gradphi(1:2)
               end if
      
               call matrixbd(sites(1:2,abs(W%e(ednext))),sites(1:2,i),W%v(1:2,ednext),gradphi,M,flag)
               if ( flag .ne. 0 ) return
  
               gval(1:2) = gval(1:2) - xi * mu * matmul( transpose( M(1:2,1:2) ),zeta(1:2) )
                     
               call matrixbd(sites(1:2,i),sites(1:2,abs(W%e(ednext))),W%v(1:2,ednext),gradphi,M,flag)
               if ( flag .ne. 0 ) return
                     
               if ( gnnz + 2 .gt. gnnzmax ) then
                  write(*,*) 'In compJ4: Increase gnnzmax to at least ',gnnz+2
                  flag = - 1
                  return
               end if
          
               gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ednext))-1, 2*abs(W%e(ednext)) /)
               gval(gnnz+1:gnnz+2) = - xi * mu * matmul( transpose( M(1:2,1:2) ), zeta(1:2) )
               gnnz = gnnz + 2
            end if
  
            ! Starting point point v_{\hat E}
            
            if ( edprev .eq. 1 ) then
               edpp = W%n - 1
            else
               edpp = edprev - 1
            end if
               
            if ( debug ) then
               write(*,*)
               write(*,*) 'Processing v_{\hat E}'
               write(*,*)
            end if
                  
          ! zeta(1:2) = ( dot_product( normaledprev(1:2), taued(1:2) ) / sizeedprev ) * normaledprev(1:2)
            zeta(1:2) = - normaledprev(1:2) / sizeedprev
  
            if ( debug ) then
               write(*,*) 'zeta = ',zeta(1:2)
            end if
                  
            if ( W%e(edprev) .lt. 0 .and. W%e(edpp) .lt. 0 ) then
               if ( debug ) then
                  write(*,*) 'We are in the case |P(w_E)|=3. ', &
                       'Neighbour cells are: ',abs(W%e(edprev)),' and ',abs(W%e(edpp))
               end if
                  
               call matrixint(sites(1:2,abs(W%e(edprev))),sites(1:2,abs(W%e(edpp))),sites(1:2,i),W%v(1:2,edprev),M,flag)
               if ( flag .ne. 0 ) return
                     
               gval(1:2) = gval(1:2) + xi * mu * matmul( transpose( M(1:2,1:2) ),zeta(1:2) )
               
               call matrixint(sites(1:2,i),sites(1:2,abs(W%e(edprev))),sites(1:2,abs(W%e(edpp))),W%v(1:2,edprev),M,flag)
               if ( flag .ne. 0 ) return
                   
               if ( gnnz + 2 .gt. gnnzmax ) then
                  write(*,*) 'In compJ4: Increase gnnzmax to at least ',gnnz+2
                  flag = - 1
                  return
               end if
          
               gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(edpp))-1, 2*abs(W%e(edpp)) /)
               gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), zeta(1:2) )
               gnnz = gnnz + 2
                     
               call matrixint(sites(1:2,abs(W%e(edpp))),sites(1:2,i),sites(1:2,abs(W%e(edprev))),W%v(1:2,edprev),M,flag)
               if ( flag .ne. 0 ) return
                   
               if ( gnnz + 2 .gt. gnnzmax ) then
                  write(*,*) 'In compJ4: Increase gnnzmax to at least ',gnnz+2
                  flag = - 1
                  return
               end if
          
               gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(edprev))-1, 2*abs(W%e(edprev)) /)
               gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), zeta(1:2) )
               gnnz = gnnz + 2
                     
            else if ( W%e(edprev) .lt. 0 ) then 
               if ( debug ) then
                  write(*,*) 'We are in the case |P(w_E)|=2 (type 1). Neighbour cell is: ',abs(W%e(edprev))
               end if
                  
               call line_exp2imp_2d(W%v(1:2,edpp),W%v(1:2,edprev),phia,phib,phic)
  
               gradphi(1:2) = (/ phia, phib /)
            
               if ( gradphi(1) * sites(1,i) + gradphi(2) * sites(2,i) + phic .gt. 0.0d0 ) then
                  gradphi(1:2) = - gradphi(1:2)
               end if
      
               call matrixbd(sites(1:2,abs(W%e(edprev))),sites(1:2,i),W%v(1:2,edprev),gradphi,M,flag)
               if ( flag .ne. 0 ) return
  
               gval(1:2) = gval(1:2) + xi * mu * matmul( transpose( M(1:2,1:2) ),zeta(1:2) )
                  
               call matrixbd(sites(1:2,i),sites(1:2,abs(W%e(edprev))),W%v(1:2,edprev),gradphi,M,flag)
               if ( flag .ne. 0 ) return
                     
               if ( gnnz + 2 .gt. gnnzmax ) then
                  write(*,*) 'In compJ4: Increase gnnzmax to at least ',gnnz+2
                  flag = - 1
                  return
               end if
          
               gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(edprev))-1, 2*abs(W%e(edprev)) /)
               gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), zeta(1:2) )
               gnnz = gnnz + 2
            
            else if ( W%e(edpp) .lt. 0 ) then 
               if ( debug ) then
                  write(*,*) 'We are in the case |P(w_E)|=2 (type 2). Neighbour cell is: ',abs(W%e(edpp))
               end if
               
               call line_exp2imp_2d(W%v(1:2,edprev),W%v(1:2,ed),phia,phib,phic)
               
               gradphi(1:2) = (/ phia, phib /)
                  
               if ( gradphi(1) * sites(1,i) + gradphi(2) * sites(2,i) + phic .gt. 0.0d0 ) then
                  gradphi(1:2) = - gradphi(1:2)
               end if
      
               call matrixbd(sites(1:2,abs(W%e(edpp))),sites(1:2,i),W%v(1:2,edprev),gradphi,M,flag)
               if ( flag .ne. 0 ) return
  
               gval(1:2) = gval(1:2) + xi * mu * matmul( transpose( M(1:2,1:2) ),zeta(1:2) )
               
               call matrixbd(sites(1:2,i),sites(1:2,abs(W%e(edpp))),W%v(1:2,edprev),gradphi,M,flag)
               if ( flag .ne. 0 ) return
                     
               if ( gnnz + 2 .gt. gnnzmax ) then
                  write(*,*) 'In compJ4: Increase gnnzmax to at least ',gnnz+2
                  flag = - 1
                  return
               end if
          
               gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(edpp))-1, 2*abs(W%e(edpp)) /)
               gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), zeta(1:2) )
               gnnz = gnnz + 2
            end if
         else
            if ( debug ) then
               write(*,*)
               write(*,*) 'The inner vertex of this edge has a fixed angle; so it will  be ignored.'
               write(*,*)
            end if
         end if
      end do
         
   end subroutine compJ3
  

   ! ******************************************************************
   ! ******************************************************************
   
   subroutine fgcompJ4(x,fden,gden,compgrad,pdataptr)

      implicit none

      ! SCALAR ARGUMENTS
      logical, intent(in) :: compgrad
      real(kind=8), intent(out) :: fden
      type(c_ptr), optional, intent(in) :: pdataptr

      ! ARRAY ARGUMENTS
      real(kind=8), intent(in) :: x(2)
      real(kind=8), intent(out) :: gden(2)

      ! LOCAL SCALARS
      type(pdata_type), pointer :: pdata
      real(kind=8) :: delta,param,a,b,scale,z1,z2,faux
      real(kind=8), parameter :: pi = 4.0d0 * atan( 1.0d0 )

      ! LOCAL ARRAY
      real(kind=8) :: z(2)

      call c_f_pointer(pdataptr,pdata)
      


      !! Function psi_1
      ! z(1:2) = -pdata%c(1:2)/5.0d0 + 0.40d0 * x(1:2)
      ! b = 0.25d0
      ! a = 0.0625d0 * (2.0d0 - b - 0.5625d0)
      ! fden = a * ( (z(2) - ( 0.25d0 * z(1) ) ** 2 ) ** 2 + ( 0.25d0 * z(1) - 1.0d0 ) ** 2 ) + b
      ! if ( compgrad ) then
      !    gden(1) = a * ( -0.1d0 * z(1) * ( z(2) - ( 0.25d0 * z(1) ) ** 2 ) + 0.2d0 * (0.25d0 * z(1) - 1.0d0) )
      !    gden(2) = a * ( 0.8d0 * ( z(2) - ( 0.25d0 * z(1) ) ** 2 ) ) 
      ! end if

      !! Function psi_2
      delta = 0.5d0*sqrt( 1000.d0 )
      param = 2.9d0/delta**2
      fden = 0.1d0 + param*( x(2) - 0.6d0*delta*sin( 2.0d0*pi*x(1)/(2.0d0*delta) ) - delta ) ** 2
      if ( compgrad ) then
          gden(1) = 2.0d0 * param * ( x(2) - 0.6d0*delta*sin( 2.0d0*pi*x(1)/(2.0d0*delta) ) ) &
                   * ( -0.6*2.0d0*pi*cos( 2.0d0*pi*x(1)/(2.0d0*delta) )/(2.0d0*delta) )
          gden(2) = 2.0d0 * param * ( x(2) - 0.6d0*delta*sin( 2.0d0*pi*x(1)/(2.0d0*delta) ) )
      end if

      !! Function 3 psi_3

      ! fden = 0.01d0 + 20.0d0 * sum( ( x(1:2) - pdata%J4c(1:2) ) ** 2 ) / pdata%J4r ** 2

      ! if ( compgrad ) gden(1:2) = 40.0d0 * ( x(1:2) - pdata%J4c(1:2) ) / pdata%J4r ** 2
      


   end subroutine fgcompJ4

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

   subroutine centroid(W,cent,volW)
      
      implicit none

      ! SCALAR ARGUMENTS
      type(polygon), intent(in) :: W
      real(kind=8), intent(out) :: volW

      ! ARRAY ARGUMENTS
      real(kind=8), intent(out) :: cent(2)

      ! LOCAL SCALARS
      integer :: j
      real(kind=8) :: re

      call polygon_area_2d(W%n-1,W%v(1:2,1:W%n-1),volW)
      
      cent = 0.0d0

      do j=1,W%n-1
         if ( j .eq. W%n-1 ) then
            re = W%v(1,j)*W%v(2,1) - W%v(1,1)*W%v(2,j)
            cent = cent + re*(W%v(1:2,j) + W%v(1:2,1))
         else
            re = W%v(1,j)*W%v(2,j+1) - W%v(1,j+1)*W%v(2,j)
            cent = cent + re*(W%v(1:2,j) + W%v(1:2,j+1)) 
         end if
      end do
      cent = cent/(6.0d0*volW)

   end subroutine centroid

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

   subroutine det(x,y,z,vdet)

        implicit none

        ! SACALAR ARGUMENTS
        real(kind=8), intent(out) :: vdet

        ! ARRAY ARGUMENTS
        real(kind=8), intent(in) :: x(2),y(2),z(2)

        vdet = ( x(1) - y(1) ) * ( z(2) - y(2) ) - ( x(2) - y(2) ) * ( z(1) - y(1) )

   end subroutine det
  
    
    ! ******************************************************************
    ! ******************************************************************

   subroutine matrixint(x,y,z,v,M,flag)

      ! It computes the matrix used in the calculation of the gradient
      ! when v belongs to T_int.
            
      implicit none

      ! SCALAR ARGUMENTS
      integer, intent(inout) :: flag

      ! ARRAY ARGUMENTS
      real(kind=8), intent(in) :: x(2),y(2),z(2),v(2)
      real(kind=8), intent(out) :: M(2,2)

      ! LOCAL SCALARS
      real(kind=8) :: det

      ! LOCAL ARRAYS
      real(kind=8) :: tmp(2)

      det = ( y(1) - x(1) ) * ( z(2) - x(2) ) - ( y(2) - x(2) ) * ( z(1) - x(1) )

      if ( abs( det ) .lt. 1.0d+02 * epsilon( 1.0d0 ) ) then
         det = sign( 1.0d+02 * epsilon( 1.0d0 ), det )
      end if

      ! This is to avoid a warning of the compiler that says flag is not
      ! being used.
      if ( .false. ) flag = - 1

      !!$    if ( det .eq. 0.0d0 ) then
      !!$       write(*,*) 'In matrixint: Division by zero in the calculation of matrix M.'
      !!$       flag = - 1
      !!$       return
      !!$       
      !!$    else if ( abs(det) .lt. 1.0d-14 ) then
      !!$       write(*,*) 'In matrixint: Warning: Determinant is too small.'
      !!$    end if
            
      tmp(1:2) = (/ y(2) - x(2), x(1) - y(1) /)

      M(1:2,1) = tmp(1:2) * ( v(1) - z(1) )
      M(1:2,2) = tmp(1:2) * ( v(2) - z(2) )
      M(1:2,1:2) = M(1:2,1:2) / det

   end subroutine matrixint

   subroutine matrixbd(y,x,v,gradphi,M,flag)
      
      ! It computes the matrix used in the calculation of the gradient
      ! when v belongs to T_bd.
          
      implicit none
          
      ! SCALAR ARGUMENTS
      integer, intent(inout) :: flag
      
      ! ARRAY ARGUMENTS
      real(kind=8), intent(in) :: gradphi(2),x(2),y(2),v(2)
      real(kind=8), intent(out) :: M(2,2)
          
      ! LOCAL SCALARS
      real(kind=8) :: det
  
      ! LOCAL ARRAYS
      real(kind=8) :: tmp(2)
      
      tmp(1:2) = (/ - gradphi(2), gradphi(1) /)
      
      det = - dot_product( tmp(1:2), y(1:2) - x(1:2) )
  
      if ( abs( det ) .lt. 1.0d+02 * epsilon( 1.0d0 ) ) then
         det = sign( 1.0d+02 * epsilon( 1.0d0 ), det )
      end if
  
      ! This is to avoid a warning of the compiler that says flag is not
      ! being used.
      if ( .false. ) flag = - 1
      
  !!$    if ( det .eq. 0.0d0 ) then
  !!$       write(*,*) 'In matrixbd: Division by zero in the calculation of matrix M.'
  !!$       flag = - 1
  !!$       return
  !!$       
  !!$    else if ( abs(det) .lt. 1.0d-14 ) then
  !!$       write(*,*) 'In matrixbd: Warning: Determinant is too small.'
  !!$    end if
          
      ! Outer product of (/ - gradphi(2), gradphi(1) /) with p(1:2) - x(1:2)
      M(1:2,1) = tmp(1:2) * ( v(1) - x(1) ) 
      M(1:2,2) = tmp(1:2) * ( v(2) - x(2) )
      M(1:2,1:2) = - M(1:2,1:2) / det
      
   end subroutine matrixbd

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

   subroutine drawsol(n,x,filename,pdataptr)
      
      implicit none
    
      ! SCALAR ARGUMENTS
      integer, intent(in) :: n
      type(c_ptr), optional, intent(in) :: pdataptr
        
      ! ARRAY ARGUMENTS
      real(kind=8), intent(in) :: x(n)
      character(len=80) :: filename
  
      ! LOCAL SCALARS
      logical :: debug
      integer :: allocerr,flag,i,ierror,j,k,maxnv,nsites,status,ntriag,v
      type(pdata_type), pointer :: pdata
      type(VoronoiCell) :: Vi
      type(Polygon) :: Wij,Pj,frame,Viframe
      
      ! LOCAL ARRAYS
      integer :: indices(n/2),sitetv(3*2*(n/2)),start((n/2)+1),til(3,2*(n/2)),tnbr(3,2*(n/2))
      real(kind=8) :: f(0:3,n/2),sites(2,n/2),vorxy(2,2*(n/2))
  
      debug = .false.
      call compallJ(n,x,f,debug,flag,pdataptr)
      if ( flag .ne. 0 ) then
         write(*,*) 'In drawsol, compallJ returned a non-null flag!'
         return
      end if
      ! write(*,*) 'ftmpscompallJ= ', f(3,1:n/2)
      ! write(*,*) 'J4compall = ', sum( f(3,1:n/2) ** 2 ) / (n/2) 
      call c_f_pointer(pdataptr,pdata)
  
      nsites = n / 2
      indices(1:nsites) = (/ (i,i=1,nsites) /)
      sites(1:2,1:nsites) = reshape( x(1:n), (/ 2, nsites /) )
    
      if ( nsites .le. 1 ) then
         write(*,*) 'In drawsol: nsites must be at least 2!'
         stop
      end if
  
      ! Compute Delaunay triangulation
  
      call dtris2(nsites,sites,indices,ntriag,til,tnbr,ierror)
  
      if ( ierror .ne. 0 .and. ierror .ne. 225 ) then
         write(*,*) 'In drawsol: dtris2 returned an error different from 0 and 225: ',ierror
         stop
      end if
    
      call voronoi_cell_patch(nsites,ntriag,til,start,sitetv)
      
      ! Compute vertices of the Voronoi diagram
            
      if ( ntriag .gt. 2 * nsites ) then
         write(*,*) 'In drawsol: 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 drawsol: triangle_circumcenter_2d returned NaN.'
            write(*,*) 'Triangle vertices: ',sites(1:2,til(1:3,k))
            stop
         end if
      end do
  
      maxnv = maxval( pdata%nvpols(1:pdata%npols) )
        
      allocate(Pj%v(2,maxnv),Pj%e(maxnv-1),stat=allocerr)
      if ( allocerr .ne. 0 ) then
         write(*,*) 'In drawsol: Allocation error.'
         stop
      end if
        
      open(unit=10,file=trim(filename))                           
  
      write(10,10) 1,pdata%drawscale
  
      ! Drawing Delaunay triangulation
      
      !! For the paper it is not necessary to draw the Delaunay triangulation
      do k = 1,ntriag
         do v = 1,3
            write(10,20) sites(1:2,til(v,k)),sites(1:2,til(mod(v,3)+1,k)),'teagreen'
         end do
      end do
  
      ! Frame
      
      frame%n = 5
      frame%deg = .false.
      
      allocate(frame%v(2,frame%n),frame%e(frame%n-1),stat=allocerr)
      if ( allocerr .ne. 0 ) then
         write(*,*) 'In drawsol: Allocation error.'
         stop
      end if
      
      frame%v(1,1:5) = (/ pdata%xmin, pdata%xmax, pdata%xmax, pdata%xmin, pdata%xmin /)
      frame%v(2,1:5) = (/ pdata%ymin, pdata%ymin, pdata%ymax, pdata%ymax, pdata%ymin /)
      frame%e(1:4) = (/ 0, 0, 0, 0 /)  
          
      ! Drawing Voronoi diagram
      
      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 drawsol: Error in voronoi_cell, status = ',status
            stop
         end if
  
         ! This is to draw the Voronoi cell restricted to the frame.
         
         if ( nsites .ge. 3 .and. ierror .eq. 0 ) then
            call VorCellInterConvPol(sites(1:2,i),frame,Vi,Viframe)
         else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then
            call VorCellInterConvPol_Collinear(i,nsites,sites,frame,Viframe)
         end if
  
         if ( .not. Viframe%deg ) then
            do v = 1,Viframe%n - 1
               write(10,20) Viframe%v(1:2,v),Viframe%v(1:2,v+1),'jonquil'
            end do
         end if
  
         call polygon_destroy(Viframe)
  
         ! This is to draw the Voronoi cell restricted to A.
  
         do j = 1,pdata%npols
            Pj%n = pdata%nvpols(j)
            Pj%deg = .false.
  
            Pj%v(1,1:Pj%n) = pdata%xpol(j,1:Pj%n)
            Pj%v(2,1:Pj%n) = pdata%ypol(j,1:Pj%n)
            Pj%e(1:Pj%n-1) = pdata%poledges(j,1:Pj%n-1)
  
            if ( nsites .ge. 3 .and. ierror .eq. 0 ) then
               call VorCellInterConvPol(sites(1:2,i),Pj,Vi,Wij)
            else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then
               call VorCellInterConvPol_Collinear(i,nsites,sites,Pj,Wij)
            end if
  
            if ( .not. Wij%deg ) then
               do v = 1,Wij%n - 1
                  if ( Wij%e(v) .lt. 1 ) then
                     write(10,20) Wij%v(1:2,v),Wij%v(1:2,v+1),'desire'
                  end if
               end do
            end if
  
            call polygon_destroy(Wij)
         end do
  
         ! This is to draw the sites
         write(10,40) sites(1:2,i)
         
         call voronoi_cell_destroy(Vi)
      end do
        
      deallocate(frame%v,frame%e,stat=allocerr)
      if ( allocerr .ne. 0 ) then
         write(*,*) 'In drawsol: Deallocation error.'
         stop
      end if
          
      ! Drawing A
      
      do j = 1,pdata%npols
         Pj%n = pdata%nvpols(j)
         Pj%deg = .false.          
  
         Pj%v(1,1:Pj%n) = pdata%xpol(j,1:Pj%n)
         Pj%v(2,1:Pj%n) = pdata%ypol(j,1:Pj%n)
         Pj%e(1:Pj%n-1) = pdata%poledges(j,1:Pj%n-1)
  
         do k = 1,Pj%n - 1
            if ( Pj%e(k) .eq. 0 ) then
               write(10,20) Pj%v(1:2,k),Pj%v(1:2,k+1),'black'
            end if
         end do
      end do
  
      ! Close figure
      
      write(10,50)
      
      write(10,*) 'Value of J4i (sum/ave/min/max) = ',sum(abs(f(3,1:n/2))),sum(abs(f(3,1:n/2)))/(n/2), &
           minval(abs(f(3,1:n/2))),maxval(abs(f(3,1:n/2))),' J4 itself = ',sum(f(3,1:n/2)**2)/(n/2)
  
      close(10)
  
      deallocate(Pj%v,Pj%e,stat=allocerr)
      if ( allocerr .ne. 0 ) then
         write(*,*) 'In drawsol: Deallocation error.'
         stop
      end if
  
      ! NON-EXECUTABLE STATEMENTS
      
  10  format('prologues := 3;',/, &
           'outputtemplate := "%j.mps";',/, &
           'input mpcolornames;',/, &
           'beginfig(',i4,');',/, &
           '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 = ',f20.10,'cm;')
  20  format('draw ( ',f20.10,'u,',f20.10,'u)--(',f20.10,'u,',f20.10,'u) withcolor',a12,';')
  40  format('drawdot(',f20.10,'u,',f20.10,'u);') 
  50  format('endfig;',/,'end;')
  
end subroutine drawsol

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

subroutine drawsolJ1(n,x,filename,pdataptr)

   implicit none
    
   ! SCALAR ARGUMENTS
   integer, intent(in) :: n
   type(c_ptr), optional, intent(in) :: pdataptr
     
   ! ARRAY ARGUMENTS
   real(kind=8), intent(in) :: x(n)
   character(len=80) :: filename

   ! LOCAL SCALARS
   integer :: allocerr,i,ierror,j,k,maxnv,nsites,status,ntriag,v
   type(pdata_type), pointer :: pdata
   type(VoronoiCell) :: Vi
   type(Polygon) :: Wij,Pj,frame,Viframe
   
   ! LOCAL ARRAYS
   integer :: indices(n/2),sitetv(3*2*(n/2)),start((n/2)+1),til(3,2*(n/2)),tnbr(3,2*(n/2))
   real(kind=8) :: sites(2,n/2),vorxy(2,2*(n/2))

   integer :: flag,find
   logical :: debug
   real(kind=8) :: ctol
   integer :: cellcolorJ1(n/2)
   real(kind=8) :: f(0:3,n/2)
   character(len=80) :: paintcolor

   call c_f_pointer(pdataptr,pdata) 

   ! The smaller the ctol for which the cell has positive f, the darker the color.
   cellcolorJ1(1:n/2) = -1       ! Means no color

   debug = .false.
   call compallJ(n,x,f,debug,flag,pdataptr)
   if ( flag .ne. 0 ) then
      write(*,*) 'In drawsol, compallJ returned a non-null flag!'
      return
   end if

   where ( abs( f(0,1:n/2) ) .gt. 1.0d-03 )
      cellcolorJ1(1:n/2) = 1
   end where

   call c_f_pointer(pdataptr,pdata)
  
   nsites = n / 2
   indices(1:nsites) = (/ (i,i=1,nsites) /)
   sites(1:2,1:nsites) = reshape( x(1:n), (/ 2, nsites /) )
 
   if ( nsites .le. 1 ) then
      write(*,*) 'In drawsol: nsites must be at least 2!'
      stop
   end if

   ! Compute Delaunay triangulation

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

   if ( ierror .ne. 0 .and. ierror .ne. 225 ) then
      write(*,*) 'In drawsol: dtris2 returned an error different from 0 and 225: ',ierror
      stop
   end if
 
   call voronoi_cell_patch(nsites,ntriag,til,start,sitetv)
   
   ! Compute vertices of the Voronoi diagram
         
   if ( ntriag .gt. 2 * nsites ) then
      write(*,*) 'In drawsol: 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 drawsol: triangle_circumcenter_2d returned NaN.'
         write(*,*) 'Triangle vertices: ',sites(1:2,til(1:3,k))
         stop
      end if
   end do

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

   allocate(Pj%v(2,maxnv),Pj%e(maxnv-1),stat=allocerr)
   if ( allocerr .ne. 0 ) then
      write(*,*) 'In drawsol: Allocation error.'
      stop
   end if
     
   open(unit=10,file=trim(filename))                           

   write(10,10) 1,pdata%drawscale

   ! Drawing Delaunay triangulation
   
   do k = 1,ntriag
      do v = 1,3
         write(10,20) sites(1:2,til(v,k)),sites(1:2,til(mod(v,3)+1,k)),'teagreen'
      end do
   end do

   ! Frame
   
   frame%n = 5
   frame%deg = .false.
   
   allocate(frame%v(2,frame%n),frame%e(frame%n-1),stat=allocerr)
   if ( allocerr .ne. 0 ) then
      write(*,*) 'In drawsol: Allocation error.'
      stop
   end if

   frame%v(1,1:5) = (/ pdata%xmin, pdata%xmax, pdata%xmax, pdata%xmin, pdata%xmin /)
   frame%v(2,1:5) = (/ pdata%ymin, pdata%ymin, pdata%ymax, pdata%ymax, pdata%ymin /)
   frame%e(1:4) = (/ 0, 0, 0, 0 /)  
          
   ! Drawing Voronoi diagram
      
   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 drawsol: Error in voronoi_cell, status = ',status
         stop
      end if
  
      ! This is to draw the Voronoi cell restricted to the frame.
      
      if ( nsites .ge. 3 .and. ierror .eq. 0 ) then
         call VorCellInterConvPol(sites(1:2,i),frame,Vi,Viframe)
      else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then
         call VorCellInterConvPol_Collinear(i,nsites,sites,frame,Viframe)
      end if

      if ( .not. Viframe%deg ) then
         do v = 1,Viframe%n - 1
            write(10,20) Viframe%v(1:2,v),Viframe%v(1:2,v+1),'jonquil'
         end do
      end if

      call polygon_destroy(Viframe)

         ! This is to paint the Voronoi cell, restricted to A, if
      ! desired conditions does not hold.

      if ( cellcolorJ1(i) .ne. - 1 ) then
         if ( cellcolorJ1(i) .eq. 1 ) then
            paintcolor = 'dollarbill'
         end if

         do j = 1,pdata%npols
            Pj%n = pdata%nvpols(j)
            Pj%deg = .false.         
            
            Pj%v(1,1:Pj%n) = pdata%xpol(j,1:Pj%n)
            Pj%v(2,1:Pj%n) = pdata%ypol(j,1:Pj%n)
            Pj%e(1:Pj%n-1)  = pdata%poledges(j,1:Pj%n-1)

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

            if ( .not. Wij%deg ) then
               do v = 1,Wij%n - 1
                  if ( v .eq. 1 ) then
                     write(10,30) Wij%v(1:2,v),Wij%v(1:2,v+1)
                  else if ( v .gt. 1 .and. v .lt. Wij%n - 1 ) then 
                     write(10,31) Wij%v(1:2,v),Wij%v(1:2,v+1)
                  else
                     write(10,32) Wij%v(1:2,v),trim(paintcolor)
                  end if
               end do
            end if

            call polygon_destroy(Wij)
         end do
      end if

      ! This is to draw the Voronoi cell restricted to A.

      do j = 1,pdata%npols
         Pj%n = pdata%nvpols(j)
         Pj%deg = .false.

         Pj%v(1,1:Pj%n) = pdata%xpol(j,1:Pj%n)
         Pj%v(2,1:Pj%n) = pdata%ypol(j,1:Pj%n)
         Pj%e(1:Pj%n-1) = pdata%poledges(j,1:Pj%n-1)

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

         if ( .not. Wij%deg ) then
            do v = 1,Wij%n - 1
               if ( Wij%e(v) .lt. 1 ) then
                  write(10,20) Wij%v(1:2,v),Wij%v(1:2,v+1),'desire'
               end if
            end do
         end if

         call polygon_destroy(Wij)
      end do



   ! This is to draw the sites
   write(10,40) sites(1:2,i)
   
   call voronoi_cell_destroy(Vi)
   end do

   deallocate(frame%v,frame%e,stat=allocerr)
   if ( allocerr .ne. 0 ) then
      write(*,*) 'In drawsol: Deallocation error.'
      stop
   end if
       
   ! Drawing A
   
   do j = 1,pdata%npols
      Pj%n = pdata%nvpols(j)
      Pj%deg = .false.          

      Pj%v(1,1:Pj%n) = pdata%xpol(j,1:Pj%n)
      Pj%v(2,1:Pj%n) = pdata%ypol(j,1:Pj%n)
      Pj%e(1:Pj%n-1) = pdata%poledges(j,1:Pj%n-1)

      do k = 1,Pj%n - 1
         if ( Pj%e(k) .eq. 0 ) then
            write(10,20) Pj%v(1:2,k),Pj%v(1:2,k+1),'black'
         end if
      end do
   end do

   ! Close figure

   write(10,50)

   close(10)

   deallocate(Pj%v,Pj%e,stat=allocerr)
   if ( allocerr .ne. 0 ) then
      write(*,*) 'In drawsol: Deallocation error.'
      stop
   end if

   ! NON-EXECUTABLE STATEMENTS
   
10  format('prologues := 3;',/, &
        'outputtemplate := "%j.mps";',/, &
        'input mpcolornames;',/, &
        'beginfig(',i4,');',/, &
        '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 = ',f20.10,'cm;')
20  format('draw ( ',f20.10,'u,',f20.10,'u)--(',f20.10,'u,',f20.10,'u) withcolor',a12,';')
30  format('fill ((',f20.10,'u,',f20.10,'u)--(',f20.10,'u,',f20.10,'u)--')
31  format('      (',f20.10,'u,',f20.10,'u)--(',f20.10,'u,',f20.10,'u)--')
32  format('      (',f20.10,'u,',f20.10,'u)--cycle) withcolor ',a30,';')              
40  format('drawdot(',f20.10,'u,',f20.10,'u);') 
50  format('endfig;',/,'end;')

end subroutine drawsolJ1

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

subroutine drawsolJ2(n,x,filename,pdataptr)

   implicit none
    
   ! SCALAR ARGUMENTS
   integer, intent(in) :: n
   type(c_ptr), optional, intent(in) :: pdataptr
     
   ! ARRAY ARGUMENTS
   real(kind=8), intent(in) :: x(n)
   character(len=80) :: filename

   ! LOCAL SCALARS
   integer :: allocerr,i,ierror,j,k,maxnv,nsites,status,ntriag,v
   type(pdata_type), pointer :: pdata
   type(VoronoiCell) :: Vi
   type(Polygon) :: Wij,Pj,frame,Viframe
   
   ! LOCAL ARRAYS
   integer :: indices(n/2),sitetv(3*2*(n/2)),start((n/2)+1),til(3,2*(n/2)),tnbr(3,2*(n/2))
   real(kind=8) :: sites(2,n/2),vorxy(2,2*(n/2))

   integer :: flag,find
   logical :: debug
   real(kind=8) :: ctol, J2ctoltmp
   integer :: cellcolorJ2(n/2)
   real(kind=8) :: f(0:3,n/2)
   character(len=80) :: paintcolor

   call c_f_pointer(pdataptr,pdata) 

   J2ctoltmp = pdata%J2ctol
   
   ! The smaller the ctol for which the cell has positive f, the darker the color.
   cellcolorJ2(1:n/2) = -1       ! Means no color

   do k = 8,1,-1
      ctol = 0.1d0 * k
      pdata%J2ctol = ctol
      debug = .false.
      call compallJ(n,x,f,debug,flag,pdataptr)
      if ( flag .ne. 0 ) then
         write(*,*) 'In drawsol, compallJ returned a non-null flag!'
         return
      end if

      where ( f(1,1:n/2) .gt. 1.0d-04 )
         cellcolorJ2(1:n/2) = k
      end where
   end do

   pdata%J2ctol = J2ctoltmp

   call c_f_pointer(pdataptr,pdata)
  
   nsites = n / 2
   indices(1:nsites) = (/ (i,i=1,nsites) /)
   sites(1:2,1:nsites) = reshape( x(1:n), (/ 2, nsites /) )
 
   if ( nsites .le. 1 ) then
      write(*,*) 'In drawsol: nsites must be at least 2!'
      stop
   end if

   ! Compute Delaunay triangulation

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

   if ( ierror .ne. 0 .and. ierror .ne. 225 ) then
      write(*,*) 'In drawsol: dtris2 returned an error different from 0 and 225: ',ierror
      stop
   end if
 
   call voronoi_cell_patch(nsites,ntriag,til,start,sitetv)
   
   ! Compute vertices of the Voronoi diagram
         
   if ( ntriag .gt. 2 * nsites ) then
      write(*,*) 'In drawsol: 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 drawsol: triangle_circumcenter_2d returned NaN.'
         write(*,*) 'Triangle vertices: ',sites(1:2,til(1:3,k))
         stop
      end if
   end do

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

   allocate(Pj%v(2,maxnv),Pj%e(maxnv-1),stat=allocerr)
   if ( allocerr .ne. 0 ) then
      write(*,*) 'In drawsol: Allocation error.'
      stop
   end if
     
   open(unit=10,file=trim(filename))                           

   write(10,10) 1,pdata%drawscale

   ! Drawing Delaunay triangulation
   
   do k = 1,ntriag
      do v = 1,3
         write(10,20) sites(1:2,til(v,k)),sites(1:2,til(mod(v,3)+1,k)),'teagreen'
      end do
   end do

   ! Frame
   
   frame%n = 5
   frame%deg = .false.
   
   allocate(frame%v(2,frame%n),frame%e(frame%n-1),stat=allocerr)
   if ( allocerr .ne. 0 ) then
      write(*,*) 'In drawsol: Allocation error.'
      stop
   end if

   frame%v(1,1:5) = (/ pdata%xmin, pdata%xmax, pdata%xmax, pdata%xmin, pdata%xmin /)
   frame%v(2,1:5) = (/ pdata%ymin, pdata%ymin, pdata%ymax, pdata%ymax, pdata%ymin /)
   frame%e(1:4) = (/ 0, 0, 0, 0 /)  
          
      ! Drawing Voronoi diagram
      
      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 drawsol: Error in voronoi_cell, status = ',status
            stop
         end if
  
         ! This is to draw the Voronoi cell restricted to the frame.
         
         if ( nsites .ge. 3 .and. ierror .eq. 0 ) then
            call VorCellInterConvPol(sites(1:2,i),frame,Vi,Viframe)
         else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then
            call VorCellInterConvPol_Collinear(i,nsites,sites,frame,Viframe)
         end if
  
         if ( .not. Viframe%deg ) then
            do v = 1,Viframe%n - 1
               write(10,20) Viframe%v(1:2,v),Viframe%v(1:2,v+1),'jonquil'
            end do
         end if
  
         call polygon_destroy(Viframe)

          ! This is to paint the Voronoi cell, restricted to A, if
         ! desired conditions does not hold.
  
         if ( cellcolorJ2(i) .ne. - 1 ) then
            if ( cellcolorJ2(i) .eq. 8 ) then
               paintcolor = 'blueaa'
            else if ( cellcolorJ2(i) .eq. 7 ) then
               paintcolor = 'bluebb'
            else if ( cellcolorJ2(i) .eq. 6 ) then
               paintcolor = 'bluecc'
            else if ( cellcolorJ2(i) .eq. 5 ) then
               paintcolor = 'bluedd'
            else if ( cellcolorJ2(i) .eq. 4 ) then
               paintcolor = 'blueee'
            else if ( cellcolorJ2(i) .eq. 3 ) then
               paintcolor = 'bluef'
            else if ( cellcolorJ2(i) .eq. 2 ) then
               paintcolor = 'blueg'
            else if ( cellcolorJ2(i) .eq. 1 ) then 
               paintcolor = 'blueh' 
            end if
            


            do j = 1,pdata%npols
               Pj%n = pdata%nvpols(j)
               Pj%deg = .false.         
               
               Pj%v(1,1:Pj%n) = pdata%xpol(j,1:Pj%n)
               Pj%v(2,1:Pj%n) = pdata%ypol(j,1:Pj%n)
               Pj%e(1:Pj%n-1)  = pdata%poledges(j,1:Pj%n-1)
  
               if ( nsites .ge. 3 .and. ierror .eq. 0 ) then
                  call VorCellInterConvPol(sites(1:2,i),Pj,Vi,Wij)
               else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then
                  call VorCellInterConvPol_Collinear(i,nsites,sites,Pj,Wij)
               end if
  
               if ( .not. Wij%deg ) then
                  do v = 1,Wij%n - 1
                     if ( v .eq. 1 ) then
                        write(10,30) Wij%v(1:2,v),Wij%v(1:2,v+1)
                     else if ( v .gt. 1 .and. v .lt. Wij%n - 1 ) then 
                        write(10,31) Wij%v(1:2,v),Wij%v(1:2,v+1)
                     else
                        write(10,32) Wij%v(1:2,v),trim(paintcolor)
                     end if
                  end do
               end if
  
               call polygon_destroy(Wij)
            end do
         end if

         ! This is to draw the Voronoi cell restricted to A.
  
         do j = 1,pdata%npols
            Pj%n = pdata%nvpols(j)
            Pj%deg = .false.
  
            Pj%v(1,1:Pj%n) = pdata%xpol(j,1:Pj%n)
            Pj%v(2,1:Pj%n) = pdata%ypol(j,1:Pj%n)
            Pj%e(1:Pj%n-1) = pdata%poledges(j,1:Pj%n-1)
  
            if ( nsites .ge. 3 .and. ierror .eq. 0 ) then
               call VorCellInterConvPol(sites(1:2,i),Pj,Vi,Wij)
            else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then
               call VorCellInterConvPol_Collinear(i,nsites,sites,Pj,Wij)
            end if
  
            if ( .not. Wij%deg ) then
               do v = 1,Wij%n - 1
                  if ( Wij%e(v) .lt. 1 ) then
                     write(10,20) Wij%v(1:2,v),Wij%v(1:2,v+1),'desire'
                  end if
               end do
            end if
  
            call polygon_destroy(Wij)
         end do



      ! This is to draw the sites
      write(10,40) sites(1:2,i)
      
      call voronoi_cell_destroy(Vi)
   end do

   deallocate(frame%v,frame%e,stat=allocerr)
   if ( allocerr .ne. 0 ) then
      write(*,*) 'In drawsol: Deallocation error.'
      stop
   end if
       
   ! Drawing A
   
   do j = 1,pdata%npols
      Pj%n = pdata%nvpols(j)
      Pj%deg = .false.          

      Pj%v(1,1:Pj%n) = pdata%xpol(j,1:Pj%n)
      Pj%v(2,1:Pj%n) = pdata%ypol(j,1:Pj%n)
      Pj%e(1:Pj%n-1) = pdata%poledges(j,1:Pj%n-1)

      do k = 1,Pj%n - 1
         if ( Pj%e(k) .eq. 0 ) then
            write(10,20) Pj%v(1:2,k),Pj%v(1:2,k+1),'black'
         end if
      end do
   end do

   ! Close figure

   write(10,50)

   close(10)

   deallocate(Pj%v,Pj%e,stat=allocerr)
   if ( allocerr .ne. 0 ) then
      write(*,*) 'In drawsol: Deallocation error.'
      stop
   end if

   ! NON-EXECUTABLE STATEMENTS
   
10  format('prologues := 3;',/, &
        'outputtemplate := "%j.mps";',/, &
        'input mpcolornames;',/, &
        'beginfig(',i4,');',/, &
        'color skyblue, dollarbill, desire, jonquil, royalblue, teagreen, ', &
        'orange, plum, bluea, blueb, bluec, blued, bluee, blueaa, bluebb, bluecc, ', &
        'bluedd, blueee, bluef, blueg, blueh, bluei, 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);',/, &
        'blueaa       = (  0/255,255/255,230/255);',/, &
        'bluebb       = (  0/255,235/255,230/255);',/, &
        'bluecc       = (  0/255,205/255,230/255);',/, &
        'bluedd       = (  0/255,175/255,230/255);',/, &
        'blueee       = (  0/255,145/255,230/255);',/, &
        'bluef        = (  0/255,115/255,230/255);',/, &
        'blueg        = (  0/255, 85/255,230/255);',/, &
        'blueh        = (  0/255, 55/255,230/255);',/, &
        'bluei        = (  0/255, 25/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 = ',f20.10,'cm;')
20  format('draw ( ',f20.10,'u,',f20.10,'u)--(',f20.10,'u,',f20.10,'u) withcolor',a12,';')
30  format('fill ((',f20.10,'u,',f20.10,'u)--(',f20.10,'u,',f20.10,'u)--')
31  format('      (',f20.10,'u,',f20.10,'u)--(',f20.10,'u,',f20.10,'u)--')
32  format('      (',f20.10,'u,',f20.10,'u)--cycle) withcolor ',a30,';')              
40  format('drawdot(',f20.10,'u,',f20.10,'u);') 
50  format('endfig;',/,'end;')

end subroutine drawsolJ2

end module optvormod2
