module optvor

  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,J0delta,J2c(2),J2r,J3ctol,J4ctol
     integer, allocatable :: nvpols(:),poledges(:,:)
     real(kind=8), allocatable :: xpol(:,:),ypol(:,:)
     real(kind=8) :: Jweight(0:5)
     ! 0: minimum distance between sites
     ! 1: cells with identical volume given by vol(A)/m
     ! 2: cells with volume proportional to density function f evaluated at their sites
     ! 3: each cell with equally sized edges
     ! 4: each cell with equally sized angles
     ! 5: each cell with edges being crossed at their middle point by their corresponding Delaunay edges
  end type pdata_type
  
  public :: evalfg,evalfc,evalgjac,compJ,proj,approxg,drawsol,InA,compscale

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 comparea(scale,area,pdataptr)
	
    implicit none
  
    ! SCALAR ARGUMENTS
    real(kind=8), intent(in) :: scale
    real(kind=8), intent(out) :: area
    type(c_ptr), optional, intent(in) :: pdataptr
	  
    ! LOCAL SCALARS
    type(pdata_type), pointer :: pdata
    
    ! LOCAL SCALARS
    integer :: allocerr,j,maxnv
    real(kind=8) :: areaPj
    
    ! LOCAL ARRAYS
    real(kind=8), allocatable :: poltmp(:,:)
    
    call c_f_pointer(pdataptr,pdata)

    maxnv = maxval( pdata%nvpols(1:pdata%npols) )
  
    allocate(poltmp(2,maxnv),stat=allocerr)
    if ( allocerr .ne. 0 ) then
       write(*,*) 'Allocation error.'
       stop
    end if
        
    area = 0.0d0
    
    do j = 1,pdata%npols
       poltmp(1,1:pdata%nvpols(j)) = scale * pdata%xpol(j,1:pdata%nvpols(j))
       poltmp(2,1:pdata%nvpols(j)) = scale * pdata%ypol(j,1:pdata%nvpols(j))
     
       call polygon_area_2d(pdata%nvpols(j)-1,poltmp(1:2,1:pdata%nvpols(j)-1),areaPj)
        
       area = area + areaPj
    end do
	
    deallocate(poltmp,stat=allocerr)
    if ( allocerr .ne. 0 ) then
       write(*,*) 'Deallocation error.'
       stop
    end if

   end subroutine comparea
  
  ! ******************************************************************
  ! ******************************************************************

  subroutine compscale(nsites,scale,area,pdataptr)

    implicit none
  
    ! SCALAR ARGUMENTS
    integer, intent(in) :: nsites
    real(kind=8), intent(out) :: area,scale
    type(c_ptr), optional, intent(in) :: pdataptr

    ! PARAMETERS
    real(kind=8), parameter :: eps = 0.01d0
    
    ! LOCAL SCALARS
    real(kind=8) :: a,b,phia,phib,tphi,tscale,tarea
    type(pdata_type), pointer :: pdata
    
    call c_f_pointer(pdataptr,pdata)

    tscale = 1.0d0

    call comparea(tscale,tarea,pdataptr)

    tphi = tarea / nsites - 1.0d0

    a = tscale
    phia = tphi
    
    write(*,*) 'scale factor (a) = ',tscale,' ideal area (a) = ',tarea / nsites
    
    if ( abs( tphi ) .le. eps ) then
       scale = tscale
       area = tarea
       return
    end if

    if ( tphi .lt. 0.0d0 ) then
       do while ( tphi .lt. 0.0d0 )
          tscale = 2.0d0 * tscale
          call comparea(tscale,tarea,pdataptr)
          tphi = tarea / nsites - 1.0d0
       end do
       
       b = tscale
       phib = tphi
       
    else
       do while ( tphi .gt. 0.0d0 )
          tscale = 0.5d0 * tscale
          call comparea(tscale,tarea,pdataptr)
          tphi = tarea / nsites - 1.0d0
       end do
       
       b = tscale
       phib = tphi
    end if

    write(*,*) 'scale factor (b) = ',tscale,' ideal area (b) = ',tarea / nsites
    
    if ( abs( tphi ) .le. eps ) then
       scale = tscale
       area = tarea
       return
    end if

