module density_profile_modulated_def_2d
implicit none
	real*8, parameter ::  R = 0.5d0						! particle radius
	real*8, parameter ::  Amp = 0.45d0					! amplitude of the force
	real*8, parameter ::  h_period = 0.45d0				! period
	integer, parameter :: n_period = 3					! not one because we need more to have a good resolution of q at low values
	real*8, parameter ::  h_box = n_period*h_period		! box size
	real*8, parameter ::  eta_in	= 0.55d0			! initial packing fraction (in the bulk)
	real*8, parameter ::  eta_fin = 0.55d0				! final packing fraction (in the bulk)
	integer, parameter :: niterations = 0				! number of interations. From one iteration to the other phi increases untill the it reaches phi_fin
	integer, parameter :: ndft = 1						! different type of density profile theory (1 == Eq.(9) of Roth et al. J. Chem. Phys. 136, 081101 (2012) )
	
	real*8, parameter ::  pi = 3.14159265358979323846264338327950288d0
	
	integer, parameter :: n = 4096								! number of points in the direction x (parallel to the modulation) must be a power of 2!
	real*8, parameter ::  dx = h_box/n							! discretization in the x direction
	real*8, parameter ::  dqx = 1.d0/h_box						! la rho avra` valori solo per: 2*pi / h_period
	integer, parameter :: nR = INT(R/dx + 0.0000000001d0)		! if -nR <= i <= nR   -->   -R <= x <= R

	real*8, parameter ::  beta = 1.D0					! inverse temperature
	real*8, parameter ::  eps_old = 1.d0				! ?
	
	
	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	!!!!!!!!!!!!!!!!   GLOBAL QUANTITES !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	real*8 x(0:n-1)
	real*8 qx(0:n-1)
	real*8 wx(4,0:n-1)
	real*8 wqx(4,0:n-1)
	real*8 rhox(0:n-1)
	real*8 rhoqx(0:n-1)
	real*8 wd(4,0:n-1)
	
	real*8 :: eta, rho_g, rho_g_in
	real*8 :: pd(0:5,0:n-1), pdw(0:5,0:n-1)
	real*8 :: p(0:n-1), f(0:n-1), p1(0:n-1), pd1(0:n-1), p2(0:n-1), pd2(0:n-1), pdwsum(0:n-1)
	real*8 :: V(0:n-1)
	complex*16 :: h1(0:n-1), h2(0:n-1),h3(0:n-1)
	real*8 :: rho_new(0:n-1)
	
	save
	

	
contains

	Subroutine print_rho ()
	implicit none
	integer i
	open (21, file='rhoV_n.dat', status='unknown')
!~ 	call contascorri(i,21)
!~ 	write (unit=21, fmt="(a,f12.8)") "# rho_bulk :", rho_g
!~ 	write (unit=21, fmt="(a,f12.8,a,f4.1,a,f5.3)") "# packing fraction :", eta, "   period :",h_period, "   Amp :",Amp
	do i = 0,n-1
		WRITE (21,'(3es14.6)') x(i), rhox(i), V(i)
	enddo
!~ 	write(21,*) '  '
	close (21)
	end Subroutine print_rho
	
	Subroutine print_x (n,ff)
	implicit none
	integer i,n
	real*8 ff(0:n-1)
	do i = 0,n-1,8
		write (*,'(I4,A1,8(1X,F10.4))') i,')', &
		 ff(i),ff(i+1),ff(i+2),ff(i+3),ff(i+4),ff(i+5),ff(i+6),ff(i+7)
	enddo
	write(*,*) '-------------------------------------------------------'
	
	end Subroutine print_x
	
	Subroutine define_grid ()
	implicit none
	integer i
	do i = 0,n-1
		x(i) = -n/2*dx + i*dx
		qx(i) = i*dqx
	enddo
	end Subroutine define_grid

	Subroutine define_weight_planar_2d ()
    implicit none
    integer i
	complex*16 y(0:n-1),ffty(0:n-1)
	! 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]
	wx(:,:) = 0.d0
	wx(1,n/2) = 2*R
	wx(2,n/2) = 2.d0
	wx(3,n/2) = 0.d0
	wx(4,n/2) = 0.d0
    do i = n/2+1,n/2+nR
      wx (1,i) = 2*dsqrt(R**2 - x(i)**2)				! real and even
      wx (2,i) = 2*R/ dsqrt(R**2 - x(i)**2)				! real and even
      wx (3,i) = 2*x(i)/ dsqrt(R**2 - x(i)**2)			! real and odd
      wx (4,i) = 2*x(i)*x(i)/R /dsqrt(R**2 - x(i)**2)	! real and even
      wx (1,n-i) =  wx(1,i)
      wx (2,n-i) =  wx(2,i)
      wx (3,n-i) = -wx(3,i)
      wx (4,n-i) =  wx(4,i)
      ! cure divergences
      if (i.eq.n/2+nR .and. x(i).ge.R) then
		wx (2,i) = 0.d0
		wx (2,n-i) = 0.d0
		wx (3,i) = 0.d0
		wx (3,n-i) = 0.d0
		wx (4,i) = 0.d0
		wx (4,n-i) = 0.d0
      endif
    enddo
    
    ! and theirs FT  (FT(f) = int f(x) exp (2 pi i q x) dx
	wqx (1,0) = pi*R/2
	wqx (2,0) = 2*pi*R										
	wqx (3,0) = 0.d0									
	wqx (4,0) = pi*R
    do i = 1,n/2-1
		wqx (1,i) = R*BESSJ1(2*pi*qx(i)*R) / qx(i)						! real and even
		wqx (2,i) = 2*pi*R*BESSJ0(2*pi*qx(i)*R)							! real and even
		wqx (3,i) = 2*pi*R*BESSJ1(2*pi*qx(i)*R)							! imaginary and odd
		wqx (4,i) = 2*pi*R*BESSJ0(2*pi*qx(i)*R)-BESSJ1(2*pi*qx(i)*R)/qx(i)	! real and even
		! we have to symmetrize to take into account symmetries
		wqx (1,n-i) = wqx (1,i)
		wqx (2,n-i) = wqx (2,i)
		wqx (3,n-i) = -wqx (3,i)
		wqx (4,n-i) = wqx (4,i)
	enddo
	wqx (1,n/2) = R*BESSJ1(2*pi*qx(n/2)*R) / qx(n/2)					! real and even
	wqx (2,n/2) = 2*pi*R*BESSJ0(2*pi*qx(n/2)*R)							! real and even
	wqx (3,n/2) = 0.													! imaginary and odd
	wqx (4,n/2) = 2*pi*R*BESSJ0(2*pi*qx(n/2)*R)-BESSJ1(2*pi*qx(n/2)*R)/qx(n/2)	! real and even
	
	! better to do the FFT because with this function there are problems...
!~ 	do i=0,n-1
!~ 		y(i) = dcmplx(wx(1,i),0.d0)						! real and even
!~ 	enddo
!~ 	call fft_mik1 (n,x,y,qx,ffty,1,1,1)
!~ 	wqx(1,:) = dreal(ffty(:))							! real and even
	
!~ 	do i=0,n-1
!~ 		y(i) = dcmplx(wx(2,i),0.d0)						! real and even
!~ 	enddo
!~ 	call fft_mik1 (n,x,y,qx,ffty,1,1,1)
!~ 	wqx(2,:) = dreal(ffty(:))							! real and even
	
!~ 	do i=0,n-1
!~ 		y(i) = dcmplx(wx(3,i),0.d0)						! real and odd
!~ 	enddo
!~ 	call fft_mik1 (n,x,y,qx,ffty,1,0,1)
!~ 	wqx(3,:) = dimag(ffty(:))							! imaginary and odd
	
!~ 	do i=0,n-1
!~ 		y(i) = dcmplx(wx(4,i),0.d0)						! real and even
!~ 	enddo
!~ 	call fft_mik1 (n,x,y,qx,ffty,1,1,1)
!~ 	wqx(4,:) = dreal(ffty(:))							! real and even
	
	end Subroutine define_weight_planar_2d

	Subroutine define_potential ()
	implicit none
	integer i
	do i = 0,n-1
		V(i) = Amp*dcos(2*pi*x(i)/h_period)
	enddo
	end Subroutine define_potential

	FUNCTION BESSJ0 (X)
	IMPLICIT NONE
	REAL*8 X,BESSJ0,AX,FR,FS,Z,FP,FQ,XX
!	This subroutine calculates the First Kind Bessel Function of
!	order 0, for any real number X. The polynomial approximation by
!	series of Chebyshev polynomials is used for 0<X<8 and 0<8/X<1.
!	REFERENCES:
!	M.ABRAMOWITZ,I.A.STEGUN, HANDBOOK OF MATHEMATICAL FUNCTIONS, 1965.
!	C.W.CLENSHAW, NATIONAL PHYSICAL LABORATORY MATHEMATICAL TABLES,
!	VOL.5, 1962.
	REAL *8 Y,P1,P2,P3,P4,P5,R1,R2,R3,R4,R5,R6,Q1,Q2,Q3,Q4,Q5,S1,S2,S3,S4,S5,S6
	DATA P1,P2,P3,P4,P5 /1.D0,-.1098628627D-2,.2734510407D-4, &
      -.2073370639D-5,.2093887211D-6 /
	DATA Q1,Q2,Q3,Q4,Q5 /-.1562499995D-1,.1430488765D-3, &
      -.6911147651D-5,.7621095161D-6,-.9349451520D-7 /
	DATA R1,R2,R3,R4,R5,R6 /57568490574.D0,-13362590354.D0, &
      651619640.7D0,-11214424.18D0,77392.33017D0,-184.9052456D0 /
	DATA S1,S2,S3,S4,S5,S6 /57568490411.D0,1029532985.D0, &
      9494680.718D0,59272.64853D0,267.8532712D0,1.D0 /
	IF(X.EQ.0.D0) GO TO 1
	AX = ABS (X)
	IF (AX.LT.8.) THEN
		Y = X*X
		FR = R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6))))
		FS = S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6))))
		BESSJ0 = FR/FS
	ELSE
		Z = 8./AX
		Y = Z*Z
		XX = AX-.785398164
		FP = P1+Y*(P2+Y*(P3+Y*(P4+Y*P5)))
		FQ = Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*Q5)))
		BESSJ0 = SQRT(.636619772/AX)*(FP*COS(XX)-Z*FQ*SIN(XX))
	ENDIF
	RETURN
	1 BESSJ0 = 1.D0
	RETURN
	END FUNCTION BESSJ0

	FUNCTION BESSJ1 (X)
	IMPLICIT NONE
	REAL *8 X,BESSJ1,AX,FR,FS,Z,FP,FQ,XX
!	This subroutine calculates the First Kind Bessel Function of
!	order 1, for any real number X. The polynomial approximation by
!	series of Chebyshev polynomials is used for 0<X<8 and 0<8/X<1.
!	REFERENCES:
!	M.ABRAMOWITZ,I.A.STEGUN, HANDBOOK OF MATHEMATICAL FUNCTIONS, 1965.
!	C.W.CLENSHAW, NATIONAL PHYSICAL LABORATORY MATHEMATICAL TABLES,
!	VOL.5, 1962.
	REAL *8 Y,P1,P2,P3,P4,P5,P6,R1,R2,R3,R4,R5,R6  &
               ,Q1,Q2,Q3,Q4,Q5,S1,S2,S3,S4,S5,S6
	DATA P1,P2,P3,P4,P5 /1.D0,.183105D-2,-.3516396496D-4,  &
      .2457520174D-5,-.240337019D-6 /,P6 /.636619772D0 /
	DATA Q1,Q2,Q3,Q4,Q5 /.04687499995D0,-.2002690873D-3,   &
      .8449199096D-5,-.88228987D-6,.105787412D-6 /
	DATA R1,R2,R3,R4,R5,R6 /72362614232.D0,-7895059235.D0, & 
      242396853.1D0,-2972611.439D0,15704.48260D0,-30.16036606D0 /
	DATA S1,S2,S3,S4,S5,S6 /144725228442.D0,2300535178.D0, &
      18583304.74D0,99447.43394D0,376.9991397D0,1.D0 /

	AX = ABS(X)
	IF (AX.LT.8.) THEN
	Y = X*X
	FR = R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6))))
	FS = S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6))))
	BESSJ1 = X*(FR/FS)
	ELSE
	Z = 8./AX
	Y = Z*Z
	XX = AX-2.35619491
	FP = P1+Y*(P2+Y*(P3+Y*(P4+Y*P5)))
	FQ = Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*Q5)))
	BESSJ1 = SQRT(P6/AX)*(COS(XX)*FP-Z*SIN(XX)*FQ)*SIGN(S6,X)
	ENDIF
	RETURN
	END FUNCTION BESSJ1

	subroutine fft_mik1 (n,x,y,q,ffty,nrev,even,reale)
	! even e reale si riferiscono sempre alla funzione nello spazio delle posizioni
	implicit none
	integer nrev,n,even,reale
	real*8 x(0:n-1)
	complex*16 y(0:n-1),ffty(0:n-1)
	real*8 q(0:n-1)
	real*8 Delta
	complex*16 g(0:n-1),h(0:n-1)
	integer i
	
	if (nrev .eq. 1) then
		Delta = x(1)-x(0)
		g(:) = y(:)
		call fft_c (n,g,h,nrev)
		do i = 0,n-1
			if (even.eq.1 .and. reale.eq.1) then
				ffty(i)=dcmplx(Delta*dreal(h(i))*(-1)**i,0.d0)
			else if (even.eq.1 .and. reale.eq.0) then
				ffty(i)=dcmplx(0.d0,Delta*dimag(h(i))*(-1)**i)
			else if (even.eq.0 .and. reale.eq.1) then
				ffty(i)=dcmplx(0.d0,Delta*dimag(h(i))*(-1)**i)
			else if (even.eq.0 .and. reale.eq.0) then
				ffty(i)=dcmplx(Delta*dreal(h(i))*(-1)**i,0.d0)
			else
				ffty (i) = Delta*h(i)*(-1)**i
			endif
		enddo
	else
		Delta = q(1)-q(0)
		do i = 0,n-1
			g(i) = ffty(i)*(-1)**i
		enddo
		call fft_c (n,g,h,nrev)		
		do i = 0,n-1
			if (even.eq.1 .and. reale.eq.1) then
				y(i)=dcmplx(Delta*dreal(h(i)),0.d0)
			else if (even.eq.1 .and. reale.eq.0) then
				y(i)=dcmplx(0.d0,Delta*dimag(h(i)))
			else if (even.eq.0 .and. reale.eq.1) then
				y(i)=dcmplx(Delta*dreal(h(i)),0.d0)
			else if (even.eq.0 .and. reale.eq.0) then
				y(i)=dcmplx(0.d0,Delta*dimag(h(i)))
			else
				y(:) = h(i)*Delta
			endif
		enddo
	endif
	
	end subroutine fft_mik1

	subroutine fft_c (n,g,h,nrev)
	implicit none
	INTEGER, INTENT (IN) :: n,nrev
	INTEGER :: i,j,m,nstep
	COMPLEX*16, DIMENSION (0:n-1) :: g
	COMPLEX*16, DIMENSION (0:n-1), INTENT (OUT) :: h
	COMPLEX*16 :: c1,ex1,ex2
	! Bit reversal section
	j = 0          ! variable for position of element to be exchanged
	do i = 0,n-1
		if (j > i) then
			c1 = g(j)
			g(j) = g(i)
			g(i) = c1
		endif
	! Determination of j: Use previous j and add bits from the "left" side
	!                     In case of overflow subtract the bits
		m = n/2
		do while (m > 1 .AND. j >= m)
			j = j - m
			m = m/2
		enddo
		j = j + m
	enddo
	! Danielson-Lanczos
	nstep = 1
	do while (nstep < n)
		do i = 0, nstep-1
			ex1 = EXP( nrev*(i)*(0,1)*pi/nstep )
			ex2 = EXP( nrev*(i+nstep)*(0,1)*pi/nstep )
			do j = i, n-1, 2*nstep
				h(j) = g(j) + g(j+nstep) * ex1
				h(j+nstep) = g(j) + g(j+nstep) * ex2
			enddo
		enddo
		nstep = 2*nstep
		g = h
	enddo
	end subroutine fft_c

	subroutine contascorri (N,a)
	implicit none
	! conta le linee N del file aperto alla unità a
	integer N,a,status
	character (len=100) line
	N=0
	do
		read (a,*,iostat=status) line
		if (status.lt.0) exit
		N = N+1
	enddo
	end subroutine contascorri

end module density_profile_modulated_def_2d
