% SIdepTh - Steady state prob. density Intdigtatd Elctrds Tu/30/5/2000
% DEP force - plot for Thesis ch. 6.4  D.Bakewell 15:30 Mo/19/2/2001 
clear all; desktop; Vo=0.1; Vthresh=0.1; d=10; du=d*10^-6; d2=2.*d; 
Xmx=d; Ymx=d; Ximax=50; Yimax=10; Xmn=0.0; Ymn=d/100.; 
Ex=[];Ey=[];tol=[]; trace=[];
% Evaluate DEP and Thermal force coefficient Kf_KbT
eo=8.854e-12; j=sqrt(-1); r = 216e-9/2; %bead radius
Ks=1.2e-9; %bead surface conductance=1.2 nS
Sigp=2*Ks/r; %bead conductivity
Sigm=2e-3; %medium conductivity S/m %Sigm=0.2e-3 fr RO water
ep=2.55*eo; em=78.4*eo; f=2e6; w = 2*pi*f; ecp=ep-j*Sigp/w; ecm=em-j*Sigm/w;
cmf=(ecp-ecm)/(ecp+2*ecm); ImCMf=imag(cmf); ReCMf=real(cmf)
ReCMf=1.;
kBT=298*1.38e-23; kf_kBT=pi*(r^3)*em*ReCMf/kBT 
depthm=kf_kBT/((pi*du)^2) % DEP thermal ratio 
for yi=1:1:Yimax+1; y=Ymn+(yi-1)*(Ymx-Ymn)/Yimax;
 for xi=1:1:Ximax+1; x=Xmn+(xi-1)*(Xmx-Xmn)/Ximax;
  xh=pi*x/d2+pi/4; yh=pi*y/d2; % h = 'hat'
 ThetaT=atan(sin(xh)/sinh(yh))-atan(cos(xh)/sinh(yh));
 ThetaL=log(cosh(yh)+cos(xh))-log(cosh(yh)-cos(xh))+...
    log(cosh(yh)+sin(xh))-log(cosh(yh)-sin(xh));
 Ex(xi,yi)=2.*(Vo/(pi*d))*ThetaT; Ey(xi,yi)=(Vo/(pi*d))*ThetaL;
 ThetaSq(xi,yi)=4.*ThetaT^2+ThetaL^2;
 end % end for xi
end % end for yi
xi=1:1:Ximax+1; x=Xmn+(xi-1)*(Xmx-Xmn)/Ximax;
yi=1:1:Yimax+1; y=Ymn+(yi-1)*(Ymx-Ymn)/Yimax;
ExpPDF=depthm*(Vo^2)*ThetaSq;
Emag=sqrt(((Vo/(pi*d))^2)*ThetaSq);
surf(x,y,Emag');dummy=menu('E magntude      ','Press to continu');close
surf(x,y,ExpPDF');dummy=menu('ExpPDF        ','Press to continu');close
if Vo<=Vthresh  
intconst=dblquad('Isidep',Xmn,Xmx,Ymn,Ymx,tol,'quad8')
SSPDF=exp(ExpPDF)/intconst; % Steady Starte p d f 
surf(x,y,SSPDF'); dummy=menu('SSPDF        ','Press to continu');close
else
Un_NrmSSPDF=exp(ExpPDF); % Steady Starte p d f   
surf(x,y,Un_NrmSSPDF'); dummy=menu('Un_normlisd SSPDF','Press to continu');close
end % end if  
%surf(x,y,Ex'); dummy=menu('Ex            ','Press to continu');close;
%surf(x,y,Ey'); dummy=menu('Ey            ','Press to continu');close;
xi=1:1:Ximax+1; XVec(xi)=Xmn+(xi-1)*Xmx/Ximax; yi=1;
EVec(xi,1)=Emag(xi,yi); plot(XVec,EVec);
dummy=menu(' E magntde           ','Press to continu');close;
ExpPDFVec(xi,1)=ExpPDF(xi,yi); plot(XVec,ExpPDFVec);
dummy=menu(' ExpPDF           ','Press to continu');close;
SspdfVec(xi,1)=SSPDF(xi,yi); plot(XVec,SspdfVec);
dummy=menu('SSpdf at Lower boundry','Press to continu');close;
SspdfVec(xi,1)=SSPDF(xi,Yimax+1); plot(XVec,SspdfVec);
dummy=menu('SSpdf at Upper boundry','Press to continu');close;
%ExVec(xi,1)=Ex(xi,yi); plot(XVec,ExVec);
%dummy=menu(' Ex           ','Press to continu');close;
%EyVec(xi,1)=Ey(xi,yi); plot(XVec,EyVec);
%dummy=menu(' Ey           ','Press to continu');close;
if Vo<=Vthresh;  
IntSSPDF=dblquad('cksidep',Xmn,Xmx,Ymn,Ymx,tol,'quad8')% should be unity!!
end % end if
cd e:
close;