10  continue
    
    tscale = 0.5d0 * ( a + b )
    
    call comparea(tscale,tarea,pdataptr)

    tphi = tarea / nsites - 1.0d0

    write(*,*) 'scale factor (c) = ',tscale,' ideal area (c) = ',tarea / nsites
    
    if ( abs( tphi ) .le. eps ) then
       scale = tscale
       area = tarea
       return
    end if

    if ( sign(1.0d0,tphi) .eq. sign(1.0d0,phia) ) then
       a = tscale
       phia = tphi
    else
       b = tscale
       phib = tphi
    end if

    go to 10
    
  end subroutine compscale
	  
  ! ******************************************************************
  ! ******************************************************************

  subroutine proj(n,x,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(inout) :: x(n)

    ! LOCAL SCALARS
    integer :: i,nsites
    real(kind=8) :: tmp
    type(pdata_type), pointer :: pdata

    call c_f_pointer(pdataptr,pdata)

    flag = 0

    nsites = n / 2

    do i = 1,nsites
       tmp = norm2( x(2*i-1:2*i) - pdata%c(1:2) )
       if ( tmp .gt. pdata%rmax ) then
          x(2*i-1:2*i) = pdata%c(1:2) + pdata%rmax * ( x(2*i-1:2*i) - pdata%c(1:2) ) / tmp
       end if
    end do
    
  end subroutine proj
    
  ! ******************************************************************
  ! ******************************************************************

  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
    logical :: wempty
    integer :: allocerr,gtmpnnz,i,ierror,j,k,maxnv,nsites,ntriag,status
    real(kind=8) :: ftmp
    type(pdata_type), pointer :: pdata
    type(VoronoiCell) :: V
    type(Polygon) :: P
    
    ! 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))
    type(Polygon), allocatable :: W(:)
		
    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(2) .gt. 0.0d0 ) then
       do i = 1,nsites
          call fgcompJ2(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(W(pdata%npols),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

    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

       do j = 1,pdata%npols
          P%n = pdata%nvpols(j)
          P%deg = .false.
         
          P%v(1,1:P%n) = pdata%xpol(j,1:P%n)
          P%v(2,1:P%n) = pdata%ypol(j,1:P%n)
          P%e(1:P%n-1) = pdata%poledges(j,1:P%n-1)
       
          if ( nsites .ge. 3 .and. ierror .eq. 0 ) then
             call VorCellInterConvPol(sites(1:2,i),P,V,W(j))
          else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then
             call VorCellInterConvPol_Collinear(i,nsites,sites,P,W(j))
          end if
       end do

       if ( pdata%Jweight(0) .gt. 0.0d0 ) then
          call compJ0(nsites,sites,i,pdata%npols,W,pdata%J0delta,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad)
          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)
       end if

       if ( pdata%Jweight(1) .gt. 0.0d0 ) then
          call compJ1(nsites,sites,i,pdata%npols,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,pdata%npols,W,pdata%volA,fden,gden,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug)
          if ( flag .ne. 0 ) return

          f = f + pdata%Jweight(2) * ftmp ** 2
          if ( compgrad ) g(gtmpind(1:gtmpnnz)) = g(gtmpind(1:gtmpnnz)) + pdata%Jweight(2) * 2.0d0 * ftmp * gtmpval(1:gtmpnnz)
       end if
          
       wempty = .true.
       do j = 1,pdata%npols
          if ( .not. W(j)%deg ) then
             wempty = .false.
          end if
       end do
       
       if ( wempty ) then
          f = f + 1.0d0

       else
          do j = 1,pdata%npols
             if ( .not. W(j)%deg ) then
                if ( debug ) then
                   write(*,*)
                   write(*,*) 'POLYGON Wij with i = ',i,' and j = ',j
                   write(*,*) 'Number of vertices: ',W(j)%n - 1
                   write(*,*)
                end if
                   
                if ( pdata%Jweight(3) .gt. 0.0d0 ) then
                   call compJ3(nsites,sites,i,W(j),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
                      
                if ( pdata%Jweight(4) .gt. 0.0d0 ) then
                   call compJ4(nsites,sites,i,W(j),pdata%J4ctol,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug)
                   if ( flag .ne. 0 ) return

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

                   f = f + pdata%Jweight(5) * ftmp
                   if ( compgrad ) g(gtmpind(1:gtmpnnz)) = g(gtmpind(1:gtmpnnz)) + pdata%Jweight(5) * gtmpval(1:gtmpnnz)
                end if
                
             else
                if ( debug ) then
                   write(*,*)
                   write(*,*) 'POLYGON Wij with i = ',i,' and j = ',j,' IS EMPTY.'
                   write(*,*)
                end if
             end if
          end do
       end if
       
       call voronoi_cell_destroy(V)

       do j = 1,pdata%npols
          call polygon_destroy(W(j))
       end do
    end do

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

  end subroutine evalfg
  
  ! ******************************************************************
  ! ******************************************************************

  subroutine evalfc(n,x,f,m,c,flag,debug,pdataptr)
	
    implicit none
  
    ! SCALAR ARGUMENTS
    logical, intent(in) :: debug
    integer, intent(in) :: m,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) :: c(m)

    ! PARAMETER
    integer, parameter :: lim = 1
    logical, parameter :: compgrad = .false.
    
    ! LOCAL SCALARS
    logical :: lmem
    integer :: jcnnz
    
    ! LOCAL ARRAYS
    integer :: jcfun(lim),jcvar(lim)
    real(kind=8) :: g(1),jcval(lim)

    call evalfcgjac(n,x,f,g,m,c,jcfun,jcvar,jcval,jcnnz,lim,lmem,flag,.false.,debug,pdataptr)

  end subroutine evalfc
  
  ! ******************************************************************
  ! ******************************************************************

  subroutine evalgjac(n,x,g,m,jcfun,jcvar,jcval,jcnnz,lim,lmem,flag,debug,pdataptr)
	
    implicit none
  
    ! SCALAR ARGUMENTS
    logical, intent(in) :: debug
    logical, intent(out) :: lmem
    integer, intent(in) :: lim,m,n
    integer, intent(out) :: flag,jcnnz
    type(c_ptr), optional, intent(in) :: pdataptr
	  
    ! ARRAY ARGUMENTS
    integer, intent(out) :: jcfun(lim),jcvar(lim)
    real(kind=8), intent(in) :: x(n)
    real(kind=8), intent(out) :: g(n),jcval(lim)

    ! LOCAL SCALARS
    real(kind=8) :: f
    
    ! LOCAL ARRAYS
    real(kind=8) :: c(m)
		
    call evalfcgjac(n,x,f,g,m,c,jcfun,jcvar,jcval,jcnnz,lim,lmem,flag,.true.,debug,pdataptr)

  end subroutine evalgjac
  
  ! ******************************************************************
  ! ******************************************************************

  subroutine evalfcgjac(n,x,f,g,m,c,jcfun,jcvar,jcval,jcnnz,lim,lmem,flag,compgrad,debug,pdataptr)
	
    implicit none
  
    ! SCALAR ARGUMENTS
    logical, intent(in) :: compgrad,debug
    logical, intent(out) :: lmem
    integer, intent(in) :: lim,m,n
    integer, intent(out) :: flag,jcnnz
    real(kind=8), intent(out) :: f
    type(c_ptr), optional, intent(in) :: pdataptr
	  
    ! ARRAY ARGUMENTS
    integer, intent(out) :: jcfun(lim),jcvar(lim)
    real(kind=8), intent(in) :: x(n)
    real(kind=8), intent(out) :: c(m),g(n),jcval(lim)

    ! LOCAL SCALARS
    logical :: wempty
    integer :: allocerr,gtmpnnz,i,ierror,j,k,maxnv,nsites,ntriag,status
    real(kind=8) :: ftmp
    type(pdata_type), pointer :: pdata
    type(VoronoiCell) :: V
    type(Polygon) :: P
    
    ! LOCAL ARRAYS
    integer :: gtmpind(max(n,lim)),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(max(n,lim)),sites(2,n/2),vorxy(2,2*(n/2))
    type(Polygon), allocatable :: W(:)
		
    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 evalfcgjac: 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 evalfcgjac: 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 evalfcgjac: triangle_circumcenter_2d returned NaN.'
          flag = - 1
          return
       end if
    end do

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

    if ( compgrad ) then
       lmem = .false.
       g(1:n) = 0.0d0
       jcnnz = 0
    end if
    
    maxnv = maxval( pdata%nvpols(1:pdata%npols) )
      
    allocate(W(pdata%npols),P%v(2,maxnv),P%e(maxnv-1),stat=allocerr)
    if ( allocerr .ne. 0 ) then
       write(*,*) 'In evalfcgjac: Allocation error.'
       flag = - 1
       return
    end if

    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)
     ! call voronoi_cell_old(i,nsites,sites,ntriag,til,tnbr,vorxy,V,status)
       if ( status .ne. 0 ) then
          write(*,*) 'In evalfcgjac: Error in voronoi_cell, status = ',status
          flag = - 1
          return
       end if

       do j = 1,pdata%npols
          P%n = pdata%nvpols(j)
          P%deg = .false.
         
          P%v(1,1:P%n) = pdata%xpol(j,1:P%n)
          P%v(2,1:P%n) = pdata%ypol(j,1:P%n)
          P%e(1:P%n-1) = pdata%poledges(j,1:P%n-1)
       
          if ( nsites .ge. 3 .and. ierror .eq. 0 ) then
             call VorCellInterConvPol(sites(1:2,i),P,V,W(j))
          else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then
             call VorCellInterConvPol_Collinear(i,nsites,sites,P,W(j))
          end if
       end do

       call compJ0(nsites,sites,i,pdata%npols,W,pdata%J0delta,ftmp,max(n,lim),gtmpnnz,gtmpind,gtmpval,flag,compgrad)
       if ( flag .ne. 0 ) return

       f = f + 10.0d0 * ftmp
       if ( compgrad ) g(gtmpind(1:gtmpnnz)) = g(gtmpind(1:gtmpnnz)) + 10.0d0 * gtmpval(1:gtmpnnz)

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

       f = f + ftmp ** 2
       if ( compgrad ) g(gtmpind(1:gtmpnnz)) = g(gtmpind(1:gtmpnnz)) + 2.0d0 * ftmp * gtmpval(1:gtmpnnz)
                      
       wempty = .true.
       do j = 1,pdata%npols
          if ( .not. W(j)%deg ) then
             wempty = .false.
          end if
       end do
       
       if ( wempty ) then
          c(i) = 1.0d0
          c(nsites+i) = 1.0d0

       else
          do j = 1,pdata%npols
             if ( .not. W(j)%deg ) then
                call compJ3(nsites,sites,i,W(j),pdata%J3ctol,ftmp,max(n,lim),gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug)
                if ( flag .ne. 0 ) return

                !write(*,*) 'For i = ',i,' Constraint value = ',ftmp
                
                c(i) = ftmp

                if ( compgrad ) then
                   if ( jcnnz + gtmpnnz .gt. lim ) then
                      lmem = .true.
                      return
                   end if

                   if ( any(isnan(gtmpval(1:gtmpnnz))) ) then
                      write(*,*) 'There is a NaN here!'
                      call compJ3(nsites,sites,i,W(j),pdata%J3ctol,ftmp,max(n,lim),gtmpnnz,gtmpind,gtmpval,flag,compgrad,.true.)
                      if ( flag .ne. 0 ) return
                      stop
                   end if
                      
                   jcfun(jcnnz+1:jcnnz+gtmpnnz) = i
                   jcvar(jcnnz+1:jcnnz+gtmpnnz) = gtmpind(1:gtmpnnz)
                   jcval(jcnnz+1:jcnnz+gtmpnnz) = gtmpval(1:gtmpnnz)
                   jcnnz = jcnnz + gtmpnnz
                end if
                
                call compJ4(nsites,sites,i,W(j),pdata%J4ctol,ftmp,max(n,lim),gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug)
                if ( flag .ne. 0 ) return

                c(nsites+i) = ftmp

                if ( compgrad ) then
                   if ( jcnnz + gtmpnnz .gt. lim ) then
                      lmem = .true.
                      return
                   end if
                
                   if ( any(isnan(gtmpval(1:gtmpnnz))) ) then
                      write(*,*) 'There is a NaN here!'
                      call compJ4(nsites,sites,i,W(j),pdata%J4ctol,ftmp,max(n,lim),gtmpnnz,gtmpind,gtmpval,flag,compgrad,.true.)
                      if ( flag .ne. 0 ) return
                      stop
                   end if
                      
                   jcfun(jcnnz+1:jcnnz+gtmpnnz) = nsites + i
                   jcvar(jcnnz+1:jcnnz+gtmpnnz) = gtmpind(1:gtmpnnz)
                   jcval(jcnnz+1:jcnnz+gtmpnnz) = gtmpval(1:gtmpnnz)
                   jcnnz = jcnnz + gtmpnnz
                end if
             end if
          end do
       end if
       
       call voronoi_cell_destroy(V)

       do j = 1,pdata%npols
          call polygon_destroy(W(j))
       end do
    end do

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

  subroutine compJ(n,x,f,g,flag,compgrad,debug,pdataptr)
	
    implicit none
  
    ! SCALAR ARGUMENTS
    logical, intent(in) :: compgrad,debug
    integer, intent(in) :: n
    integer, intent(out) :: flag
    real(kind=8), intent(out) :: f(0:5)
    type(c_ptr), optional, intent(in) :: pdataptr
	  
    ! ARRAY ARGUMENTS
    real(kind=8), intent(in) :: x(n)
    real(kind=8), intent(out) :: g(n,0:5)

    ! PARAMETERS
    integer, parameter :: gtmpnnzmax = 100000
    
    ! LOCAL SCALARS
    logical :: wempty
    integer :: allocerr,gtmpnnz,i,ierror,j,k,maxnv,nsites,ntriag,status
    real(kind=8) :: ftmp
    type(pdata_type), pointer :: pdata
    type(VoronoiCell) :: V
    type(Polygon) :: P
    
    ! 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))
    type(Polygon), allocatable :: W(:)
		
    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 fgcompJ2(sites(1:2,i),fden(i),gden(1:2,i),compgrad,pdataptr)
    end do
    
    f(0:5) = 0.0d0

    if ( compgrad ) g(1:n,0:5) = 0.0d0
    
    maxnv = maxval( pdata%nvpols(1:pdata%npols) )
      
    allocate(W(pdata%npols),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

    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)
     ! call voronoi_cell_old(i,nsites,sites,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

       do j = 1,pdata%npols
          P%n = pdata%nvpols(j)
          P%deg = .false.
         
          P%v(1,1:P%n) = pdata%xpol(j,1:P%n)
          P%v(2,1:P%n) = pdata%ypol(j,1:P%n)
          P%e(1:P%n-1) = pdata%poledges(j,1:P%n-1)
       
          if ( nsites .ge. 3 .and. ierror .eq. 0 ) then
             call VorCellInterConvPol(sites(1:2,i),P,V,W(j))
          else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then
             call VorCellInterConvPol_Collinear(i,nsites,sites,P,W(j))
          end if
       end do

       call compJ0(nsites,sites,i,pdata%npols,W,pdata%J0delta,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad)
       if ( flag .ne. 0 ) return

       f(0) = f(0) + ftmp
       if ( compgrad ) g(gtmpind(1:gtmpnnz),0) = g(gtmpind(1:gtmpnnz),0) + gtmpval(1:gtmpnnz)

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

       f(1) = f(1) + ftmp ** 2
       if ( compgrad ) g(gtmpind(1:gtmpnnz),1) = g(gtmpind(1:gtmpnnz),1) + 2.0d0 * ftmp * gtmpval(1:gtmpnnz)
                      
       call compJ2(nsites,sites,i,pdata%npols,W,pdata%volA,fden,gden,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug)
       if ( flag .ne. 0 ) return

       f(2) = f(2) + ftmp ** 2
       if ( compgrad ) g(gtmpind(1:gtmpnnz),2) = g(gtmpind(1:gtmpnnz),2) + 2.0d0 * ftmp * gtmpval(1:gtmpnnz)
          
       wempty = .true.
       do j = 1,pdata%npols
          if ( .not. W(j)%deg ) then
             wempty = .false.
          end if
       end do
       
       if ( wempty ) then
          f = f + 1.0d0

       else
          do j = 1,pdata%npols
             if ( .not. W(j)%deg ) then
                call compJ3(nsites,sites,i,W(j),pdata%J3ctol,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug)
                if ( flag .ne. 0 ) return

                f(3) = f(3) + ftmp
                if ( compgrad ) g(gtmpind(1:gtmpnnz),3) = g(gtmpind(1:gtmpnnz),3) + gtmpval(1:gtmpnnz)
                      
                call compJ4(nsites,sites,i,W(j),pdata%J4ctol,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug)
                if ( flag .ne. 0 ) return

                f(4) = f(4) + ftmp
                if ( compgrad ) g(gtmpind(1:gtmpnnz),4) = g(gtmpind(1:gtmpnnz),4) + gtmpval(1:gtmpnnz)
                      
                call compJ5(nsites,sites,i,W(j),ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug)
                if ( flag .ne. 0 ) return

                f(5) = f(5) + ftmp
                if ( compgrad ) g(gtmpind(1:gtmpnnz),5) = g(gtmpind(1:gtmpnnz),5) + gtmpval(1:gtmpnnz)
             end if
          end do
       end if
       
       call voronoi_cell_destroy(V)

       do j = 1,pdata%npols
          call polygon_destroy(W(j))
       end do
    end do

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

  end subroutine compJ
  
  ! ******************************************************************
  ! ******************************************************************

  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:5,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
    logical :: wempty
    integer :: allocerr,gtmpnnz,i,ierror,j,k,maxnv,nsites,ntriag,status
    real(kind=8) :: ftmp
    type(pdata_type), pointer :: pdata
    type(VoronoiCell) :: V
    type(Polygon) :: P
    
    ! 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))
    type(Polygon), allocatable :: W(:)
		
    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 fgcompJ2(sites(1:2,i),fden(i),gden(1:2,i),compgrad,pdataptr)
    end do
    
    maxnv = maxval( pdata%nvpols(1:pdata%npols) )
      
    allocate(W(pdata%npols),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

    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)
     ! call voronoi_cell_old(i,nsites,sites,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

       do j = 1,pdata%npols
          P%n = pdata%nvpols(j)
          P%deg = .false.
         
          P%v(1,1:P%n) = pdata%xpol(j,1:P%n)
          P%v(2,1:P%n) = pdata%ypol(j,1:P%n)
          P%e(1:P%n-1) = pdata%poledges(j,1:P%n-1)
       
          if ( nsites .ge. 3 .and. ierror .eq. 0 ) then
             call VorCellInterConvPol(sites(1:2,i),P,V,W(j))
          else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then
             call VorCellInterConvPol_Collinear(i,nsites,sites,P,W(j))
          end if
       end do

       call compJ0(nsites,sites,i,pdata%npols,W,pdata%J0delta,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad)
       if ( flag .ne. 0 ) return

       f(0,i) = ftmp

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

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

       f(2,i) = ftmp
          
       wempty = .true.
       do j = 1,pdata%npols
          if ( .not. W(j)%deg ) then
             wempty = .false.
          end if
       end do
       
       if ( wempty ) then
          f(3,i) = 1.0d0
          f(4,i) = 1.0d0
          f(5,i) = 1.0d0

       else
          do j = 1,pdata%npols
             if ( .not. W(j)%deg ) then
                call compJ3(nsites,sites,i,W(j),pdata%J3ctol,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug)
                if ( flag .ne. 0 ) return

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

                f(4,i) = ftmp
                      
                call compJ5(nsites,sites,i,W(j),ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug)
                if ( flag .ne. 0 ) return

                f(5,i) = ftmp
             end if
          end do
       end if
       
       call voronoi_cell_destroy(V)

       do j = 1,pdata%npols
          call polygon_destroy(W(j))
       end do
    end do

    deallocate(W,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 compJ0(nsites,sites,i,npols,W,delta,f,gnnzmax,gnnz,gind,gval,flag,compgrad)
	
    implicit none
  
    ! SCALAR ARGUMENTS
    logical, intent(in) :: compgrad
    integer, intent(in) :: gnnzmax,i,npols,nsites
    integer, intent(out) :: gnnz
    integer, intent(inout) :: flag
    real(kind=8), intent(in) :: delta
    real(kind=8), intent(out) :: f
	  
    ! ARRAY ARGUMENTS
    integer, intent(out) :: gind(gnnzmax)
    real(kind=8), intent(in) :: sites(2,nsites)
    real(kind=8), intent(out) :: gval(gnnzmax)
    type(Polygon), intent(in) :: W(npols)

    ! LOCAL SCALARS
    integer :: ed,j
    real(kind=8) :: dist2
    
    ! LOCAL ARRAYS
    real(kind=8) :: diff(2)

    f = 0.0d0

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

       gnnz = 2
       gind(1:2) = (/ 2*i-1, 2*i /)
       gval(1:2) = 0.0d0
    end if

    do j = 1,npols
       if ( .not. W(j)%deg ) then
          do ed = 1,W(j)%n - 1
             if ( W(j)%e(ed) .lt. 0 ) then
                diff(1:2) = sites(1:2,i) - sites(1:2,abs(W(j)%e(ed)))
                dist2 = sum( diff(1:2) ** 2 )
                if ( dist2 .lt. delta ** 2 ) then
                   f = f + ( delta ** 2 - dist2 ) ** 2

                   if ( compgrad ) then
                      gval(1:2) = gval(1:2) - 4.0d0 * ( delta ** 2 - dist2 ) * diff(1:2)

                      if ( gnnz + 2 .gt. gnnzmax ) then
                         write(*,*) 'In compJ0: Increase gnnzmax to at least ',gnnz+2
                         flag = - 1
                         return
                      end if
        
                      gind(gnnz+1:gnnz+2) = (/ 2*abs(W(j)%e(ed))-1, 2*abs(W(j)%e(ed)) /)
                      gval(gnnz+1:gnnz+2) = 4.0d0 * ( delta ** 2 - dist2 ) * diff(1:2)
                      gnnz = gnnz + 2
                   end if
                end if
             end if
          end do
       end if
    end do
       
  end subroutine compJ0

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

  subroutine compJ1alternative(nsites,sites,i,npols,W,volA,f,gnnzmax,gnnz,gind,gval,flag,compgrad,debug)
	
    implicit none
  
    ! SCALAR ARGUMENTS
    logical, intent(in) :: compgrad,debug
    integer, intent(in) :: gnnzmax,i,npols,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(npols)
	  
    ! 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,j
    real(kind=8) :: mu,phia,phib,phic,vol,volj,xi
	    
    ! LOCAL ARRAYS
    real(kind=8) :: gradphi(2),M(2,2),normal(2)

    vol = 0.0d0
    do j = 1,npols
       if ( .not. W(j)%deg ) then
!!$          volj = 0.0d0
!!$          do ed = 1,W(j)%n - 1
!!$             volj = volj + 0.5d0 * ( W(j)%v(1,ed) * W(j)%v(2,ed+1) - W(j)%v(1,ed+1) * W(j)%v(2,ed) )
!!$          end do
          
          call polygon_area_2d(W(j)%n-1,W(j)%v(1:2,1:W(j)%n-1),volj)

          vol = vol + volj
       end if
    end do
    
    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
    
    gnnz = 0
       
    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*i-1, 2*i /)
    gval(gnnz+1:gnnz+2) = 0.0d0
    gnnz = gnnz + 2
    
    do j = 1,npols
       if ( .not. W(j)%deg ) then
          if ( debug ) then
             write(*,*)
             write(*,*) 'POLYGON Wij with i = ',i,' and j = ',j
             write(*,*) 'Number of vertices: ',W(j)%n - 1
             write(*,*)
          end if
          
          do ed = 1,W(j)%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(j)%v(1:2,ed)
                write(*,*) '===================='
                write(*,*)
             end if
             
             ! Processing edge (v_E,w_E)
             
             if ( W(j)%e(ed) .lt. 0 ) then
                
                if ( ed .eq. 1 ) then
                   edprev = W(j)%n - 1
                   ednext = ed + 1
                else if ( ed .eq. W(j)%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(j)%v(1:2,edprev),' and ',W(j)%v(1:2,ednext),' respectively.'
                end if
                
                ! Computing the outward unit normal vector to the edge ed
                call line_exp_normal_2d(W(j)%v(1:2,ednext),W(j)%v(1:2,ed),normal)
                   
                if ( debug ) then
                   write(*,*) 'Outward unit normal vector to the edge ed: ',normal(1:2)
                end if
                
                mu = 0.5d0 * norm2( W(j)%v(1:2,ed) - W(j)%v(1:2,ednext) )
                   
                if ( debug ) then
                   write(*,*) 'mu = ',mu
                end if
                
                ! Starting point v_E

                if ( debug ) then
                   write(*,*)
                   write(*,*) 'Processing v_E'
                   write(*,*)
                end if
          
                if ( W(j)%e(edprev) .lt. 0 ) then
                   if ( debug ) then
                      write(*,*) 'We are in the case |P(v_E)|=3. ', &
                           'Neighbour cells are: ',abs(W(j)%e(edprev)),' and ',abs(W(j)%e(ed))
                   end if
                
                   call matrixint(sites(1:2,abs(W(j)%e(edprev))),sites(1:2,abs(W(j)%e(ed))),sites(1:2,i),W(j)%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) ), normal(1:2) )
                   
                   call matrixint(sites(1:2,i),sites(1:2,abs(W(j)%e(edprev))),sites(1:2,abs(W(j)%e(ed))),W(j)%v(1:2,ed),M,flag)
                   if ( flag .ne. 0 ) return
                   
                   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(j)%e(ed))-1, 2*abs(W(j)%e(ed)) /)
                   gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) )
                   gnnz = gnnz + 2
                   
                   call matrixint(sites(1:2,abs(W(j)%e(ed))),sites(1:2,i),sites(1:2,abs(W(j)%e(edprev))),W(j)%v(1:2,ed),M,flag)
                   if ( flag .ne. 0 ) return

                   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(j)%e(edprev))-1, 2*abs(W(j)%e(edprev)) /)
                   gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) )
                   gnnz = gnnz + 2

                else
                   if ( debug ) then
                      write(*,*) 'We are in the case |P(v_E)|=2. Neighbour cell is: ',abs(W(j)%e(ed))
                   end if
                
                   call line_exp2imp_2d(W(j)%v(1:2,ed),W(j)%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(j)%e(ed))),sites(1:2,i),W(j)%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) ), normal(1:2) )
                   
                   call matrixbd(sites(1:2,i),sites(1:2,abs(W(j)%e(ed))),W(j)%v(1:2,ed),gradphi,M,flag)
                   if ( flag .ne. 0 ) return
                      
                   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(j)%e(ed))-1, 2*abs(W(j)%e(ed)) /)
                   gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) )
                   gnnz = gnnz + 2
                end if
                   
                ! Ending point w_E
                  
                if ( debug ) then
                   write(*,*)
                   write(*,*) 'Processing w_E'
                   write(*,*)
                end if
                
                if ( W(j)%e(ednext) .lt. 0 ) then
                   if ( debug ) then
                      write(*,*) 'We are in the case |P(w_E)|=3. ', &
                           'Neighbour cells are: ',abs(W(j)%e(ed)),' and ',abs(W(j)%e(ednext))
                   end if
                   
                   call matrixint(sites(1:2,abs(W(j)%e(ed))),sites(1:2,abs(W(j)%e(ednext))),sites(1:2,i),W(j)%v(1:2,ed+1),M,flag)
                   if ( flag .ne. 0 ) return
                   
                   gval(1:2) = gval(1:2) + xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) )

                   call matrixint(sites(1:2,i),sites(1:2,abs(W(j)%e(ed))),sites(1:2,abs(W(j)%e(ednext))),W(j)%v(1:2,ed+1),M,flag)
                   if ( flag .ne. 0 ) return
             
                   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(j)%e(ednext))-1, 2*abs(W(j)%e(ednext)) /)
                   gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) )
                   gnnz = gnnz + 2
             
                   call matrixint(sites(1:2,abs(W(j)%e(ednext))),sites(1:2,i),sites(1:2,abs(W(j)%e(ed))),W(j)%v(1:2,ed+1),M,flag)
                   if ( flag .ne. 0 ) return
             
                   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(j)%e(ed))-1, 2*abs(W(j)%e(ed)) /)
                   gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) )
                   gnnz = gnnz + 2
             
                else 
                   if ( debug ) then
                      write(*,*) 'We are in the case |P(w_E)|=2. Neighbour cell is: ',abs(W(j)%e(ed))
                   end if
                
                   call line_exp2imp_2d(W(j)%v(1:2,ednext),W(j)%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(j)%e(ed))),sites(1:2,i),W(j)%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) ), normal(1:2) )
                   
                   call matrixbd(sites(1:2,i),sites(1:2,abs(W(j)%e(ed))),W(j)%v(1:2,ednext),gradphi,M,flag)
                   if ( flag .ne. 0 ) return
                   
                   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(j)%e(ed))-1, 2*abs(W(j)%e(ed)) /)
                   gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) )
                   gnnz = gnnz + 2
                end if
                
             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 Wij with i = ',i,' and j = ',j,' IS EMPTY.'
             write(*,*)
          end if
       end if
    end do
    
  end subroutine compJ1alternative

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

  subroutine compJ1(nsites,sites,i,npols,W,volA,f,gnnzmax,gnnz,gind,gval,flag,compgrad,debug)
	
    implicit none
  
    ! SCALAR ARGUMENTS
    logical, intent(in) :: compgrad,debug
    integer, intent(in) :: gnnzmax,i,npols,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(npols)
	  
    ! 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,volj,xi
	    
    ! LOCAL ARRAYS
    real(kind=8) :: p(2)

    vol = 0.0d0
    do j = 1,npols
       if ( .not. W(j)%deg ) then
