%% 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 in M (100 mM)
R=50*10^-9;                 %Radius of Sphere in meters (50 nm)

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
Q=4*pi*R^2*sigma;           %Total charge on sphere

%make r in nm
en=60*(10^-9);
h=.01*(10^-9);
r=(R:h:en)';

%solve as two first order ODEs
%W=dVdr
%dWdr=1/lD^2*sinh(V(x))-2*W/r

%This is the Debye-Hueckel Result
A=Q/(4*pi*epsilon*(1+R/lD)); %*exp(R/lD) which is combined below
V_DH=A./r.*exp((R-r)/lD);
Vbar_DH=C/KbT*V_DH;

figure
plot(r/10^-9,V_DH,'linewidth',3)
hold on

V=zeros(length(r),1);
W=zeros(length(r),1);
V(1)=Vbar_DH(1);
V(end)=Inf; %Run once to get this value
W(1)=-(C/KbT)*Q/(4*pi*epsilon*R^2);
tol=Vbar_DH(1)/10000;
Vendold=Inf; 
dV1=V(1)/10;
V1old=V(1)-dV1;


while abs(V(end))>tol
    
    %shooting method: http://en.wikipedia.org/wiki/Shooting_method
    for ii=2:length(r)
        %Euler method
        V(ii)=V(ii-1)+h*W(ii-1);
        W(ii)=W(ii-1)+h*(1/lD^2*sinh(V(ii-1))-2*W(ii-1)/r(ii-1) ); %Spherical PB 
        %W(ii)=W(ii-1)+h*(1/lD^2*(V(ii-1))-2*W(ii-1)/r(ii-1) ); %Spherical DH
    end
    
    if isnan(V(end)) %This sets NaN values to the extreme value-needed for updates to V(1)
        V(end)=max(V)+min(V);   
    end
    
    V1old=V(1);

    %adjust V(1) to converge to V(end)=0
    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;

plot(r/10^-9,V_PB,'r','linewidth',2)
s=sprintf('R=%g nm \\sigma=%g e/\\mum^2',R*10^9,sigma_um)
title(s,'fontsize',14)
xlabel('Distance from center of sphere (nm)','fontsize',14)
ylabel('Electric Potential (V)','fontsize',14)
legend('V(r) Deybe-Hueckel','V(r) Poisson-Boltzmann')

%% 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(r/10^-9,cplus_PB,'r','linewidth',2)
hold on
semilogy(r/10^-9,cminus_PB,'r--','linewidth',2)

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

s=sprintf('R=%g nm \\sigma=%g e/\\mum^2',R*10^9,sigma_um)
title(s,'fontsize',14)
xlabel('Distance from center of sphere (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)')


