program density_profile_modulated_main_2d
use density_profile_modulated_def_2d
implicit none

integer i

	print*, 'Attention: the density profile must be an even function!'
	call define_grid ()				! define x and q vectors
	call define_potential ()		! define potential on the grid points
	call define_weight_planar_2d ()	! define the characteristic function and their FT in geometry where only 
	
	do i = 0,niterations
		call main_iteration (i)
		print*, '------------------------------------------------------'
		
	enddo


contains
	
	subroutine compute_weighted_densities_2d ()
	implicit none
	! weighted densities
	!  wd (1,:) = n2 == n[2,2]*R
	!  wd (2,:) = n(0) = n0 * 2*pi*R
	!  wd (3,:) = n(1)
	!  wd (4,:) = n[1,1]
	complex*16 y(0:n-1),ffty(0:n-1)
	integer i
	
	do i=0,n-1
		ffty(i) = dcmplx(rhoqx(i)*wqx(1,i),0.d0)		! real and even
	enddo
	call fft_mik1 (n,x,y,qx,ffty,-1,1,1)
	wd(1,:) = dreal(y(:))								! real and even
	
	do i=0,n-1
		ffty(i) = dcmplx(rhoqx(i)*wqx(2,i),0.d0)		! real and even
	enddo
	call fft_mik1 (n,x,y,qx,ffty,-1,1,1)
	wd(2,:) = dreal(y(:))								! real and even
	
	do i=0,n-1
		ffty(i) = dcmplx(0.d0,rhoqx(i)*wqx(3,i))		! immaginary and odd
	enddo
	call fft_mik1 (n,x,y,qx,ffty,-1,0,1)
	wd(3,:) = dreal(y(:))								! real and odd
	
	do i=0,n-1
		ffty(i) = dcmplx(rhoqx(i)*wqx(4,i),0.d0)		! real and even
	enddo
	call fft_mik1 (n,x,y,qx,ffty,-1,1,1)
	wd(4,:) = dreal(y(:))								! real and even
	
!~ 	write(*,'(4(F8.5,1X,F8.5,3X))')  wd(1,0), eta, wd (2,0), 2*eta/R, wd (3,0),0.0, wd (4,0), eta/R
!~ 	write(*,'(4(F8.5,1X,F8.5,3X))')  wd(1,10), eta, wd (2,10), 2*eta/R, wd (3,10),0.0, wd (4,10), eta/R
!~ 	write(*,'(4(F8.5,1X,F8.5,3X))')  wd(1,20), eta, wd (2,20), 2*eta/R, wd (3,20),0.0, wd (4,20), eta/R
!~ 	write(*,'(4(F8.5,1X,F8.5,3X))')  wd(1,30), eta, wd (2,30), 2*eta/R, wd (3,30),0.0, wd (4,30), eta/R
!~ 	stop
	end subroutine compute_weighted_densities_2d

	subroutine compute_pd_p_2d ()
	implicit none
	! weighted densities
	!  wd (1,:) = n2 == n[2,2]*R
	!  wd (2,:) = n(0) = n0 * 2*pi*R
	!  wd (3,:) = n(1)
	!  wd (4,:) = n[1,1]
	integer i1
	
	select case (ndft)
		case(1)
			! According to Eq.(9) of Roth et al. J. Chem. Phys. 136, 081101 (2012)
			pd(0,:) = - dlog( 1.d0 - wd(1,:) )
			pd(1,:) = wd(2,:)/(2*pi*R) * 1.d0/(1.d0 - wd(1,:)) + &
			        1.d0/(4*pi*(1-wd(1,:))**2) * &
					(19.d0/12 *wd(2,:)**2 - 5.d0/12 *wd(3,:)**2  - 7.d0/6 *(wd(1,:)/R)**2  - 7.d0/6 *wd(4,:)**2)
			pd(2,:) = 1.d0/(4*pi*(1-wd(1,:))) * 19.d0/6 *wd(2,:)
			pd(3,:) = -1.d0/(4*pi*(1-wd(1,:))) * 5.d0/6 *wd(3,:)
			pd(4,:) = -1.d0/(4*pi*(1-wd(1,:))) * 7.d0/3 *wd(1,:)/R
			pd(5,:) = -1.d0/(4*pi*(1-wd(1,:))) * 7.d0/3 *wd(4,:)
			! Rosenfeld function itself
			p (:) = -wd(2,:)/(2*pi*R) * dlog( 1.d0 - wd(1,:) ) + &
					1.d0/(4*pi*(1-wd(1,:))) * &
					(19.d0/12 *wd(2,:)**2 - 5.d0/12 *wd(3,:)**2  - 7.d0/6 *(wd(1,:)/R)**2  - 7.d0/6 *wd(4,:)**2)
		case(2)
			stop 'not ready yet'			
	end select