!!$          volj = 0.0d0
!!$          do ed = 1,W(j)%n - 1
!!$             volj = volj + 0.5d0 * ( W(j)%v(1,ed) * W(j)%v(2,ed+1) - W(j)%v(1,ed+1) * W(j)%v(2,ed) )
!!$          end do
          
          call polygon_area_2d(W(j)%n-1,W(j)%v(1:2,1:W(j)%n-1),volj)

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

                p(1:2) = 0.5d0 * ( W(j)%v(1:2,ed) + W(j)%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(j)%e(ed))-1, 2*abs(W(j)%e(ed)) /)
                gval(gnnz+1:gnnz+2) = - xi * mu * ( p(1:2) - sites(1:2,abs(W(j)%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 Wij with i = ',i,' and j = ',j,' IS EMPTY.'
             write(*,*)
          end if
       end if
    end do
    
  end subroutine compJ1

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

  subroutine compJ2(nsites,sites,i,npols,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,npols,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(npols)
	  
    ! 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,volj,xi
	    
    ! LOCAL ARRAYS
    real(kind=8) :: p(2)

    vol = 0.0d0
    do j = 1,npols
       if ( .not. W(j)%deg ) then
!!$          volj = 0.0d0
!!$          do ed = 1,W(j)%n - 1
!!$             volj = volj + 0.5d0 * ( W(j)%v(1,ed) * W(j)%v(2,ed+1) - W(j)%v(1,ed+1) * W(j)%v(2,ed) )
!!$          end do
          
          call polygon_area_2d(W(j)%n-1,W(j)%v(1:2,1:W(j)%n-1),volj)

          vol = vol + volj
       end if
    end do
    
    xi = nsites / volA

    if ( debug ) then
       write(*,*)
       write(*,*) 'xi = ',xi
       write(*,*)
    end if
       
  ! f = xi * vol - 1.0d0
    f = xi * vol - fden(i)

    if ( .not. compgrad ) then
       return
    end if
    
    if ( 2 .gt. gnnzmax ) then
       write(*,*) 'In compJ2: Increase gnnzmax to at least ',2
       flag = - 1
       return
    end if
    
    gind(1:2) = (/ 2*i-1, 2*i /)
  ! gval(1:2) = 0.0d0
    gval(1:2) = - gden(1:2,i)
    gnnz = 2
    
    do j = 1,npols
       if ( .not. W(j)%deg ) then
          if ( debug ) then
             write(*,*)
             write(*,*) 'POLYGON Wij with i = ',i,' and j = ',j
             write(*,*) 'Number of vertices: ',W(j)%n - 1
             write(*,*)
          end if
          
          do ed = 1,W(j)%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(j)%v(1:2,ed)
                write(*,*) 'Its outer vertex, named ',ednext,' as well, corresponds to the point ',W(j)%v(1:2,ednext)
                write(*,*) '===================='
                write(*,*)
             end if
             
             if ( W(j)%e(ed) .lt. 0 ) then
                
                if ( ed .eq. W(j)%n - 1 ) then
                   ednext = 1
                else
                   ednext = ed + 1
                end if
                
                if ( debug ) then
                   write(*,*) 'Next: ',ednext
                   write(*,*) 'It corresponds to the point ',W(j)%v(1:2,ednext)
                end if
                
                mu = norm2( W(j)%v(1:2,ed) - W(j)%v(1:2,ednext) ) / norm2( sites(1:2,i) - sites(1:2,abs(W(j)%e(ed))) )
                   
                if ( debug ) then
                   write(*,*) 'mu = ',mu
                end if

                p(1:2) = 0.5d0 * ( W(j)%v(1:2,ed) + W(j)%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 compJ2: Increase gnnzmax to at least ',gnnz+2
                   flag = - 1
                   return
                end if
        
                gind(gnnz+1:gnnz+2) = (/ 2*abs(W(j)%e(ed))-1, 2*abs(W(j)%e(ed)) /)
                gval(gnnz+1:gnnz+2) = - xi * mu * ( p(1:2) - sites(1:2,abs(W(j)%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 Wij with i = ',i,' and j = ',j,' IS EMPTY.'
             write(*,*)
          end if
       end if
    end do
    
  end subroutine compJ2

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

  subroutine compJ2alternative1(nsites,sites,i,npols,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,npols,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(npols)
	  
    ! 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,edprev,ednext,j
    real(kind=8) :: mu,phia,phib,phic,vol,volj,xi
	    
    ! LOCAL ARRAYS
    real(kind=8) :: gradphi(2),M(2,2),normal(2)

    vol = 0.0d0
    do j = 1,npols
       if ( .not. W(j)%deg ) then
          call polygon_area_2d(W(j)%n-1,W(j)%v(1:2,1:W(j)%n-1),volj)

          vol = vol + volj
       end if
    end do
    
    xi = nsites / volA

    if ( debug ) then
       write(*,*)
       write(*,*) 'xi = ',xi
       write(*,*)
    end if
       
  ! f = xi * vol - 1.0d0
    f = xi * vol - fden(i)

    if ( .not. compgrad ) then
       return
    end if
    
    gnnz = 0
       
    if ( gnnz + 2 .gt. gnnzmax ) then
       write(*,*) 'In compJ2: 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
    gval(gnnz+1:gnnz+2) = - gden(1:2,i)
    gnnz = gnnz + 2
    
    do j = 1,npols
       if ( .not. W(j)%deg ) then
          if ( debug ) then
             write(*,*)
             write(*,*) 'POLYGON Wij with i = ',i,' and j = ',j
             write(*,*) 'Number of vertices: ',W(j)%n - 1
             write(*,*)
          end if
          
          do ed = 1,W(j)%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(j)%v(1:2,ed)
                write(*,*) '===================='
                write(*,*)
             end if
             
             ! Processing edge (v_E,w_E)
             
             if ( W(j)%e(ed) .lt. 0 ) then
                
                if ( ed .eq. 1 ) then
                   edprev = W(j)%n - 1
                   ednext = ed + 1
                else if ( ed .eq. W(j)%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(j)%v(1:2,edprev),' and ',W(j)%v(1:2,ednext),' respectively.'
                end if
                
                ! Computing the outward unit normal vector to the edge ed
                call line_exp_normal_2d(W(j)%v(1:2,ednext),W(j)%v(1:2,ed),normal)
                   
                if ( debug ) then
                   write(*,*) 'Outward unit normal vector to the edge ed: ',normal(1:2)
                end if
                
                mu = 0.5d0 * norm2( W(j)%v(1:2,ed) - W(j)%v(1:2,ednext) )
                   
                if ( debug ) then
                   write(*,*) 'mu = ',mu
                end if
                
                ! Starting point v_E

                if ( debug ) then
                   write(*,*)
                   write(*,*) 'Processing v_E'
                   write(*,*)
                end if
          
                if ( W(j)%e(edprev) .lt. 0 ) then
                   if ( debug ) then
                      write(*,*) 'We are in the case |P(v_E)|=3. ', &
                           'Neighbour cells are: ',abs(W(j)%e(edprev)),' and ',abs(W(j)%e(ed))
                   end if
                
                   call matrixint(sites(1:2,abs(W(j)%e(edprev))),sites(1:2,abs(W(j)%e(ed))),sites(1:2,i),W(j)%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) ), normal(1:2) )
                   
                   call matrixint(sites(1:2,i),sites(1:2,abs(W(j)%e(edprev))),sites(1:2,abs(W(j)%e(ed))),W(j)%v(1:2,ed),M,flag)
                   if ( flag .ne. 0 ) return
                   
                   if ( gnnz + 2 .gt. gnnzmax ) then
                      write(*,*) 'In compJ2: Increase gnnzmax to at least ',gnnz+2
                      flag = - 1
                      return
                   end if
        
                   gind(gnnz+1:gnnz+2) = (/ 2*abs(W(j)%e(ed))-1, 2*abs(W(j)%e(ed)) /)
                   gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) )
                   gnnz = gnnz + 2
                   
                   call matrixint(sites(1:2,abs(W(j)%e(ed))),sites(1:2,i),sites(1:2,abs(W(j)%e(edprev))),W(j)%v(1:2,ed),M,flag)
                   if ( flag .ne. 0 ) return

                   if ( gnnz + 2 .gt. gnnzmax ) then
                      write(*,*) 'In compJ2: Increase gnnzmax to at least ',gnnz+2
                      flag = - 1
                      return
                   end if
        
                   gind(gnnz+1:gnnz+2) = (/ 2*abs(W(j)%e(edprev))-1, 2*abs(W(j)%e(edprev)) /)
                   gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) )
                   gnnz = gnnz + 2

                else
                   if ( debug ) then
                      write(*,*) 'We are in the case |P(v_E)|=2. Neighbour cell is: ',abs(W(j)%e(ed))
                   end if
                
                   call line_exp2imp_2d(W(j)%v(1:2,ed),W(j)%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(j)%e(ed))),sites(1:2,i),W(j)%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) ), normal(1:2) )
                   
                   call matrixbd(sites(1:2,i),sites(1:2,abs(W(j)%e(ed))),W(j)%v(1:2,ed),gradphi,M,flag)
                   if ( flag .ne. 0 ) return
                      
                   if ( gnnz + 2 .gt. gnnzmax ) then
                      write(*,*) 'In compJ2: Increase gnnzmax to at least ',gnnz+2
                      flag = - 1
                      return
                   end if
        
                   gind(gnnz+1:gnnz+2) = (/ 2*abs(W(j)%e(ed))-1, 2*abs(W(j)%e(ed)) /)
                   gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) )
                   gnnz = gnnz + 2
                end if
                   
                ! Ending point w_E
                  
                if ( debug ) then
                   write(*,*)
                   write(*,*) 'Processing w_E'
                   write(*,*)
                end if
                
                if ( W(j)%e(ednext) .lt. 0 ) then
                   if ( debug ) then
                      write(*,*) 'We are in the case |P(w_E)|=3. ', &
                           'Neighbour cells are: ',abs(W(j)%e(ed)),' and ',abs(W(j)%e(ednext))
                   end if
                   
                   call matrixint(sites(1:2,abs(W(j)%e(ed))),sites(1:2,abs(W(j)%e(ednext))),sites(1:2,i),W(j)%v(1:2,ed+1),M,flag)
                   if ( flag .ne. 0 ) return
                   
                   gval(1:2) = gval(1:2) + xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) )

                   call matrixint(sites(1:2,i),sites(1:2,abs(W(j)%e(ed))),sites(1:2,abs(W(j)%e(ednext))),W(j)%v(1:2,ed+1),M,flag)
                   if ( flag .ne. 0 ) return
             
                   if ( gnnz + 2 .gt. gnnzmax ) then
                      write(*,*) 'In compJ2: Increase gnnzmax to at least ',gnnz+2
                      flag = - 1
                      return
                   end if
        
                   gind(gnnz+1:gnnz+2) = (/ 2*abs(W(j)%e(ednext))-1, 2*abs(W(j)%e(ednext)) /)
                   gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) )
                   gnnz = gnnz + 2
             
                   call matrixint(sites(1:2,abs(W(j)%e(ednext))),sites(1:2,i),sites(1:2,abs(W(j)%e(ed))),W(j)%v(1:2,ed+1),M,flag)
                   if ( flag .ne. 0 ) return
             
                   if ( gnnz + 2 .gt. gnnzmax ) then
                      write(*,*) 'In compJ2: Increase gnnzmax to at least ',gnnz+2
                      flag = - 1
                      return
                   end if
        
                   gind(gnnz+1:gnnz+2) = (/ 2*abs(W(j)%e(ed))-1, 2*abs(W(j)%e(ed)) /)
                   gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) )
                   gnnz = gnnz + 2
             
                else 
                   if ( debug ) then
                      write(*,*) 'We are in the case |P(w_E)|=2. Neighbour cell is: ',abs(W(j)%e(ed))
                   end if
                
                   call line_exp2imp_2d(W(j)%v(1:2,ednext),W(j)%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(j)%e(ed))),sites(1:2,i),W(j)%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) ), normal(1:2) )
                   
                   call matrixbd(sites(1:2,i),sites(1:2,abs(W(j)%e(ed))),W(j)%v(1:2,ednext),gradphi,M,flag)
                   if ( flag .ne. 0 ) return
                   
                   if ( gnnz + 2 .gt. gnnzmax ) then
                      write(*,*) 'In compJ2: Increase gnnzmax to at least ',gnnz+2
                      flag = - 1
                      return
                   end if
        
                   gind(gnnz+1:gnnz+2) = (/ 2*abs(W(j)%e(ed))-1, 2*abs(W(j)%e(ed)) /)
                   gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) )
                   gnnz = gnnz + 2
                end if
                
             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 Wij with i = ',i,' and j = ',j,' IS EMPTY.'
             write(*,*)
          end if
       end if
    end do
    
  end subroutine compJ2alternative1
  
  ! ******************************************************************
  ! ******************************************************************

  subroutine compJ2alternative2(nsites,sites,i,npols,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,npols,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(npols)
	  
    ! 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,edprev,ednext,ind,j,k
    real(kind=8) :: mu,phia,phib,phic,sumfden,vol,volj,xi
	    
    ! LOCAL ARRAYS
    real(kind=8) :: gradphi(2),M(2,2),normal(2)

    vol = 0.0d0
    do j = 1,npols
       if ( .not. W(j)%deg ) then
          call polygon_area_2d(W(j)%n-1,W(j)%v(1:2,1:W(j)%n-1),volj)
          vol = vol + volj
       end if
    end do

    ! The code below corresponds to the formulation
    !
    ! ( vol - desired_vol )^2;
    !
    ! while the one that follows corresponds to the formulation
    !
    ! ( vol / desired_vol - 1 )^2.
    
!!$    sumfden = sum( fden(1:nsites) )
!!$    
!!$    f = f + ( vol - fden(i) * volA / sumfden ) ** 2
!!$
!!$    xi = 2.0d0 * ( vol - fden(i) * volA / sumfden )
!!$
!!$    ind = 2 * i - 1
!!$    g(ind:ind+1) = g(ind:ind+1) - ( xi * volA / sumfden ) * gden(1:2,i)
!!$
!!$    do k = 1,nsites
!!$       ind = 2 * k - 1
!!$       g(ind:ind+1) = g(ind:ind+1) + ( xi * volA * fden(i) / sumfden ** 2 ) * gden(1:2,k) 
!!$    end do

    sumfden = sum( fden(1:nsites) )

    xi = sumfden / ( fden(i) * 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 * nsites .gt. gnnzmax ) then
       write(*,*) 'In compJ2: Increase gnnzmax to at least ',2 * nsites
       flag = - 1
       return
    end if

    gnnz = 2 * nsites
    gind(1:gnnz) = (/ (ind,ind=1,gnnz) /)
    gval(1:gnnz) = 0.0d0
    
    ind = 2 * i - 1
    gval(ind:ind+1) = gval(ind:ind+1) - ( sumfden / ( fden(i) ** 2 * volA ) ) * vol * gden(1:2,i)

    do k = 1,nsites
       ind = 2 * k - 1
       gval(ind:ind+1) = gval(ind:ind+1) + ( 1.0d0 / ( fden(i) * volA ) ) * vol * gden(1:2,k) 
    end do

    do j = 1,npols
       if ( .not. W(j)%deg ) then
          if ( debug ) then
             write(*,*)
             write(*,*) 'POLYGON Wij with i = ',i,' and j = ',j
             write(*,*) 'Number of vertices: ',W(j)%n - 1
             write(*,*)
          end if
          
          do ed = 1,W(j)%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(j)%v(1:2,ed)
                write(*,*) '===================='
                write(*,*)
             end if
             
             ! Processing edge (v_E,w_E)
             
             if ( W(j)%e(ed) .lt. 0 ) then
                
                if ( ed .eq. 1 ) then
                   edprev = W(j)%n - 1
                   ednext = ed + 1
                else if ( ed .eq. W(j)%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(j)%v(1:2,edprev),' and ',W(j)%v(1:2,ednext),' respectively.'
                end if
                
                ! Computing the outward unit normal vector to the edge ed
                call line_exp_normal_2d(W(j)%v(1:2,ednext),W(j)%v(1:2,ed),normal)
                   
                if ( debug ) then
                   write(*,*) 'Outward unit normal vector to the edge ed: ',normal(1:2)
                end if
                
                mu = 0.5d0 * norm2( W(j)%v(1:2,ed) - W(j)%v(1:2,ednext) )
                   
                if ( debug ) then
                   write(*,*) 'mu = ',mu
                end if
                
                ! Starting point v_E

                if ( debug ) then
                   write(*,*)
                   write(*,*) 'Processing v_E'
                   write(*,*)
                end if
          
                if ( W(j)%e(edprev) .lt. 0 ) then
                   if ( debug ) then
                      write(*,*) 'We are in the case |P(v_E)|=3. ', &
                           'Neighbour cells are: ',abs(W(j)%e(edprev)),' and ',abs(W(j)%e(ed))
                   end if
                
                   call matrixint(sites(1:2,abs(W(j)%e(edprev))),sites(1:2,abs(W(j)%e(ed))),sites(1:2,i),W(j)%v(1:2,ed),M,flag)
                   if ( flag .ne. 0 ) return
                      
                   ind = 2 * i - 1
                   gval(ind:ind+1) = gval(ind:ind+1) + xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) )
                   
                   call matrixint(sites(1:2,i),sites(1:2,abs(W(j)%e(edprev))),sites(1:2,abs(W(j)%e(ed))),W(j)%v(1:2,ed),M,flag)
                   if ( flag .ne. 0 ) return
                   
                   ind = 2 * abs(W(j)%e(ed)) - 1
                   gval(ind:ind+1) = gval(ind:ind+1) + xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) )
                   
                   call matrixint(sites(1:2,abs(W(j)%e(ed))),sites(1:2,i),sites(1:2,abs(W(j)%e(edprev))),W(j)%v(1:2,ed),M,flag)
                   if ( flag .ne. 0 ) return

                   ind = 2 * abs(W(j)%e(edprev)) - 1
                   gval(ind:ind+1) = gval(ind:ind+1) + xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) )

                else
                   if ( debug ) then
                      write(*,*) 'We are in the case |P(v_E)|=2. Neighbour cell is: ',abs(W(j)%e(ed))
                   end if
                
                   call line_exp2imp_2d(W(j)%v(1:2,ed),W(j)%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(j)%e(ed))),sites(1:2,i),W(j)%v(1:2,ed),gradphi,M,flag)
                   if ( flag .ne. 0 ) return

                   ind = 2 * i - 1
                   gval(ind:ind+1) = gval(ind:ind+1) + xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) )
                   
                   call matrixbd(sites(1:2,i),sites(1:2,abs(W(j)%e(ed))),W(j)%v(1:2,ed),gradphi,M,flag)
                   if ( flag .ne. 0 ) return
                      
                   ind = 2 * abs(W(j)%e(ed)) - 1
                   gval(ind:ind+1) = gval(ind:ind+1) + xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) )
                end if
                   
                ! Ending point w_E
                  
                if ( debug ) then
                   write(*,*)
                   write(*,*) 'Processing w_E'
                   write(*,*)
                end if
                
                if ( W(j)%e(ednext) .lt. 0 ) then
                   if ( debug ) then
                      write(*,*) 'We are in the case |P(w_E)|=3. ', &
                           'Neighbour cells are: ',abs(W(j)%e(ed)),' and ',abs(W(j)%e(ednext))
                   end if
                   
                   call matrixint(sites(1:2,abs(W(j)%e(ed))),sites(1:2,abs(W(j)%e(ednext))),sites(1:2,i),W(j)%v(1:2,ed+1),M,flag)
                   if ( flag .ne. 0 ) return
                   
                   ind = 2 * i - 1
                   gval(ind:ind+1) = gval(ind:ind+1) + xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) )

                   call matrixint(sites(1:2,i),sites(1:2,abs(W(j)%e(ed))),sites(1:2,abs(W(j)%e(ednext))),W(j)%v(1:2,ed+1),M,flag)
                   if ( flag .ne. 0 ) return
             
                   ind = 2 * abs(W(j)%e(ednext)) - 1
                   gval(ind:ind+1) = gval(ind:ind+1) + xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) )
             
                   call matrixint(sites(1:2,abs(W(j)%e(ednext))),sites(1:2,i),sites(1:2,abs(W(j)%e(ed))),W(j)%v(1:2,ed+1),M,flag)
                   if ( flag .ne. 0 ) return
             
                   ind = 2 * abs(W(j)%e(ed)) - 1
                   gval(ind:ind+1) = gval(ind:ind+1) + xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) )
             
                else 
                   if ( debug ) then
                      write(*,*) 'We are in the case |P(w_E)|=2. Neighbour cell is: ',abs(W(j)%e(ed))
                   end if
                
                   call line_exp2imp_2d(W(j)%v(1:2,ednext),W(j)%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(j)%e(ed))),sites(1:2,i),W(j)%v(1:2,ednext),gradphi,M,flag)
                   if ( flag .ne. 0 ) return

                   ind = 2 * i - 1
                   gval(ind:ind+1) = gval(ind:ind+1) + xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) )
                   
                   call matrixbd(sites(1:2,i),sites(1:2,abs(W(j)%e(ed))),W(j)%v(1:2,ednext),gradphi,M,flag)
                   if ( flag .ne. 0 ) return
                   
                   ind = 2 * abs(W(j)%e(ed)) - 1
                   gval(ind:ind+1) = gval(ind:ind+1) + xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) )
                end if
             
             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 Wij with i = ',i,' and j = ',j,' IS EMPTY.'
             write(*,*)
          end if
       end if
    end do
    
  end subroutine compJ2alternative2

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

  subroutine fgcompJ2(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

    ! LOCAL ARRAY
    real(kind=8) :: z(2)
    
    call c_f_pointer(pdataptr,pdata)

    z(1:2) = (/ 2.0d0, 2.0d0 /) + 0.40d0 * x(1:2)
    if ( ( z(2) - ( 0.25d0 * z(1) ) ** 2 ) ** 2 + ( 0.25d0 * z(1) - 1.0d0 ) ** 2 - 1.0d0 .le. 0.0d0 ) then
       fden = 0.25d0
    else
       fden = 1.075d0
    end if
    
    if ( compgrad ) gden(1:2) = 0.0d0

!!$    fden = 2.5d0 - ( 2.5d0 - 0.5d0 ) * sum( ( x(1:2) - pdata%J2c(1:2) ) ** 2 ) / pdata%J2r ** 2
!!$    
!!$    if ( compgrad ) gden(1:2) = - ( 2.5d0 - 0.5d0 ) * 2.0d0 * ( x(1:2) - pdata%J2c(1:2) ) / pdata%J2r ** 2

  end subroutine fgcompJ2
    
  ! ******************************************************************
  ! ******************************************************************

  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,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 compJ3

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

  subroutine compJ4(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 compJ4

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

  subroutine compJ5(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,nmov
    real(kind=8) :: phia,phib,phic,sizeed,xi
	    
    ! LOCAL ARRAYS
    real(kind=8) :: gradphi(2),M(2,2),mu(2),p(2),q(2),diff(2)
		
    nmov = count( W%e(1:W%n-1) .lt. 0 )

    if ( debug ) then
       write(*,*)
       write(*,*) 'The polygon has nmov = ',nmov,' moving edges.'
       write(*,*)
    end if
    
    xi = 1.0d0 / nmov

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

    f = 0.0d0

    if ( compgrad ) then
       gnnz = 0
    
       if ( gnnz + 2 .gt. gnnzmax ) then
          write(*,*) 'In compJ5: 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 ( W%e(ed) .lt. 0 ) then
          p(1:2) = 0.5d0 * ( W%v(1:2,ed) + W%v(1:2,ednext) )
          
          q(1:2) = 0.5d0 * ( sites(1:2,i) + sites(1:2,abs(W%e(ed))) )
                
          diff(1:2) = p(1:2) - q(1:2)

          sizeed = norm2( W%v(1:2,ednext) - W%v(1:2,ed) )

          f = f + xi * sum( diff(1:2) ** 2 ) / sizeed ** 2

          if ( .not. compgrad ) then
             cycle
          end if
          
          gval(1:2) = gval(1:2) - xi * diff(1:2) / sizeed ** 2
                
          if ( gnnz + 2 .gt. gnnzmax ) then
             write(*,*) 'In compJ5: 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 * diff(1:2) / sizeed ** 2
          gnnz = gnnz + 2

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

          mu(1:2) = diff(1:2) / sizeed ** 2 + ( 2.0d0 * sum( diff(1:2) ** 2 ) / sizeed ** 4 ) * ( W%v(1:2,ednext) - W%v(1:2,ed) )
          
          if ( all( mu(1:2) .eq. 0.0d0 ) ) then
             cycle
          end if
       
          if ( 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) ), mu(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 compJ5: 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) ), mu(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 compJ5: 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) ), mu(1:2) )
             gnnz = gnnz + 2

          else 
             if ( debug ) then
                write(*,*) 'We are in the case |P(v_E)|=2. 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) ), mu(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 compJ5: 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) ), mu(1:2) )
             gnnz = gnnz + 2
          end if
                   
          ! Ending point w_E
                  
          if ( debug ) then
             write(*,*)
             write(*,*) 'Processing w_E'
             write(*,*)
          end if
                
          mu(1:2) = diff(1:2) / sizeed ** 2 - ( 2.0d0 * sum( diff(1:2) ** 2 ) / sizeed ** 4 ) * ( W%v(1:2,ednext) - W%v(1:2,ed) )
          
          if ( 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) ), mu(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 compJ5: 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) ), mu(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 compJ5: 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) ), mu(1:2) )
             gnnz = gnnz + 2
                   
          else 
             if ( debug ) then
                write(*,*) 'We are in the case |P(w_E)|=2. 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) ), mu(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 compJ5: 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) ), mu(1:2) )
             gnnz = gnnz + 2
          end if
       end if
    end do
       
  end subroutine compJ5

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

  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 drawsol1(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))

    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 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,';')
