% filename rxnOC.m
% This code solves the equations for reaction-diffusion
% in a porous catalyst using the orthogonal collocation method.
% a = 1, 2, 3 for planar, cylindrical, and spherical geometry
% need coll.m, coll_interp.m
% If the problem is linear, set linear = 1 to avoid a second iteration.
global xxcoll aacoll bbcoll qinvcoll wwcoll ncol

geometry = 3
linear = 0
% set ncol before calling %ncol = 1
coll(ncol,geometry,1); % spherical geometry
% set phi before calling %phi = 10
phi2 = phi*phi;

% iterate
tol = 1.e-12;
change = 1.;
iter=0;
% initial guess
for i=1:ncol+1
c(i) = 0.3;
end

while change>tol
iter=iter+1;

% Set up the matrices
% set the matrices to zero
for i=1:ncol+1
f(i) = 0.;
for j=1:ncol+1
aa(i,j) = 0.;
end
end
% set the matrices for the differential equation
for i=1:ncol
for j=1:ncol+1
aa(i,j) = bbcoll(i,j);
end
ans = rate(c(i));
r = ans(1);
dr = ans(2);
aa(i,i) = aa(i,i) - phi2*dr;
f(i) = phi2*(r - dr*c(i));
end
aa(ncol+1,ncol+1) = 1;
f(ncol+1) = 1;

% Solve one iteration.
c11 = aa\f';

% Calculate the criterion to stop the iterations.
sum=0.;
for i=1:ncol+1
sum = sum + abs(c11(i) - c(i));
c(i) = c11(i);
end
if linear==1,break,end
change = sum
end

% Plot the solution.
% Find the coefficients
dd1 = qinvcoll*c';
xp = 0:0.02:1; %plotting points
cinterp = coll_interp(ncol,dd1,xp);
plot(xxcoll,c,'*b',xp,cinterp,'-b')
xlabel('r')
ylabel('c')
c

% calculate the effectiveness factor
sum1 = 0.;
sum2 = 0.;
for i=1:ncol+1
ans = rate(c(i));
r = ans(1);
sum1 = sum1 + wwcoll(i)*r;
sum2 = sum2 + wwcoll(i);
end
eta = sum1/sum2


************************************************
% rate.m
% calculate the rate and derivative for a given concentration

function ans = rate(c)
% use for linear reaction
%ans(1) = c;
%ans(2) = 1.;

% use for second order reaction
%ans(1) = c*c;
%ans(2) = 2.*c;

% use for non-isothermal reaction
%ggg = 18.; bbb = 0.3;
ggg = 30; bbb = 0.4;
t = 1 + bbb - bbb*c;
eee = exp(ggg*(1 - 1/t));
ans(1) = c*eee;
ans(2) = eee*(1 - bbb*ggg*c/(t*t));