disp('partcl');
d=2;
so=1;
zo=.001;
u=10;
us=.035*u;
k=.4;
r=15E-6;
D=1.3E-11/(r*100)*1e-4;
if r<4.5e-5;
	w=1.38E8*(r)^2;
else;
	w=(r-1e-5)*8000;
end;
b=w/(k*us);
Sc=1E10*(r*100);
St=4.4E3*(r*100*us*100)^2;
St=max(.10,St);
vt=us/log(10/zo)/10^(3/St);
kd=us/log(10/zo)/sqrt(Sc);
no=so/(vt+kd+w);
nd=(so/w+(no-so/w)/(d/zo)^b);
for i=1:50;
	j=(i)/10-3;
	z=10^(j);
	zz(i)=z;
	if z<d;
		n(i)=so/w+(no-so/w)/(z/zo)^b;
	else;
		n(i)=(so/w+(no-so/w)/(d/zo)^b)/(z/d)^b;
	end;
	if i>1;
		dn(i)=(n(i)-n(i-1))/(zz(i)-zz(i-1))*k*(zz(i)+zz(i-1))/2*us+w*(n(i)+n(i-1))/2+.01;
	else;
		dn(i)=0;
	end;
end;
loglog(n,zz);
%axis([.0,4,.01,10]);
end;