40  format('drawdot(',f20.10,'u,',f20.10,'u);') 
50  format('endfig;',/,'end;')
    
  end subroutine drawsol1

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

  subroutine drawsol2(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, parameter :: nstats = 1000
    integer :: flag,find
    logical :: debug
    real(kind=8) :: ctol,J3ctoltmp,J4ctoltmp
    integer :: cellcolorJ3(n/2)
    real(kind=8) :: f(0:5,n/2),stats(nstats,3)
    character(len=80) :: paintcolor,statsfn

    call c_f_pointer(pdataptr,pdata)

    J3ctoltmp = pdata%J3ctol
    J4ctoltmp = pdata%J4ctol
       
    do k = 1,nstats
       ctol = 1.0d0 - ( nstats - real( k ) ) / ( nstats - 1 )
       pdata%J3ctol = ctol
       pdata%J4ctol = 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

       do j = 1,3
          if ( j .eq. 1 ) then
             find = 1
             stats(k,j) = real( count( min( 0.0d0, f(find,1:n/2) + 1.0d0 - ctol ) ** 2 .le. 1.0d-04 ) ) / ( n / 2 )
          else if ( j .eq. 2 ) then
             find = 3
             stats(k,j) = real( count( f(find,1:n/2) .le. 1.0d-04 ) ) / ( n / 2 )
          else ! if ( j .eq. 3 ) then
             find = 4
             stats(k,j) = real( count( f(find,1:n/2) .le. 1.0d-04 ) ) / ( n / 2 )
          end if
       end do
    end do

    write(statsfn,"(A6,F3.1,A1,F3.1,A4)") 'stats-',J3ctoltmp,'-',J4ctoltmp,'.dat'
    open(unit=10,file=trim(statsfn))                           
    do k = 1,nstats
       write(10,1000) 1.0d0 - ( nstats - real( k ) ) / ( nstats - 1 ),stats(k,1:3)
    end do
    close(10)

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

    do k = 5,1,-1
       ctol = 0.1d0 * k
       pdata%J3ctol = 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(3,1:n/2) .gt. 1.0d-04 )
          cellcolorJ3(1:n/2) = k
       end where
    end do

    pdata%J3ctol = J3ctoltmp
    pdata%J4ctol = J4ctoltmp
       
    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 ( cellcolorJ3(i) .ne. - 1 ) then
          if ( cellcolorJ3(i) .eq. 5 ) then
             paintcolor = 'bluea'
          else if ( cellcolorJ3(i) .eq. 4 ) then
             paintcolor = 'blueb'
          else if ( cellcolorJ3(i) .eq. 3 ) then
             paintcolor = 'bluec'
          else if ( cellcolorJ3(i) .eq. 2 ) then
             paintcolor = 'blued'
          else if ( cellcolorJ3(i) .eq. 1 ) then
             paintcolor = 'bluee'
          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 adding a few statistics
    
    write(10,*) '% Area of A = ',pdata%volA
    write(10,*) '% Number of polygons that define A = ',pdata%npols
    write(10,*) '% Number of sites = ',nsites
    
    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;')
    
