%% Solve the Poisson Boltzmann Equation Using the Shooting Method

sigma_um=10^6;              %Surface charge density: number per micron^2 (adjust)

Cinf=   10^-1;              %Salt concentration at infinity (M)

Na=6.03*10^23;              %Avagadro's number
C=1.6*10^-19;               %Charge on electron in Coulomb
KbT=1.38*10^-23*300         %KbT at 300 K
epsilon=80*8.85*10^-12      %Permittivity of water
lB=C^2/(4*pi*epsilon*KbT)   %Bjerrum length
Csi=Cinf*Na*10^3;           %Salt concentration in number/m^3
lD=(8*pi*lB*Csi)^-(1/2)     %Deybe screening length
sigma=sigma_um*10^12*C;     %Surface charge density: charge per m^2

%make x in nm
en=10*(10^-9);
h=.01*(10^-9);
x=(0:h:en)';

%solve as two first order ODE
%W=dVdx
%dWdx=1/lD^2*sinh(V(x))

%This is the Debye-Hueckel Result
V_DH=sigma*lD/epsilon*exp(-x/lD);
Vbar_DH=C/KbT*sigma*lD/epsilon*exp(-x/lD);

figure
plot(x*10^9,V_DH,'linewidth',3)
hold on

V=zeros(length(x),1);
W=zeros(length(x),1);
V(1)=Vbar_DH(1);
V(end)=Inf; %Run once to get this value
W(1)=-4*pi*lB*sigma/C;
tol=Vbar_DH(1)/10000;
Vendold=Inf; 
dV1=V(1)/10;
V1old=V(1)-dV1;

while abs(V(end))>tol
    %shoot
    for ii=2:length(x)
        %Euler method
        V(ii)=V(ii-1)+h*W(ii-1);
        W(ii)=W(ii-1)+h*(lD^-2*sinh(V(ii-1))); %Planar PB 
        %W(ii)=W(ii-1)+h*(lD^-2*(V(ii-1))); %Planar DH
    end
    
    V1old=V(1);
    V(end)
    
    %adjust V(1)
    if (V(end)>0)&(Vendold>0)
        dV1=-abs(dV1);
        V(1)=V(1) + dV1;
        Vendold=V(end);
        continue
    end
    
    if (V(end)>0)&(Vendold<0)
        dV1=-0.5*dV1;
        V(1)=V(1) + dV1;
        Vendold=V(end);
        continue
    end
    
    if (V(end)<0)&(Vendold<0)
       dV1=abs(dV1);
        V(1)=V(1) + dV1;
        Vendold=V(end);
        continue
    end
    
    if (V(end)<0)&(Vendold>0)
        dV1=-0.5*dV1;
        V(1)=V(1) + dV1;
        Vendold=V(end);
        continue
    end
    
end

V_PB=V*KbT/C;

%This is the analytical solution to the PB equation
xp=lD*log(C/(2*pi*lB*sigma*lD)*(1+sqrt(1+(2*pi*lB*sigma*lD/C)^2)));
Vtheory=2*log( (1+exp(-(x+xp)/lD))./(1-exp(-(x+xp)/lD)));

plot(x/10^-9,V_PB,'r','linewidth',3)
plot(x/10^-9,Vtheory*KbT/C,'g--','linewidth',2)

s=sprintf('Plane with \\sigma=%g e/\\mum^2',sigma_um)
title(s,'fontsize',14)
xlabel('Distance from Plane (nm)','fontsize',14)
ylabel('Electric Potential (V)','fontsize',14)
legend('V(r) Deybe-Hueckel','V(r) Poisson-Boltzmann','V(r) Poisson-Boltzmann (Analytical)')

%% Plot the ion concentrations

%calculate ion concentrations
cplus_PB=Cinf*exp(-C/KbT*V_PB);
cminus_PB=Cinf*exp(C/KbT*V_PB);

cplus_DH=Cinf*exp(-C/KbT*V_DH);
cminus_DH=Cinf*exp(C/KbT*V_DH);

%plot
figure
semilogy(x/10^-9,cplus_PB,'r','linewidth',2)
hold on
semilogy(x/10^-9,cminus_PB,'r--','linewidth',2)

semilogy(x/10^-9,cplus_DH,'b','linewidth',2)
semilogy(x/10^-9,cminus_DH,'b--','linewidth',2)

s=sprintf('Plane with \\sigma=%2.1g e/\\mum^2',sigma_um)
title(s,'fontsize',14)
xlabel('Distance from Plane (nm)','fontsize',14)
ylabel('Electric Potential (V)','fontsize',14)
legend('c_+(r) Poisson-Boltzmann (M)','c_-(r) Poisson-Boltzmann (M)','c_+(r) Deybe-Hueckel (M)','c_-(r) Deybe-Hueckel (M)')









