clear; clc;
c=1e-4;er=80;eo=8.85e-12;T=300;q=1.6e-19;k=1.38e-23;Na=6.02e23;
AE=1e-4; % cm^2
AS=1e-6; % mm^2
alpha=q*Na/(er*eo);betasq=2*q*c/(k*T);
d=0.1; % 10 cm
vo=1; % 1 m/s
LS=0.001; % 1 mm
B=1.0; % Tesla
basqrt=(betasq*alpha)^0.5;
x=linspace(-d,d,1000);
xhi=linspace(0,0.95*d,100);
xlo=linspace(-0.95*d,0,100);
xs=linspace(0.9999*d,d,200);
DVE=-(4/3)*vo*B*d/(1++2*AS/(LS*er*AE*basqrt));
sgE=DVE*eo*AS/(LS*AE);
dcs=2*vo*B.*xs/(d^2*alpha)+(betasq*sgE/(eo*er)+2*vo*B/(d^2*alpha))/basqrt*((exp(-basqrt.*(xs+d))-exp(basqrt.*(xs-d)))/(1+exp(-2*basqrt*d)));
dchi=2*vo*B.*xhi/(d^2*alpha)+(betasq*sgE/(eo*er)+2*vo*B/(d^2*alpha))/basqrt*((exp(-basqrt.*(xhi+d))-exp(basqrt.*(xhi-d)))/(1+exp(-2*basqrt*d)));
dclo=2*vo*B.*xlo/(d^2*alpha)+(betasq*sgE/(eo*er)+2*vo*B/(d^2*alpha))/basqrt*((exp(-basqrt.*(xlo+d))-exp(basqrt.*(xlo-d)))/(1+exp(-2*basqrt*d)));
Exs=-sgE/(eo*er)+vo*B*(xs.^2/d^2-1)+(sgE/(eo*er)+2*vo*B/(basqrt*d)^2)*(1-(exp(basqrt.*(xs-d))+exp(-basqrt.*(xs+d)))/(1+exp(-2*basqrt*d)));
Ex=-sgE/(eo*er)+vo*B*(x.^2/d^2-1)+(sgE/(eo*er)+2*vo*B/(basqrt*d)^2)*(1-(exp(basqrt.*(x-d))+exp(-basqrt.*(x+d)))/(1+exp(-2*basqrt*d)));

figure(1);
subplot(4,1,1);
plot(xs,dcs/c)
xlabel('d (m)')
ylabel('deltaC/cbar')
subplot(4,1,2);
plot(xhi,dchi/c,xlo,dclo/c)
xlabel('d (m)')
ylabel('deltaC/cbar')
subplot(4,1,3);
plot(xs,Exs)
xlabel('d (m)')
ylabel('Ex (V/m)')
subplot(4,1,4);
plot(x,Ex)
xlabel('d (m)')
ylabel('Ex (V/m)')

S=['DeltaVE = ',num2str(DVE)];
disp(S)