!~ 	write(*,'(6(F8.5,1X,F8.5,3X))')  pd(0,0),-dlog(1.d0-eta), &
!~ 	                                 pd(1,0), rho_g/(1.d0-eta) + eta**2/(pi*R**2*(1.d0-eta)**2), &
!~ 	                                 pd(2,0), 19*eta/(12*pi*R*(1.d0-eta)), &
!~ 	                                 pd(3,0),0.d0, &
!~ 	                                 pd(4,0), -7*eta/(12*pi*R*(1.d0-eta))  , &
!~ 	                                 pd(5,0),-7*eta/(12*pi*R*(1.d0-eta))
!~ 	write(*,'(6(F8.5,1X,F8.5,3X))')  pd(0,10),-dlog(1.d0-eta), &
!~ 	                                 pd(1,10), rho_g/(1-eta) + eta**2/(pi*R**2*(1.d0-eta)**2), &
!~ 	                                 pd(2,10), 19*eta/(12*pi*R*(1.d0-eta)), &
!~ 	                                 pd(3,10),0.d0, &
!~ 	                                 pd(4,10),-7*eta/(12*pi*R*(1.d0-eta))  , &
!~ 	                                 pd(5,10),-7*eta/(12*pi*R*(1.d0-eta))
!~ 	write(*,'(6(F8.5,1X,F8.5,3X))')  pd(0,20),-dlog(1.d0-eta), &
!~ 	                                 pd(1,20), rho_g/(1-eta) + eta**2/(pi*R**2*(1.d0-eta)**2), &
!~ 	                                 pd(2,20), 19*eta/(12*pi*R*(1.d0-eta)), &
!~ 	                                 pd(3,20),0.d0, &
!~ 	                                 pd(4,20),-7*eta/(12*pi*R*(1.d0-eta))  , &
!~ 	                                 pd(5,20),-7*eta/(12*pi*R*(1.d0-eta))
!~ 	stop
	end subroutine compute_pd_p_2d

	subroutine compute_pdw_2d ()
	implicit none
	! Define weight functions
	!    wx(1,:) == w2 == w[2,2]*R
	!    wx(2,:) == w(0) == w0*2*pi*R
	!    wx(3,:) == w(1)
	!	 wx(4,:) == w[1,1]
	complex*16 y(0:n-1),ffty(0:n-1)
	integer i,i1,i2
	complex*16 c1, c2, c3, c4
	integer j1, j2
	
	! d Phi/ d n0
	do i=0,n-1
		y(i) = dcmplx(pd(0,i),0.d0)								! real and even
	enddo
	call fft_mik1 (n,x,y,qx,ffty,1,1,1)
	do i=0,n-1
		ffty(i) = dcmplx(dreal(ffty(i))*wqx(2,i)/(2*pi*R),0.d0)		! real and even
	enddo
	call fft_mik1 (n,x,y,qx,ffty,-1,1,1)
	pdw(0,:) = dreal(y(:))										! real and even
	
	! d Phi/ d n2
	do i=0,n-1
		y(i) = dcmplx(pd(1,i),0.d0)								! real and even
	enddo
	call fft_mik1 (n,x,y,qx,ffty,1,1,1)
	do i=0,n-1
		ffty(i) = dcmplx(dreal(ffty(i))*wqx(1,i),0.d0)					! real and even
	enddo
	call fft_mik1 (n,x,y,qx,ffty,-1,1,1)
	pdw(1,:) = dreal(y(:))										! real and even
	
	! d Phi/ d n(0)
	do i=0,n-1
		y(i) = dcmplx(pd(2,i),0.d0)								! real and even
	enddo
	call fft_mik1 (n,x,y,qx,ffty,1,1,1)
	do i=0,n-1
		ffty(i) = dcmplx(dreal(ffty(i))*wqx(2,i),0.d0)			! real and even
	enddo
	call fft_mik1 (n,x,y,qx,ffty,-1,1,1)
	pdw(2,:) = dreal(y(:))										! real and even
	
	! d Phi/ d n(1)   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! riflettere su questo pezzo
	do i=0,n-1
		y(i) = dcmplx(pd(3,i),0.d0)								! real and odd
	enddo
	call fft_mik1 (n,x,y,qx,ffty,1,0,1)
	do i=0,n-1
		ffty(i) = dcmplx(dimag(ffty(i))*wqx(3,i),0.d0)			! real and even   un segno meno perche` e` prodotto di immaginari ed un segno meno perche` devo invertire la convoluzione
	enddo
	call fft_mik1 (n,x,y,qx,ffty,-1,1,1)
	pdw(3,:) = dreal(y(:))										! real and even
	
	! d Phi/ d n[2,2]
	do i=0,n-1
		y(i) = dcmplx(pd(4,i),0.d0)								! real and even
	enddo
	call fft_mik1 (n,x,y,qx,ffty,1,1,1)
	do i=0,n-1
		ffty(i) = dcmplx(dreal(ffty(i))*wqx(1,i)/R,0.d0)		! real and even
	enddo
	call fft_mik1 (n,x,y,qx,ffty,-1,1,1)
	pdw(4,:) = dreal(y(:))										! real and even
	
	! d Phi/ d n[1,1]
	do i=0,n-1
		y(i) = dcmplx(pd(5,i),0.d0)								! real and even
	enddo
	call fft_mik1 (n,x,y,qx,ffty,1,1,1)
	do i=0,n-1
		ffty(i) = dcmplx(dreal(ffty(i))*wqx(4,i),0.d0)		! real and even
	enddo
	call fft_mik1 (n,x,y,qx,ffty,-1,1,1)
	pdw(5,:) = dreal(y(:))										! real and even
	
	pdwsum(:) = pdw(0,:) + pdw(1,:) + pdw(2,:) + pdw(3,:) + pdw(4,:) + pdw(5,:)
	
	end subroutine compute_pdw_2d
	
	subroutine main_iteration (i)
	implicit none
	integer i, i1, i2, nit
	real*8 eps, eps_old, eps_eta, alpha, alpha0, alpha_max, alpha_min
	real*8 mu_eq, p_eq
	real*8 coverage, sigma, eta_eff,rhonew,rescalfactor
	complex*16 y(0:n-1),ffty(0:n-1)
	
	if (i == 0) then
		eta = eta_in
		rhox(:) = eta/(pi*R**2)	  !!!!!!!!!!! EVENTUALLY I HAVE TO MODIFY      = dexp(-beta*V(:) + mu_eq * beta -  pdwsum(:))
		! to add if one wants something oscillating     - 0.5*0.1*dcos(2*pi*x(i)/h_period)
		rho_g = eta/(pi*R**2)
	else
		eta = eta_in + i*(eta_fin-eta_in)/niterations
		rho_g = eta/(pi*R**2)
	endif
	
	print*, 'eta = ', eta
	print*, 'rho_g = ', rho_g
	print*, 'n, nR, dx : ', n,nR,dx
	
	eps_eta = 1.d0
	mu_eq = mu_2d(rho_g)
	
	alpha0 = (1.1-rho_g)**3/1000
	alpha_min = alpha0
	alpha_max = 0.05d0
	
	f = 0.d0
	eps = 1.d0
	nit = 0
	alpha = alpha0
		
		do while (eps > 0.0000001d0)
			nit = nit + 1
			
			! FFT of the density
			do i1=0,n-1
				y(i1) = dcmplx(rhox(i1),0.d0)		! real and even
			enddo
			call fft_mik1 (n,x,y,qx,ffty,1,1,1)
			rhoqx(:) = dreal(ffty(:))				! real and even

			
			! antiFFT to compute the weighted densities
			call compute_weighted_densities_2d ()
!~ 			call print_x (n,wd(1,:))
!~ 			call print_x (n,wd(2,:))
!~ 			call print_x (n,wd(3,:))
!~ 			call print_x (n,wd(4,:))
!~ 			stop
			
			! compute derivative of the FMT functional
			call compute_pd_p_2d ()
!~ 			call print_x (n,pd(0,:))
!~ 			call print_x (n,pd(1,:))
!~ 			call print_x (n,pd(2,:))
!~ 			call print_x (n,pd(3,:))
!~ 			call print_x (n,pd(4,:))
!~ 			call print_x (n,pd(5,:))
			
			! calculate  d Phi^a/ d n^a  *  w_a
			! comes out as \int dx   d Phi^a/ d n^a (x)  w_a (x-r)
			call compute_pdw_2d ()
			
			! surface tension
!~ 			do i1 = nzw,2*nzw-1
!~ 				f(i1) = p(i1)/beta
!~ 			enddo
!~ 			do i1 = 2*nzw, 2*nzw+n_slit
!~ 				f(i1) = (rho(i1)*(dlog(rho(i1)) - 1.d0) + p(i1) )/beta - mu_eq*rho(i1) + p_eq
!~ 			enddo
!~ 			do i1 = 2*nzw+n_slit+1,3*nzw+n_slit
!~ 				f(i1) = p(i1)/beta 
!~ 			enddo
!~ 			sigma = SUM( f(:) ) * dz
			
			rho_new(:) = dexp( -beta*V(:) + beta * mu_eq - pdwsum(:) )  !!!!!!!!!!! HERE I HAVE MODIFIED   old:   = dexp(  beta * mu_eq - pdwsum(:) )
			rhonew = (sum(rho_new(:))*dx)/h_box   !(h_box+1.d0)
			rho_new(:) = rho_new(:) * rho_g/rhonew
			rescalfactor = rho_g/rhonew
			eps = SUM( abs( (rho_new(:) - rhox(:)) / (rho_new(:) + rhox(:)) ))/ n
			
			! Eq. (A) of R. Roth., J.Phys.Cond.Matt. 22, (2010)
			rhox = alpha * rho_new + (1.d0 - alpha) * rhox
			
			if (MOD(nit,2) == 0 .and. nit > 20 ) then
				if (eps_old > eps) then
					alpha = MIN(1.1 * alpha,alpha_max)
				else
					alpha = MAX(alpha / 5.D0, alpha_min)
				endif
			endif
				

			eps_old = eps
			coverage = ( SUM( rhox(:) - rho_g) ) * dx
			
			if (mod(nit,10) == 0) then
				write(*,'(I5,7D15.7)') nit,eps,alpha,sigma,coverage,rhox(0),rhox(n/2),rescalfactor
!~ 				call print_rho ()
!~ 				if (nit.eq.200) stop
			endif
			
			if ( MOD(nit,100000) == 0 ) then
				stop 'too many iterations nit'
			endif

		enddo
		
!~ 		stop
		
		eta_eff = pi*R**2 * (sum(rhox(:))*dx)/h_box   !(h_box+1.d0)
!~ 		print*, eta_eff
!~ 		eps_eta = dabs( (eta_eff-eta)/eta ) 
!~ 		write(*,'(A,F14.6,2es12.4)') "eta_eff, eps_eta: ", eta_eff, eps_eta
!~ 		rho_g = (rho_g/eta_eff)*eta
		
	
	call print_rho ()
	end subroutine main_iteration

	function mu_2d (rr)
	implicit none
	real*8 mu_2d, rr, xx
	xx = rr*pi*R**2
	mu_2d = ( dlog(rr) + 3*xx/(1-xx) + xx*xx/(1-xx)**2 - dlog(1-xx) ) /beta
	end function mu_2d


	
	
end program density_profile_modulated_main_2d