clear
clc
addpath('./src')
%% choose calculation method
global Rosenfeld WhiteBear LJ LR
Rosenfeld = 1;
WhiteBear = 0;
LJ = 1; % turn on attractive tail of interact potential
LR = 0;
LATTICE = 1;
% turn on atractive 9-3(lj) wall-fluid potential
inirhoflag = 0; % 0 without initial density file; 1 with presetted density file
%% definition of global variables
wsig = 0; weps = 0;
global N M 
global R
global LJCUT
global dz
global epsilon
global drho
global TOL
global MAXITER
global BARRIER
BARRIER = 500; %hard wall potential - infinite 
N = 2000; M = 2000; %total number of grid points of x,y axises.
R = 0.5; %particle size sigma =1, raius = 0.5
LJCUT = 2.5; %cutoff length of lj potential
dz = 0.005; %grid size
epsilon = 1.0;% energy unit
drho = 0.0000001; %step size of density for Picard iteration
TOL = 1e-10; %tolerance on convergence of density profile
MAXITER = 10000;     %Maximum number of iterations

global Vext
global rho
global rhonew
global rhokeep
rho=zeros(N,M); rhonew = zeros(N,M); rhokeep = zeros(N,M); Vext = zeros(N,M);

global w0 w1 w2 w3 w1vx w1vy w2vx w2vy %weight functions(omiga_alpha) in FMT, v - vector
w0 = zeros(N,M); w1= zeros(N,M); w2= zeros(N,M); w3= zeros(N,M); 
w1vx= zeros(N,M); w1vy= zeros(N,M); w2vx= zeros(N,M); w2vy= zeros(N,M); 

global n0 n1 n2 n3 n1vx n1vy n2vx n2vy %weighted density in FMT, v - vector
n0 = zeros(N,M);n1= zeros(N,M); n2= zeros(N,M); n3= zeros(N,M); 
n1vx= zeros(N,M); n1vy= zeros(N,M); n2vx= zeros(N,M); n2vy= zeros(N,M);

global dn0 dn1 dn2 dn3 dn1vx dn1vy dn2vx dn2vy 
dn0 = zeros(N,M);dn1= zeros(N,M); dn2= zeros(N,M); dn3= zeros(N,M); 
dn2vx= zeros(N,M); dn2vy= zeros(N,M);
global cder1 cder2x cder2y cder2z c2 c3 c2v dcf
cder1 = zeros(N,M); cder2x = zeros(N,M); cder2y = zeros(N,M); cder2z = zeros(N,M); 
c2 = zeros(N,M);  c3 = zeros(N,M); c2v = zeros(N,M);  dcf = zeros(N,M);

global phiatt cphiatt phi phiid planepot
phiatt = zeros(N,M); cphiatt= zeros(N,M); phi= zeros(N,M); phiid= zeros(N,M); planepot= zeros(N,M);
global mu new_mu dmu alpha rhob etab ew

global p T invT rmin dev old_gamma new_gamma
p=0; T=0; invT =0; rmin=0; dev=0; old_gamma = 0; new_gamma = 0;
global Pi4R2 Pi4R
global iter isweep iend
iter = 0; isweep = 0; iend = 0;
global NiR  
NiR = R/dz; %Number of grid points within particle radius
global NiW
NiW = R/dz;           % Number of grid points within wall
global NiRCUT
NiRCUT=LJCUT/dz;    % Number of grid points within LJ cutoff
global d
d = zeros(N,M);
% main function
%% setting input parameters 
if LR == 1
    in = input('please putin: T alpha ew rhob \n');
    T = in(1);
    alpha = in(2);
    ew = in(3);
    rhob = in(4);% Setting this also sets the chemical potential and pressure
    if WhiteBear == 1
        filename = ['LR-WB wall_sigma',num2str(wsig),'wall_eps',num2str(weps),'.txt'];
    else
        filename = ['LR-RS wall_sigma',num2str(wsig),'wall_eps',num2str(weps),'.txt'];
    end
    fpout = fopen(filename,'w');
else
    in = input('please input: T alpha rhob\n');
    T = in(1);
    alpha = in(2);
    rhob = in(3);% Setting this also sets the chemical potential and pressure
    filename = ['HW.txt'];
    fpout = fopen(filename,'w');
end
mu = chempot();%calculate bulk chemical potential and pressure according to pressure.
p = pressure();
invT = 1/T;
Pi4R2 = pi*4*R^2;
Pi4R = pi*4*R;
etab = rhob*pi/6.;
fprintf(fpout,"\nDFT for a fluid in planar geometry\n");
fprintf(fpout,"sysmtem parameters\nrho_b= %12.10f\neta_b= %12.10f\nmu= %12.10f\nPressure= %12.10f\nTemperature= %10.8f\nInverse Temperature= %10.8f\n",rhob,etab,mu,p,T,invT);
fprintf(fpout,"\n Model parameters:\n");
if WhiteBear ==1 
    fprintf(fpout,"White-Bear Functional adopted\n");
end
if Rosenfeld == 1
    fprintf(fpout,"Rosenfeld FMT Functional adopted\n");
end
if LJ == 1
    fprintf(fpout," LJ interaction turned on\n");
else
    fprintf(fpout," Hard Sphere interaction turned on\n");
end
if LR == 1
    fprintf(fpout,"LJ 9-3 wall-fluids turned on\n");
else
    fprintf(fpout,"hard wall-fluids turned on\n");
end
fprintf(fpout,"dz= %f\nN= %i\nSystem size= %4.2f\nNiR= %i\nNiW= %i\nmixing (alpha)= %4.2f\nTolerance=%e\n",dz,N,dz*N,NiR,NiW,alpha,TOL);

