% listing of main file "kin2.m"
% this file calls the user function file "kin2F.m"
% and uses the standard matlab function lsqcurvefit (least square curve fit)

% The input data are the weight of catalyst (W), volumetric flow
% rate (v), input concentrations (cAin cBin) and concentrations
% in the reactor (cA cB) for a CSTR.
% These data are fit to a rate equation that is specified in the
% file "kin2F.m" that contains three rate parameters whose
% final values are returned in the array x below.
% In the file "kin2F.m", the rate equation specified is:
%     rcalc = (kf * cA - kr * cB) ./ (1 + KA * cA + KB * cB)

% run separator
fprintf('------------------------------------------\n') 

% input data
W = 100; % 100 g of catalyst
v = 10;  % 10 liter/hr flow rate
cAin = [0.559;1.54;1.13;0.755;0.164;0.326;1.36;0.515;0.596;0.907];
cBin = [0.761;0.035;0.435;0.185;0.186;0.794;0.0488;0.385;0.724;0.0633];
cA = [0.42;0.72;0.62;0.37;0.11;0.32;0.64;0.31;0.43;0.41];
cB = [0.9;0.85;0.95;0.57;0.24;0.8;0.77;0.59;0.89;0.56];

% calculate "experimental" rates = -rA(mol/g hr)
% dot . NOT needed before the * here since (v/W) is scalar

rexp = (v/W) * (cAin - cA);

% put data in form for input to matlab function lsqcurvefit and user function kin2F
xdata = [cA cB];
ydata = rexp;

x0 = [1;1;1;1]; % initial guesses of kf, kr, KA, KB

[x, resnorm] = lsqcurvefit('kin2F',x0,xdata,ydata)

fprintf('row 1 of x is kf value, row 2 is kr, row 3 is KA, row 4 is KB \n')
fprintf('and "resnorm" is sum of squared errors (squared 2-norm of the residual at x) \n\n')

fprintf('SEE PLOT \n\n')

% plot

rcalc = kin2F(x,xdata);

% get indices of max and min values of rexp
% to use in drawing diagonal line
% (or just plot rexp vs. rexp)
mini = find(rexp == min(rexp)); % indices of min value
maxi = find(rexp == max(rexp)); % indices of max value
rdiag = [rexp(mini),rexp(maxi)];

plot(rexp,rcalc,'bo',rdiag,rdiag,'k') 
title('kin2 problem')
ylabel('rcalc (mol/g/hr)')
xlabel('rexp (mol/g/hr)')