% rxnOCFE.m
% This program runs a code to solve the reaction-diffusion
% problem using the method of orthogonal collocation on finite elements.

global xoc woc aoc boc qinv np ne nt delx Aocfe Bocfe

% set linear to 1 if problem is linear - this avoids one unnecessary iteration
linear = 0
% set the type of reaction rate
iwhich = 3
%set Thiele modulus and geometry (geometry = 1, 2, 3 for planar, cylindrical, spherical)
phi2 = phi^2
geometry = 3

% np is the total number of collocation points in one finite element (z direction)
np = 4
% ne is the number of elements in the z-direction
ne = 10
% This can be changed to produce a variable grid
for k=1:ne
	delx(k) = 1/ne;
end
delx(1)

% nt is the total number of points
nt = (np-1)*ne + 1
% number of interior points is np-2 per element
nint=(np-2)*ne;

% create the collocation matrices for the fluid
planar(np,0)

% set up the z values for plotting purposes
num=0;
xstart = 0.;
for k=1:ne
	for j=1:np-1
		num=num+1;
		x(num)=xstart+delx(k)*xoc(j);
	end
xstart = xstart+delx(k);
end
x(nt) = 1.

% iterate
tol = 1.e-10;
change = 1.;
iter=0;
% initial guess
%for i=1:nt
%	c(i) = 1;
%end
%c(nt) = 1.;

while change>tol
	iter=iter+1;

% set the matrices to zero	
		Bocfe = zeros(nt,1);
		Aocfe = zeros(np,np,ne);
		
% set the matrices for the differential equation
	for k=1:ne
		for j=2:np-1
			index = (np-1)*(k-1) + j;
			term = (geometry-1)/x(index);
			for i=1:np
				Aocfe(j,i,k) = boc(j,i) + term*delx(k)*aoc(j,i);
			end
			ans = rate1(c(index),iwhich);
			raterxn = ans(1);
			drate = ans(2);
			term=phi2*delx(k)*delx(k);
			Aocfe(j,j,k) = Aocfe(j,j,k) - term*drate;
			Bocfe(index) = term*(raterxn - drate*c(index));
		end		
	end

% set the continuity conditions
	for k=1:ne-1
		for i=1:np
			Aocfe(np,i,k) = aoc(np,i)/delx(k);
		end
	end
	for k=2:ne
		for i=1:np
			Aocfe(1,i,k) = -aoc(1,i)/delx(k);
		end
	end
	for k=2:ne
		term = Aocfe(np,np,k-1)+Aocfe(1,1,k);
		Aocfe(np,np,k-1) = term;
		Aocfe(1,1,k) = term; % this is redundant,since it isn't used
		Bocfe((np-1)*(k-1)+1) = 0.;
	end
	for i=1:np
		Aocfe(1,i,1) = aoc(1,i);
	end
	Aocfe(np,np,ne) = 1.;
	Bocfe(nt) = 1;

% Solve one iteration.
	OCFElud(np,ne);
	OCFEfas(np,ne,nt);
	c11 = Bocfe;
	
% Calculate the criterion to stop the iterations.
	sum=0.;
	for i=1:nt
		sum = sum + abs(c11(i) - c(i));
		c(i) = c11(i);
	end
	if linear==1,break,end
	change = sum
end

% Plot the solution.
plot(x,c,'-r*')
% Find the coefficients
dd1 = qinv*c';
xp = 0:0.02:1; %plotting points
cinterp = planar_interp(nx,dd1,xp,0,1);
plot(xoc,c,'*b',xp,cinterp,'-b')
xlabel('z')
ylabel('c')
c

sum = 0;
for k=1:ne
	for i=1:np
		index = (np-1)*(k-1) + i;
		ans = rate1(c(index),iwhich);
		raterxn = ans(1);
		sum = sum + delx(k)*woc(i)*raterxn*x(index)^(geometry-1);
	end
end
efffectiveness = geometry*sum



function planar(nx,kon)                                                
%   this subroutine computes the matrices for collocation without       
%   symmetry, using collocation points in table 4-3.                    
%   input variable                                                     
%     nx = number of collocation points, including two end points                     
%                                                                       
%   output variables, must be declared global                                                    
%     aoc = matrix for first derivative, eq.  4-103                       
%     boc = matrix for second derivative, eq. 4-103                       
%     qinv = matrix for q inverse, eq. 4-101                               
%     xoc = vector of collocations points, from table 4-3                 
%     woc = vector of weights, eq. 4-106                                  
% 
global xoc woc aoc boc qinv                                                                      
if(nx<3) 
  nx=3;
  fprintf('\n ***************************************** \n')
  fprintf('The number of collocation points must be between \n')
  fprintf('3 and 7 (inclusive).  nx has been changed to %g\n\n',nx)
end
if(nx>7)
  nx=7;
  fprintf('\n ***************************************** \n')
  fprintf('The number of collocation points must be between \n')
  fprintf('3 and 7 (inclusive).  nx has been changed to %g\n\n',nx)
end
% now 3<=nx<=7
if kon==0 
	fprintf('The total number of collocation points is %g\n\n',nx)
end

xoc(1)=0.;
if (nx==3)
  xoc(2)=0.5;
  elseif (nx==4)
  xoc(2)=0.21132486540519;
  xoc(3)=1-xoc(2);
  elseif (nx==5)
  xoc(2)=0.11270166537926;
  xoc(3)=0.5;
  xoc(4)=1-xoc(2);
  elseif (nx==6)
  xoc(2)=0.06943184420297;
  xoc(3)=0.33000947820757;
  xoc(4)=1-xoc(3);
  xoc(5)=1-xoc(2);
  else
  xoc(2)=0.04691007703067;
  xoc(3)=0.23076534494716;
  xoc(4)=0.5;
  xoc(5)=1-xoc(3);
  xoc(6)=1-xoc(2);