%% preparation for iteration
%formalizm of external potential, initialize density profile and weighted densities
if LJ == 1
    rmin = 2*R*2^(1/6); % with radius=0.5sigma, calculate minimum position of LJ
    fprintf(fpout, "rmin=%f\n Number of grids with in LJcut=%i\n\n",rmin,NiRCUT);
end
if LATTICE == 1
    Vext = importdata("lattice.mat");
else
setVext(wsig,weps);  % Set the external (wall) potential
end
initrho(inirhoflag);  % Set the initial density distribution and print it out
setwhts;%initWeightFunctions;  % Set the weights in fourier space.
w3 = w3*2;
w2 = w2*2;
w2vx = w2vx*2;
w2vy = w2vy*2;
%% LJ attractive potential formalizm
planepot = planepot*0;
if LJ ==1 
%     for i =0:N-1
%         for j = 0:M-1
%             r = sqrt(i^2+j^2)*dz;
%             if r <= LJCUT
%                 planepot(i+1,j+1) = calclinepot(r);
%             if i> 0 && j>0
%                 planepot(N-i+1,M-j+1) = planepot(i+1,j+1);
%             end
%             end
%         end
%     end
%     planepot(find(isnan(planepot)==1)) = 0; 
planepot =  importdata("ljpot.mat");

    %planepot(NiRCUT+1)=planepot(NiRCUT+1)*3/8; planepot(NiRCUT)=planepot(NiRCUT)*7/6;    planepot(NiRCUT-2+1)=planepot(NiRCUT-2+1)*23/24; % Apply the extended quadrature Numerical Recipes Eq. (eq 4.1.14) to deal with the truncation
    %planepot(N-NiRCUT+1)=planepot(N-NiRCUT+1)*3/8; planepot(N-NiRCUT+1+1)=planepot(N-NiRCUT+1+1)*7/6;  planepot(N-NiRCUT+2+1)=planepot(N-NiRCUT+2+1)*23/24;
end
if LJ == 1
    iend = M - 1*NiRCUT; % define ending calculation point to avoid a second wall.
else
    iend = M - 3*NiR; %for hard sphere fluid
end

%% Solvation of functional derivative by Iteration of density profile
disp("Starting Iteration progress...\n");
converged = 0; % flag for iteration, 0 -- not converge;1-- converged.
while converged ==0 && iter   <  MAXITER
    
    %convolution to calculate weighted densities.

    n2=convl(rho,w2);
    n3=convl(rho,w3);
    n2vx=convl(rho,w2vx);
    n2vy=convl(rho,w2vy);

    
    n0 = n2*(Pi4R2^-1);
    n1 = n2*(Pi4R^-1);
    n1vx = n2vx*(Pi4R^-1);
    n1vy = n2vy*(Pi4R^-1);
    %n1vz = n2vz*(Pi4R^-1);
    
    if WhiteBear == 1 && LR ==1 %make sure none zero n3 to prevent WB version phi derivatives from failing
        for i = 1:N
            for j = 1:M
                
                    if n3(i,j)<TOL 
                        n3(i,j) = 1e-12;
                    end
               
            end
        end
    end
    
    %calculate direct correlation function by phi derivatives
    make_dn();
    
    for i = 1:N
        for j = 1:M
          
                cder1(i,j) = dn0(i,j)/Pi4R2 + dn1(i,j)/Pi4R + dn2(i,j); % combined derivatives
                cder2x(i,j) = dn1vx(i,j)/Pi4R + dn2vx(i,j);
                cder2y(i,j) = dn1vy(i,j)/Pi4R + dn2vy(i,j);
                
          
        end
    end
    
    c2=convl(cder1,w2); %componnents of dcf from the convolution with w2
    c2vx = convl(cder2x,w2vx);
    c2vy = convl(cder2y,w2vy);
    %c2vz = convl(cder2z,w2vz);
    c3 =convl(dn3,w3);
    
    for i = 1:N
        for j = 1:M
           
                dcf(i,j) =  -( c2(i,j) + c3(i,j) - (c2vx(i,j) +c2vy(i,j)));%dcf -- direct correlation function.
           
        end
    end
    %calculate LJ acctraction contribution to phi(cphiatt)
    if LJ == 1
        cphiatt=convl(rho,planepot);
    end
    Picardupdate(); %call the Picard-Ng update progress.
    
    if(dev<TOL) 
        converged = 1;
    end
    if(mod(iter,5)==0)
        disp(['Iteration ',num2str(iter),' Deviation=',num2str(dev)]);
        %write_rho();
    end
    
    for i = 1:N
        for j = 1:iend
           
                rho(i,j)=rhonew(i,j);
         
        end
    end
    iter = iter+1;
end%iteration loop ends.

if(converged==1) 
    disp(['Successively converged within tolerance in' ,num2str(iter), 'iterations!']);
else
    disp(["Failed to converge after ",MAXITER," iterations"]);
end

%% Results output
%omega(1);%claculate grand potential with mode 0, or excess potential i.e. surface tension, gamma, with mode 1
fprintf(fpout,"       z        rho(z)      rho(z)/rhob      eta(z)         d[z]  \n\n");
for i = 1:N
    for j = 1:M
        x = (i-NiW-1)*dz;
        y = (j-NiW-1)*dz; 
        fprintf(fpout,"%f  %f  %12.10e  %12.10e  %12.10e  %12.10e\n",x,y,rho(i,j),rho(i,j)/rhob,rho(i,j)*pi/6,d(i,j));
    end
end

save("2d_box_data.mat",'-mat');
fclose(fpout);