% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt%% This program is free software; you can redistribute it and/or% modify it under the terms of the GNU General Public License as% published by the Free Software Foundation; either version 2, or (at% your option) any later version.%% This program is distributed in the hope that it will be useful, but% WITHOUT ANY WARRANTY; without even the implied warranty of% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU% General Public License for more details.%% You should have received a copy of the GNU General Public License% along with this program; see the file COPYING. If not, write to% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,% MA 02111-1307, USA.globalthetabarcafcbfalphakmakmbk1ncolptABQzzrpcabarcbbark1=1e3;%cm^3/mol minkma=0.01/60;%cm/minkmb=0.01/60;%cm/minalpha=1;r=5;%cmcaf=1e-3;%mol/cm^3cbf=1e-3;%mol/cm^3thetabar=10;%minncolpt=20;[z,A,B,Q]=colloc(ncolpt-1,'left');theta=z./(1-z);p=exp(-theta/thetabar)/thetabar;Qz=Q./((1-z).^2);r0=1;rfin=1e-5;nrs=50;rvec=logspace(log10(r0),log10(rfin),nrs)';cabar=zeros(nrs,1);cbbar=zeros(nrs,1);catot=zeros(nrs,1);cbtot=zeros(nrs,1);fori=1:nrsr=rvec(i);x0=[caf*ones(ncolpt,1);cbf*ones(ncolpt,1);caf;cbf];[x,fval,info]=fsolve('particles',x0);fsolve_failed=info<=0;if(fsolve_failed)fprintf('warning, fsolve failed\n');info,x0endca=x(1:ncolpt);cb=x(ncolpt+1:2*ncolpt);cabar(i)=x(2*ncolpt+1);cbbar(i)=x(2*ncolpt+2);inta=Qz'*(ca.*p);intb=Qz'*(cb.*p);catot(i)=(inta*alpha+cabar(i))/(1+alpha);cbtot(i)=(intb*alpha+cbbar(i))/(1+alpha);endy0=[caf;cbf];opts=optimset('MaxFunEvals',2000*numel(x0),...'MaxIter',500*numel(x0));[y,fval,info]=fsolve('cstr',y0,opts);fsolve_failed=info<=0;if(fsolve_failed)fprintf('warning, fsolve failed\n');info,y0endcacstr=y(1);cbcstr=y(2);table1=[rvec*1e41000*catot1000*cacstr*ones(nrs,1)1000*caf*alpha/(1+alpha)*ones(nrs,1)];rvec=[1e-4;1e-3;1e-2];table=[theta];fori=1:length(rvec)r=rvec(i);x0=[caf*ones(ncolpt,1);cbf*ones(ncolpt,1);caf;cbf];[x,fval,info]=fsolve('particles',x0);fsolve_failed=info<=0;if(fsolve_failed)fprintf('warning, fsolve failed\n');info,x0endca=x(1:ncolpt);cb=x(ncolpt+1:2*ncolpt);cabar(i)=x(2*ncolpt+1);cbbar(i)=x(2*ncolpt+2);inta=Qz'*(ca.*p);intb=Qz'*(cb.*p);catot(i)=(inta*alpha+cabar(i))/(1+alpha);cbtot(i)=(intb*alpha+cbbar(i))/(1+alpha);table=[table1000*ca1000*cb];endtable2=[tablezp];savesuspension.dattable1table2;if(~strcmp(getenv('OMIT_PLOTS'),'true'))% PLOTTINGsubplot(2,1,1);semilogx(table1(:,1),table1(:,2:4));% TITLE suspension_1subplot(2,1,2);plot(table2(:,1),table2(:,2:7));% TITLE suspension_2end% PLOTTING