%% This is the driver program to solve the simple plug flow
%% reactor example (nonisothermal).
%% Reaction
%%  
%% A = B     ... r1 = k1*cA    
%% B = C     ... r2 = k2*cA
%% B + C = D ... r3 = k3*cB*cC 
%% PUTENV('GNUTERM', 'dumb');
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                   Main program                    %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;clc;clf;
%%.... global (COMMON) blocks
global Rg
global nc nr ni pwe kr

%%.... Physical constants (single values)
Rg    = 8.314;      % J/mol/K, universal gas constant
%%-----------------------------------
%% .... Read imput form file
%%-----------------------------------
fid=fopen('genrk-input.dat','r');
%% Read comment  # nc nr j
trash=fgetl(fid);
[inpt,count] = fscanf(fid,'%i\n',Inf);
nc  = inpt(1); nr = inpt(2); 
key = inpt(3);

%% Read time
trash=fgetl(fid);
[inpt, count] = fscanf(fid,'%f\n',Inf);
tend = inpt(1);

%% Chem compounds (name)
trash = fgetl(fid);
idl = fgetl(fid);
i = 1;
while (any(idl))
  [temp,idl] = strtok(idl,' ');
  id{i} = temp;
  i = i + 1;
end

%% mol/m3/s, pre-exp. factor r(1)
trash = fgetl(fid);
kr    = fscanf(fid,'%f\n',Inf);

%% Inlet mol. conc.
trash = fgetl(fid);
c0    = fscanf(fid,'%f\n',Inf);

%% Read stech. coef. matrix
trash = fgetl(fid);
ni = zeros(nr,nc);
for i = 1 : nr
  tmp = fscanf(fid,'%f\n', [nc]);
  for j = 1 : nc
      ni(i,j) = tmp(j);
  end
end

%% Read matrix of pwr. law exponents
trash = fgetl(fid);
pwe = zeros(nr,nc);
for i = 1 : nr
  tmp = fscanf(fid, '%f\n', [nc]);
  for j = 1 : nc;
      pwe(i,j) = tmp(j,1);
  end
end
fclose(fid);

%%-----------------------------------
%%************ ODE SOLVER ***********
%%-----------------------------------
tspan =[0.0 tend]';
[t,c]=ode45(@dydz_genrk_M, tspan, c0);

%%.... Output
%fprintf('Conversion X(%s): %9.4f \n', id(key){:}, (y(key,key)-y(nsc,key))/y(key,key)); 
%fprintf('Reactor volume:   %9.3e m3 \n', Vr);
fprintf('Reaction system \n');
wrgs = zeros(nr);
prtrcell(nr, nc, ni, id, wrgs)
%%
%% .... plot the result
%%
plot(t,c)
legend (id)
title ('Molar concentration of components as function of time')
grid on
xlabel('t / s');
ylabel('c_{i} / (mol.m^{-3)}');
