% MATLAB file used in the solution of Activity 13.1 of 
% "A First Course on Kinetics and Reaction Engineering."

% Ideal gas constant
R = 8.314e-3; % kJ/mol/K

% Data provided in the problem statement, in consistent units
V = 100; % L
% Enter the data as a big matrix, then break it out
dat = [305	20	1	1	0	0	307	0.3838
305	100	1	1	0	0	306	0.1473
305	300	1	1	0	0	306	0.0599
305	600	1	1	0	0	303	0.0315
305	600	2	5	0.5	0.5	304	0.7735
305	400	2	5	0.5	0.5	305	0.8772
305	200	2	5	0.5	0.5	306	1.1115
305	50	2	5	0.5	0.5	306	1.7045
305	5	2	5	0.5	0.5	304	2.3506
315	5	5	5	0.5	0.5	314	4.7473
315	40	5	5	0.5	0.5	314	3.6643
315	80	5	5	0.5	0.5	317	3.1304
315	120	5	5	0.5	0.5	314	2.7927
315	160	5	5	0.5	0.5	316	2.5465
315	160	0.5	0.25	0	0	316	0.0247
315	140	0.5	0.25	0	0	316	0.0257
315	100	0.5	0.25	0	0	315	0.0354
315	50	0.5	0.25	0	0	314	0.0605
315	20	0.5	0.25	0	0	315	0.1078
315	5	0.5	0.25	0	0	314	0.1759
325	1000	1	1	1	1	324	1.0615
325	1000	2	1	1	1	326	1.1118
325	1000	1	2	1	1	326	1.1126
325	1000	1	1	2	1	325	2.059
325	1000	1	1	1	2	326	1.0613
325	5	2	1	0.2	0.5	327	1.1362
325	50	2	1	0.2	0.5	324	0.8449
325	200	2	1	0.2	0.5	326	0.558
325	600	2	1	0.2	0.5	324	0.3711
325	600	0.5	3	0.5	0.1	326	0.6237
325	300	0.5	3	0.5	0.1	326	0.691
325	20	0.5	3	0.5	0.1	326	0.9494
340	5	4	2	1	1	341	2.9659
340	25	4	2	1	1	340	2.8588
340	100	4	2	1	1	341	2.5767
340	200	4	2	1	1	340	2.3394
340	200	2	4	2	0	338	3.3409
340	150	2	4	2	0	341	3.4456
340	100	2	4	2	0	340	3.5761
340	75	2	4	2	0	340	3.6547
340	25	2	4	2	0	339	3.8572
340	5	2	4	2	0	342	3.9668
360	5	3	1	0.1	1	360	1.0944
360	50	3	1	0.1	1	362	1.0427
360	100	3	1	0.1	1	361	0.9946
360	200	3	1	0.1	1	360	0.9156
360	200	3	5	1	0.1	359	3.5076
360	100	3	5	1	0.1	361	3.7089
360	50	3	5	1	0.1	358	3.8373
360	25	3	5	1	0.1	360	3.9166
360	5	3	5	1	0.1	359	3.9809];
T0 = dat(:,1);
VFR = dat(:,2);
CA0 = dat(:,3);
CB0 = dat(:,4);
CY0 = dat(:,5);
CZ0 = dat(:,6);
T = dat(:,7);
CY = dat(:,8);

% Calcuate x and y
nY = VFR.*CY;
nY0 = VFR.*CY0;
y = nY - nY0;
nA = (VFR.*CA0) - nY + nY0;
nB = (VFR.*CB0) - nY + nY0;
CA = nA./VFR;
CB = nB./VFR;
xx = V*(CA.*CB);

% allocate storage for the the fitting results of part 1 for 5 data blocks
T_block = zeros(5,1);
corr = zeros(5,1);
k = zeros(5,1);
k_u = zeros(5,1);

% step 1, process blocks separately
% 305 K block, data points 1 through 9
T_block(1) = 305;
x = xx(1:9);
y_hat = y(1:9);
figure;
FitLinmSR;
corr(1) = r_squared;
k(1) = m;
k_u(1) = m_u;

% 315 K block, data points 10 through 20
T_block(2) = 315;
x = xx(10:20);
y_hat = y(10:20);
figure;
FitLinmSR;
corr(2) = r_squared;
k(2) = m;
k_u(2) = m_u;

% 325 K block, data points 21 through 32
T_block(3) = 325;
x = xx(21:32);
y_hat = y(21:32);
figure;
FitLinmSR;
corr(3) = r_squared;
k(3) = m;
k_u(3) = m_u;

% 340 K block, data points 33 through 42
T_block(4) = 340;
x = xx(33:42);
y_hat = y(33:42);
figure;
FitLinmSR;
corr(4) = r_squared;
k(4) = m;
k_u(4) = m_u;

% 360 K block, data points 43 through 51
T_block(5) = 360;
x = xx(43:51);
y_hat = y(43:51);
figure;
FitLinmSR;
corr(5) = r_squared;
k(5) = m;
k_u(5) = m_u;

% Show the results of step 1
disp('Step 1');
disp('T_block r_squared    k    k_confid');
res = [T_block corr k k_u]

% Show all 5 model plots on one graph
y_model = zeros(51,1);
y_model(1:9) = k(1)*xx(1:9);
y_model(10:20) = k(2)*xx(10:20);
y_model(21:32) = k(3)*xx(21:32);
y_model(33:42) = k(4)*xx(33:42);
y_model(43:51) = k(5)*xx(43:51);
figure;
plot(xx,y,'o',xx,y_model,'-');
xlabel('x');
ylabel('y');
legend('data','model');
title('Combined Model Plots');

% Step 2, fit the Arrhenius expression to the k vs. T data
x = 1./T_block;
y_hat = log(k);
figure;
FitLinmbSR;
xlabel('1/T');
ylabel('ln(k)');
title('Arrhenius Plot');

% Calculate k0 and E and their uncertainties
E = -m*R
E_u = R*m_u
k0 = exp(b)
k0_u = b_u*exp(b)