%% PART C

%Concept of c:
%C = [XV XA XP1 XP2 XP3 YO2];

%Clears data
clear;

%Holds the graphs
hold all;

%Initial conditions
Ka = 3120;
M = 250;
Qa = 4200;
Qb = 5000;
Va = 1500;
Vp = 20;
Vv = 40000;
sigma = .99;
Yl = .198;
%Vl = @(t) 2720 + 260*cos((2*pi*t)/(tau));
%dVl = @(t) -104*pi*sin((2*pi*t)/(tau));

% Y1star = @(c) (0.00032*c(2))/(0.20155 - c(2)) + (11.27 - c(2))*(0.055*c(2));
% Y2star = @(c) (0.00032*c(3))/(0.20155 - c(3)) + (11.27 - c(3))*(0.055*c(3));
% Y3star = @(c) (0.00032*c(4))/(0.20155 - c(4)) + (11.27 - c(4))*(0.055*c(4));
%Vmax = 2960;

%Creats a matrix of dxdt to be solved by ode45, sets equations as part of
%the matrix c with distinct elements of XC, XA, XP1, XP2, XP3, and YO2.
%ode45 will solve the matrix into the distinct elements as matrices of
%points. This time a new function is put in from number 2 of part A to
%solve for a new concentration
dxdt = @(t, c) [(1/Vv)*(c(2)*Qb-c(1)*Qb-M);
                (1/Va)*(sigma*Qb*c(5)+(1-sigma)*Qb*c(1)-c(2)*Qb);
                (1/Vp)*(sigma*Qb*c(1)-sigma*Qb*c(3)+Ka*(c(6)-(((.00032*c(3))/(.20155 - c(3))) + (11.27 - c(3))*(.055*c(3)))));
                (1/Vp)*(sigma*Qb*c(3)-sigma*Qb*c(4)+Ka*(c(6)-(((.00032*c(4))/(.20155 - c(4))) + (11.27 - c(4))*(.055*c(4)))));
                (1/Vp)*(sigma*Qb*c(4)-sigma*Qb*c(5)+Ka*(c(6)-(((.00032*c(5))/(.20155 - c(5))) + (11.27 - c(5))*(.055*c(5)))));
                %(1/(2720+260*cos(24*pi*t)))*(Qa*Yl-Qa*c(6)-(-6240*pi*sin(24*pi*t))*c(6)-3120*(3*c(6)-(((.00032*c(3))/(.20155 - c(3))) + (11.27 - c(3))*(.055*c(3)))-(((.00032*c(4))/(.20155 - c(4))) + (11.27 - c(4))*(.055*c(4)))-(((.00032*c(5))/(.20155 - c(5))) + (11.27 - c(5))*(.055*c(5)))))];
                -(1/(2720+260*cos(24*pi*0)))*((Ka)*(3*c(6) - (((.00032*c(3))/(.20155 - c(3))) + (11.27 - c(3))*(.055*c(3)))-(((.00032*c(4))/(.20155 - c(4))) + (11.27 - c(4))*(.055*c(4)))-(((.00032*c(5))/(.20155 - c(5))) + (11.27 - c(5))*(.055*c(5)))))];
                
%Sets the proper tolerance (very tight tolerance in this case) to reduce error                
options = odeset('RelTol',1e-6);

%Solves matrix dxdt with respect to t and c(six variables) and puts the
%values into a matrix of t coordinates and each varible points.
%Inputs initial conditions as well
[jared,baller] = ode45 (dxdt, [0 5] ,[.145 .196 .164 .181 .195 .135],options);

%Concentration to solve for on the graph, found by zooming in, not
%analytically
NewKa = 0.7*0.196;

%Xv
plot(jared,baller(:,1));
%Xa
plot(jared,baller(:,2));
%XP1
plot(jared,baller(:,3));
%XP2
plot(jared,baller(:,4));
%XP3
plot(jared,baller(:,5));
%YO2
plot(jared,baller(:,6));

vandame = 2.9174;
plot(vandame,NewKa,'rs');

%Adds a legend to the plot
legend('XV','XA','XP1','XP2','XP3','YO2');
xlabel('Time (Minutes)');
ylabel('Concentration');

%% Part C
%The concentration occurs at t = 2.9174, solved by looking at the graph
%with a very tight tol.