1000 format(4(1X,F10.4))
    
  end subroutine drawsol2

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

  subroutine drawsol3(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, parameter :: nstats = 1000
    integer :: flag,find
    logical :: debug
    real(kind=8) :: ctol,J3ctoltmp,J4ctoltmp
    integer :: cellcolorJ4(n/2)
    real(kind=8) :: f(0:5,n/2),stats(nstats,3)
    character(len=80) :: paintcolor,statsfn

    call c_f_pointer(pdataptr,pdata)

    J3ctoltmp = pdata%J3ctol
    J4ctoltmp = pdata%J4ctol
       
    do k = 1,nstats
       ctol = 1.0d0 - ( nstats - real( k ) ) / ( nstats - 1 )
       pdata%J3ctol = ctol
       pdata%J4ctol = 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

       do j = 1,3
          if ( j .eq. 1 ) then
             find = 1
             stats(k,j) = real( count( min( 0.0d0, f(find,1:n/2) + 1.0d0 - ctol ) ** 2 .le. 1.0d-04 ) ) / ( n / 2 )
          else if ( j .eq. 2 ) then
             find = 3
             stats(k,j) = real( count( f(find,1:n/2) .le. 1.0d-04 ) ) / ( n / 2 )
          else ! if ( j .eq. 3 ) then
             find = 4
             stats(k,j) = real( count( f(find,1:n/2) .le. 1.0d-04 ) ) / ( n / 2 )
          end if
       end do
    end do

    write(statsfn,"(A6,F3.1,A1,F3.1,A4)") 'stats-',J3ctoltmp,'-',J4ctoltmp,'.dat'
    open(unit=10,file=trim(statsfn))                           
    do k = 1,nstats
       write(10,1000) 1.0d0 - ( nstats - real( k ) ) / ( nstats - 1 ),stats(k,1:3)
    end do
    close(10)

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

    do k = 7,5,-1
       ctol = 0.1d0 * k
       pdata%J4ctol = 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(4,1:n/2) .gt. 1.0d-04 )
          cellcolorJ4(1:n/2) = k
       end where
    end do

    pdata%J3ctol = J3ctoltmp
    pdata%J4ctol = J4ctoltmp
       
    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)
     ! call voronoi_cell_old(i,nsites,sites,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 ( cellcolorJ4(i) .ne. - 1 ) then
          if ( cellcolorJ4(i) .eq. 7 ) then
             paintcolor = 'corn'
          else if ( cellcolorJ4(i) .eq. 6 ) then
             paintcolor = 'deepsaffron'
          else if ( cellcolorJ4(i) .eq. 5 ) then
             paintcolor = 'ferrarired'
          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 adding a few statistics
    
    write(10,*) '% Area of A = ',pdata%volA
    write(10,*) '% Number of polygons that define A = ',pdata%npols
    write(10,*) '% Number of sites = ',nsites
    
    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;')
    
1000 format(4(1X,F10.4))
    
  end subroutine drawsol3

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

  subroutine drawsol4(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:5,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
    
    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 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
                ! This is to draw the distance between the midpoint of
                ! Delaunay and Voronoy corresponding edges
                if ( Wij%e(v) .lt. 0 ) then
                   write(10,20) 0.5d0 * ( sites(1:2,i) + sites(1:2,abs(Wij%e(v))) ), &
                                0.5d0 * ( Wij%v(1:2,v) + Wij%v(1:2,v+1) ),'royalblue'
                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 J0i (sum/ave/min/max) = ',sum(abs(f(0,1:n/2))),sum(abs(f(0,1:n/2)))/(n/2), &
         minval(abs(f(0,1:n/2))),maxval(abs(f(0,1:n/2)))
    write(10,*) 'Value of J1i (sum/ave/min/max) = ',sum(abs(f(1,1:n/2))),sum(abs(f(1,1:n/2)))/(n/2), &
         minval(abs(f(1,1:n/2))),maxval(abs(f(1,1:n/2))),' J1 itself = ',sum(f(1,1:n/2)**2)/(n/2)
    write(10,*) 'Value of J5i (sum/ave/min/max) = ',sum(abs(f(5,1:n/2))),sum(abs(f(5,1:n/2)))/(n/2), &
         minval(abs(f(5,1:n/2))),maxval(abs(f(5,1: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 drawsol4

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

  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:5,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
    
    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 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 J2i (sum/ave/min/max) = ',sum(abs(f(2,1:n/2))),sum(abs(f(2,1:n/2)))/(n/2), &
         minval(abs(f(2,1:n/2))),maxval(abs(f(2,1:n/2))),' J2 itself = ',sum(f(2,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

end module optvor
