% 	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


function coll(nx,aaa,kon)
% coll - for orthogonal collocation with symmetric polynomials
% nx+1 = total number of collocation points
% aaa = geometry factor (1 = planar, 2 = cylindrical, 3 = spherical)
% kon = 0, will print matrices
% xxcoll(I) = location of collocation points, X(1) >0, X(N+1) = 1.0
% The collocation points for nx = 1 are for W = 1-x^2
% The collocation points for nx > 1 are for W = 1
% returned output is:
%	xxcoll(i), the collocation points
%	aacoll(i), the matrix for first derivative
%	bbcoll(i), the matrix for the Laplacian
%	qinvcoll(i), the matrix for the inverse of q

global xxcoll aacoll bbcoll qinvcoll wwcoll

if(nx<1) 
  nx=1;
  fprintf('\n ***************************************** \n')
  fprintf('The number of interior collocation points must be between \n')
  fprintf('1 and 6 (inclusive).  nx has been changed to %g\n\n',nx)
end
if(nx>6)
  nx=6;
  fprintf('\n ***************************************** \n')
  fprintf('The number of interior collocation points must be between \n')
  fprintf('1 and 6 (inclusive).  nx has been changed to %g\n\n',nx)
end

% now the code
n1 = nx+1;
%fprintf('The total number of collocation points (symmetric) is %g\n\n',n1)
%fprintf('The geometry factor is (1=planar, 2=cyl., 3=sph.) %g\n\n',aaa)
xxcoll(n1) = 1.0;

% find the collocation points
if (aaa==3)    %spherical
	if nx == 1
		xxcoll(1) = .6546536707;
	elseif nx==2
		xxcoll(1) = .5384693101;
		xxcoll(2) = .9061798459;
	elseif nx==3
		xxcoll(1) = .4058451514;
		xxcoll(2) = .7415311856;
		xxcoll(3) = .9491079123;
	elseif nx==4
		xxcoll(1) = .3242534234;
		xxcoll(2) = .6133714327;
		xxcoll(3) = .8360311073;
		xxcoll(4) = .9681602395;
	elseif nx==5
		xxcoll(1) = .2695431560;
		xxcoll(2) = .5190961292;
		xxcoll(3) = .7301520056;
		xxcoll(4) = .8870625998;
		xxcoll(5) = .9782286581;
	elseif nx==6
		xxcoll(1) = .2304583160;
		xxcoll(2) = .4484927510;
		xxcoll(3) = .6423493394;
		xxcoll(4) = .8015780907;
		xxcoll(5) = .9175983992;
		xxcoll(6) = .9841830547;
	end	
elseif (aaa==2)   %cylindrical
	if nx == 1
		xxcoll(1) = .5773502692;
	elseif nx==2
		xxcoll(1) = .4597008434;
		xxcoll(2) = .8880738340;
	elseif nx==3
		xxcoll(1) = .3357106870;
		xxcoll(2) = .7071067812;
		xxcoll(3) = .9419651451;
	elseif nx==4
		xxcoll(1) = .2634992300;
		xxcoll(2) = .5744645143;
		xxcoll(3) = .8185294874;
		xxcoll(4) = .9646596062;
	elseif nx==5
		xxcoll(1) = .2165873427;
		xxcoll(2) = .4803804169;
		xxcoll(3) = .7071067812;
		xxcoll(4) = .8770602346;
		xxcoll(5) = .9762632447;
	elseif nx==6
		xxcoll(1) = .1837532119;
		xxcoll(2) = .4115766111;
		xxcoll(3) = .6170011402;
		xxcoll(4) = .7869622564;
		xxcoll(5) = .9113751660;
		xxcoll(6) = .9829724091;
	end	
else   %planar
	if nx == 1
		xxcoll(1) = .4472135955;
	elseif nx==2
		xxcoll(1) = .3399810436;
		xxcoll(2) = .8611363116;
	elseif nx==3
		xxcoll(1) = .2386191861;
		xxcoll(2) = .6612093865;
		xxcoll(3) = .9324695142;
	elseif nx==4
		xxcoll(1) =.1834346425;
		xxcoll(2) =.5255324099; 
		xxcoll(3) =.7966664774;
		xxcoll(4) =.9602898565;
	elseif nx==5
		xxcoll(1) = .1488743390;
		xxcoll(2) = .4333953941;
		xxcoll(3) = .6794095683;
		xxcoll(4) = .8650633667;
		xxcoll(5) = .9739065285;
	elseif nx==6
		xxcoll(1) = .1252334085;
		xxcoll(2) = .3678314990;
		xxcoll(3) = .5873179543;
		xxcoll(4) = .7699026742;
		xxcoll(5) = .9041172564;
		xxcoll(6) = .9815606342;
	end
end

% form Q matrix
for i=1:n1
    z(i) = 1./(2.*i+aaa-2.);
end
for i=1:n1
	q(i,1) = 1.;
	for j=2:n1
		q(i,j) = xxcoll(i)^(2*j-2);
	end
end

% form c and d matrices
for j=1:n1
	ca = 2.*j-2;
	da = (2.*j-2)*(2*j+aaa-4);
	for i=1:n1
		c(i,j) = ca*xxcoll(i)^(2*j-3);
		d(i,j) = da*xxcoll(i)^(2*j-4);
	end
end

% calculate the inverse of q
qinvcoll = inv(q);

% form the aa, bb, and ww matrices
for i=1:n1
	wwcoll(i) = 0.0;
	for j=1:n1
		aacoll(i,j) = 0.;
		bbcoll(i,j) = 0.;
		for k=1:n1
			aacoll(i,j) = aacoll(i,j)+c(i,k)*qinvcoll(k,j);
			bbcoll(i,j) = bbcoll(i,j)+d(i,k)*qinvcoll(k,j);
		end
		wwcoll(i) = wwcoll(i) + z(j)*qinvcoll(j,i);
	end
end

clear c d

if kon == 0
fprintf('A matrix (aacoll) \n')
disp(aacoll)
fprintf('B matrix (bbcoll) \n')
disp(bbcoll) 
fprintf('Qinv matrix (qinvcoll) \n')
disp(qinvcoll) 
fprintf('Collocation points (xxcoll) \n')
disp(xxcoll)
fprintf('W matrix (wwcoll) \n')
disp(wwcoll)
end

		
			
% filename coll_interp.m
function y = coll_interp(ncol,dd,xb)
% ncol = number of internal collocation points, symmetric polynomials
% dd = the coefficients in the polynomial (vector)
% xb is a vector of the x values to evaluate the polynomial at

nn = size(xb);
for j=1:nn(2)
	x = xb(j);
	sum = 0;
	for i=1:ncol+1
		sum = sum*x*x + dd(ncol+2-i);
	end
	y(j) = sum;
end


%	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));