%% PART A

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

%Clears data
clear;

%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 = 2720 + 260*cos((2*pi*0)/(tau));

%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.
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/Vmax)*(Ka)*(3*c(6) - Y1star(c) - Y2star(c) - Y3star(c))];

%Sets the proper tolerance (very tight tolerance in this case) to reduce error
x = 0.1;
y = 0.099999999;

error = 100;

tol = 1;

while error > tol

    options = odeset('RelTol',x);
    options2 = odeset('RelTol',y);

    %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 0.5] ,[.145 .196 .164 .181 .195 .135],options);
    [kareem,jabbar] = ode45 (dxdt, [0 0.5] ,[.145 .196 .164 .181 .195 .135],options2);


    [base,stupidval1] = size(baller);
    [ketball,stupidval2] = size(jared);
    [santana,stupidval1] = size(kareem);
    [johan,stupidval2] = size(jabbar);

    %Formats the output LONG
    format long;

    %Extracts the output matrices from Part A into new vectors
    for i=1:ketball
        gunit(i) = jared(i,1);
    end

    for i=1:base
        easye(i) = baller(i,6);
    end

    for i=1:johan
        jj(i) = kareem(i,1);
    end

    for i=1:santana
        putz(i) = jabbar(i,6);
    end


%Computes a trapezoid rule on the new vectors to find the average YO2
%concentration
abridgetoofar = (trapz(gunit,easye))/(0.5);
abridgetoonear = (trapz(jj,putz))/(0.5);

error = abs(((abridgetoonear - abridgetoofar) / (abridgetoonear))*100);

x = x - 0.00001;
y = y - 0.00001;

end

format long;

disp(['The desired accuracy is: (%) ',num2str(tol)]);
disp(['The associated value is: ',num2str(abridgetoonear)]);
disp(['The associated previous value is: ',num2str(abridgetoofar)]);
disp(['The percent relative error is: (%) ',num2str(error)]);
disp(['The tolerance associated with this accuracy is: ',num2str(y)]);
disp(['The previous tolerance associated with this accuracy is: ',num2str(x)]);
