% 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.
% Try to solve for the dimensionless pellet profile using
% a Langmuir Hinshelwood mechanism.
% d^2 c/ d r^2 -Phi^2 c / (1 + phi c) = 0
% c(1) = 1
% dc/dr (0) = 0
% goal is to not use a shooting method to avoid numerical problems
% try it as a pure ODE
% jbr, 6/12/01
% first effort: consider dimensional problem
% x1 = c_A
% x2 = dc_A/dr
% dx2/dr = p1 p2 c_A / (1 + p2 c_A)
% dx1/dr = x2
% x1(0) = c0
% x2(0) = 0
% solve this IVP to r=r1 and x1=c1
% p1 = k cm/De
% p2 = KA
% Phi^2 = p1*p2*r1^2
% phi = p2*c1
% problem up to here is both Phi and phi will change with r1,c1
% make p2 change with time as an ODE as well, to keep phi=constant
% dp2/dr = -p2*x2/x1
% p2(0) = phi/x1(0)
% OK, don't think the above idea of changing p2 is solving the right problem.
% Just march the ODE to desired phi, change p2 and solve ODE again
% to same phi. The two values of Phi should be different, however,
% and that allows a family of Phi at constant phi. Drawback is that
% we are solving the ODE a lot.
% Revised 8/22/2018
p1 = 1;
p2 = 1;
c0 = 0.1;
p.p1 = p1; p.p2 = p2;
x0 = [c0; 0; 0; 0];
npts = 200;
rstop = 10;
rsteps = linspace(0.001,rstop,npts)';
phivec = [0.1 10 100 1000];
nphi = length(phivec);
results = ;
for j = 1:nphi
% March ODE solver out to phi stopping condition for a range of p2
% next pairs go together for smooth log log plot of etavec vs. Phivec
phi = phivec(j);
p.phi = phi;
p2vec = [1e-40 1e-10 .1 1 10 20 30 40 50 55 60 65 70 75 80 90 95 99 99.5];
if (phi >= 100)
p2vec(1) = 1e-20;
np2 = length(p2vec);
c0 = 1e-2*phi;
Phivec = zeros(np2,1);
etavec = zeros(np2,1);
for i = 1:np2
p2 = p2vec(i);
p.p2 = p2;
rstop = max(100,100*sqrt(1/p2));
rsteps = [0; linspace(0.001,rstop,npts)'];
x0 = [c0; 0; 0; 0];
opts = odeset('Events', @(t, x) g(t,x,p), ...
'AbsTol', sqrt(eps), 'RelTol', sqrt(eps));
[rout,x] = ode15s (@(t,x) hougensrt(t,x,p), rsteps, x0, opts);
if ( length(rout) == npts + 1)
fprintf ('hey, did not reach requested phi value, increase rstop\n');
nout = length(rout);
cbar = x(:,1)/x(nout,1);
rbar = rout/rout(nout);
% plot (rbar,cbar);
Phivec(i) = sqrt(p1*p2)*rout(nout);
etavec(i) = (1 + phi)*x(nout,2)*rout(nout)/x(nout,1)/Phivec(i)^2;
% check, should be zero if solved the dim'less model
% hmm, apparently cannot trust the last xdot. Could make a call
% to hougensrt to correct this last one if this were important
% check passed, 6/13/01
% ode15s doesn't return xdot at all the solution times, so would
% need to compute it here to check.
% check = rout(nout)^2/x(nout,1)*xdot(:,2) - ...
% fprintf ('norm of error in solution: %g\n', max(abs(check(1:nout-1))));
Phiprimevec = phi/(1 + phi) / sqrt( 2*(phi - log(1+phi)) ) * Phivec;
results = [results Phivec Phiprimevec etavec];
save -ascii etahougwat.dat results;
if (~ strcmp (getenv ('OMIT_PLOTS'), 'true')) % PLOTTING
subplot (2, 1, 1);
loglog (results(:,1), results(:,3), ...
results(:,4), results(:,6), ...
results(:,7), results(:,9), ...
axis ([0.07, 100, .08, 1.5]);
% TITLE etahougwat
subplot (2, 1, 2);
loglog (results(:,2), results(:,3), ...
results(:,5), results(:,6), ...
results(:,8), results(:,9), ...
axis ([0.07, 10, .08, 1.5]);
% TITLE etahougwatasy
end % PLOTTING