end
xoc(nx)=1.;
format long                                                         
for i=1:nx,                                                        
    r(i,i) = 0.0;                                                     
    aoc(i,i) = 0.0;                                                     
    s(i) = 1.0;                                                       
    boc(i,i) = 0.0 ;                                                    
     for j=1:nx,                                                     
        if (i~=j)                                         
           r(i,j) = 1.0/(xoc(i)-xoc(j)) ;                                     
           s(i) = s(i)*r(i,j);
        end                                            
     end                                                         
     for j=1:nx,                                                     
         jx = nx-j+1;                                                   
         if (jx<j), break, end                                         
         if (jx==j) 
            aoc(i,i) = aoc(i,i)+r(i,j);
         end                           
         if (jx>j) 
            aoc(i,i) = aoc(i,i)+r(i,j)+r(i,jx);
         end                   
     end                                                        
end                                                            
for i=1:nx                                                        
   for j=1:nx                                                     
       if (i~=j)                                           
          aoc(i,j) = s(j)*r(i,j)/s(i);                                    
          boc(i,j) = 2.0*aoc(i,j)*(aoc(i,i)-r(i,j));                           
          boc(i,i) = boc(i,i)+r(i,j)*(aoc(i,i)-r(i,j)); 
       end                       
   end                                                        
end                                                           
for i=1:nx                                                        
   qinv(1,i) = s(i);                                                    
   k = 1;                                                            
   woc(i) = 0.0 ;                                                      
   for j=1:nx                                                     
      if (j~=i)                                           
         l = k;                                                         
         k = k+1;                                                       
         qinv(k,i) = qinv(l,i);                                               
         while (l>1)                                          
            m = l-1;                                                       
            qinv(l,i) = qinv(m,i)-xoc(j)*qinv(l,i);                                  
            l = m;                                                        
         end                                                                                                                  
           qinv(1,i) = -xoc(j)*qinv(1,i);
      end                                       
   end                                                         
   for j=1:nx                                                     
      woc(i) = woc(i)+qinv(j,i)/j;
   end                                     
end

if kon == 0
fprintf('A matrix (aoc) \n')
disp(aoc)
fprintf('B matrix (boc) \n')
disp(boc) 
fprintf('Q matrix (qinv) \n')
disp(qinv) 
fprintf('Collocation points (xoc) \n')
disp(xoc)
fprintf('W matrix (woc) \n')
disp(woc)
end



% filename planar_interp.m
function y = planar_interp(nx,dd,x,x0,x1)
% nx = number of collocation points, including end points
% dd = the coefficients in the polynomial (vector)
% x is a vector of the x values to evaluate the polynomial at
% The points start at x0 and go to x1, and the xb below has to 
% be scaled to go from 0 to 1.

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





% This subroutine does the LU decomposition for the
% orthogonal collocation method on finite elements.
% This version is for a single unknown, giving a block
% diagonal matrix with a single entry of overlap.
% This subroutine is in the Appendix
% of "Nonlinear Analysis in Chemical Engineering".
% It assumes diagonal pivoting is OK.
% The matrix Aocfe is stored as Aocfe(np,np,ne), where
% np is the number of collocation points in one element
% and ne is the number of elements.

function OCFElud(np,ne)
global Aocfe

n1 = np - 1;
for l=1:ne
	for k=1:n1
		k1=k+1;
		for i=k1:np
			s = Aocfe(i,k,l)/Aocfe(k,k,l);
			Aocfe(i,k,l) = s;
			for j=k1:np
				Aocfe(i,j,l) = Aocfe(i,j,l) - s*Aocfe(k,j,l);
			end
		end
	end
	if l<ne
		Aocfe(1,1,l+1) = Aocfe(np,np,l);
	end
end


% This subroutine does the fore and aft sweep for the
% orthogonal collocation method on finite elements.
% This version is for a single unknown, giving a block
% diagonal matrix with a single entry of overlap.
% This subroutine is in the Appendix
% of "Nonlinear Analysis in Chemical Engineering".
% It assumes diagonal pivoting is OK.
% The matrix Aocfe is stored as Aocfe(np,np,ne), where
% np is the number of collocation points in one element
% ne is the number of elements.
% nt is the total number of unknowns, nt = (np-1)*ne + 1
% To use OCFEfas, you must first run OCFElud.
% The right-hand side is in Bocfe(nt)

function OCFEfas(np,ne,nt)
global Aocfe Bocfe

np1=np;
% forward sweep
for l=1:ne
	for i=2:np1
		i2=(l-1)*(np-1)+i;
		s = 0.;
		i1=i-1;
		for j=1:i1
			j2=i2-i+j;
			s=s+Aocfe(i,j,l)*Bocfe(j2);
		end
		Bocfe(i2)=Bocfe(i2)-s;
	end
end
% back substitution
for l1=1:ne
	l=ne-l1+1;
	if (l==ne)
		Bocfe(nt)=Bocfe(nt)/Aocfe(np,np,ne);
	end
	n1 = np-1;
	for k=1:n1
		i=n1+1-k;
		i2=(l-1)*(np-1)+i;
		m=i+1;
		n2=n1+1;
		s=0.;
		for j=m:n2
			j2=(l-1)*(np-1)+j;
			s=s+Aocfe(i,j,l)*Bocfe(j2);
		end
		Bocfe(i2)=(Bocfe(i2) - s)/Aocfe(i,i,l);
	end
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));