clc;clear;

%% Task 3(A)
% Details:
% This script is for task 3(A), consisting of the original "Task_3" script 
% and function "mydissfuniii" with a tic-toc and iteration counter added. 

% The bisection method (my script) has been added after the orignial
% script, containing the following functions "nEquations",
% "BisectionErrSolve" and "GibbsEngCalcs". Again a tic-toc and iteration
% counter have been added. With Comparisons of the two methods at the end.

% IMPORTANT:
% This file requires matlab "gibbseng.mat" file and following functions to
% work: "GibbsEngCalcs", "BisectionErrSolve", "nEquations", "mydissfuniii".  

%% Original Task 3 Script

tic
c=8;
h=18;
o=0;
airs=12.5;
temp=(100:25:3000)';
press=(20:5:120)';
feq=(0.8:0.02:1.3)';
n5ntDis=zeros(49,4,26);
n5Dis=zeros(49,4,26);
e=zeros(49,4,26);
n1Dis=zeros(49,4,26);
n2Dis=zeros(49,4,26);
n3Dis=zeros(49,4,26);
n4Dis=zeros(49,4,26);
n6Dis=zeros(49,4,26);
ItterationsDiss=zeros(49,4,26);


for i=1:49
    for j=1:21
        for b=1:26
            [n5ntDis(i,j,b),n5Dis(i,j,b),e(i,j,b),n1Dis(i,j,b),n2Dis(i,j,b),n3Dis(i,j,b),n4Dis(i,j,b),n6Dis(i,j,b),ItterationsDiss(i,j,b)]=mydissfuniii(c,h,o,temp(i,1),feq(b,1),airs,press(j,1));
        end
    end
end
tendDis = toc;
fprintf('Time Elapsed Originial Script = %.6f\n',tendDis)


%% Bisection Method (my script and functions)

% Define variables
tic; load gibbseng.mat;
pressure = (20:5:120)';         % Define pressure variables
temperature = (100:25:3000)';   % Define temperature variables
fuelequiv = (0.8:0.02:1.3)';    % Define Fuel equivalance ratio variables
carbon=8;
hydrogen=18;
oxygen=0;
molar_air_stoich=12.5;
n5ntvalue=zeros(49,4,26);

% Populating arrays with zeros 
n1=zeros(49,4,26); n2=zeros(49,4,26); n3=zeros(49,4,26);
n4=zeros(49,4,26); n5=zeros(49,4,26); n6=zeros(49,4,26);
Error=zeros(49,4,26); Itterations=zeros(49,4,26);

% Bisection method variables
a = 0;                   % Lower bound (x1)
b = 2;                   % Upper bound (x3)
tolerance = 1e-12;       % Error tolerance.
max_iterations = 100;    % Maximum number of iterations

% For loops for pressure (j), temperature(i) and fuelequiv (q)
for i=1:49
    for j=1:21
        for q=1:26
[Error(i,j,q),n5ntvalue(i,j,q),n1(i,j,q),n2(i,j,q),n3(i,j,q),n4(i,j,q),n5(i,j,q),n6(i,j,q),Itterations(i,j,q)] = BisectionErrSolve(a,b,max_iterations...
 ,tolerance,pressure(j,1),temperature(i,1),fuelequiv(q,1),carbon,oxygen,hydrogen,molar_air_stoich,gibbseng);
        end
    end
end
tend = toc;
fprintf('Time Elapsed Bisection Method = %.6f\n',tend)


%% Comparing Different Methods

Method = ["Original";"Bisection"];
Time_Taken = [tendDis;tend];

% Total cumulative sum of each methods iterations
ItterationsCumsumDiss = cumsum(ItterationsDiss);
ItterationsCumsum = cumsum(Itterations);
Total_Iters = [ItterationsCumsumDiss(end);ItterationsCumsum(end)];

% Total cumulative sum of each methods error
ErrorCumsumDiss = cumsum(e);
ErrorCumsum = cumsum(Error);
Total_Error = [ErrorCumsumDiss(end);ErrorCumsum(end)];

% Total cumulative sum of difference in method errors
Error_Diff = abs(e-Error);
ErrorDiffCumsumDiss = cumsum(Error_Diff);
Total_Error_Diff = [ErrorDiffCumsumDiss(end);ErrorDiffCumsumDiss(end)];

Comparing_Methods = table(Method,Time_Taken,Total_Iters,Total_Error,...
    Total_Error_Diff);
display(Comparing_Methods)