%University of Wales TSD
%Thomas Holland
%ACHA 5011 - Assignment 2 Task 2
%LOAD "RPM_VARIABLES" MATLAB DATA BEFORE RUNNING

readtable("Sigma 1.6L.xls")

Bore = 0.079;       %Bore (m)
Stroke = 0.0814;    %Stroke (m)
Rod = 0.137;        %Conrod Length (m)
CompR = 11.16;      %Compression Ratio

Crank = RPM_6500{:,1};          %Crank Angle (deg)
Volume = RPM_6500{:,2};          %Cylinder Volume (cc)


%Put Pressure for each RPM Range

Pressure_6500 = RPM_6500{:,3};  %Mean Cylinder Pressure @ 6500 RPM
Pressure_6000 = RPM_6000{:,3};  %Mean Cylinder Pressure @ 6000 RPM
Pressure_5500 = RPM_5500{:,3};  %Mean Cylinder Pressure @ 5500 RPM
Pressure_5000 = RPM_5000{:,3};  %Mean Cylinder Pressure @ 5000 RPM
Pressure_4500 = RPM_4500{:,3};  %Mean Cylinder Pressure @ 4500 RPM
Pressure_4250 = RPM_4250{:,3};  %Mean Cylinder Pressure @ 4250 RPM
Pressure_4000 = RPM_4000{:,3};  %Mean Cylinder Pressure @ 4000 RPM
Pressure_3500 = RPM_3500{:,3};  %Mean Cylinder Pressure @ 3500 RPM
Pressure_3000 = RPM_3000{:,3};  %Mean Cylinder Pressure @ 3000 RPM
Pressure_2500 = RPM_2500{:,3};  %Mean Cylinder Pressure @ 2500 RPM
Pressure_2000 = RPM_2000{:,3};  %Mean Cylinder Pressure @ 2000 RPM
Pressure_1750 = RPM_1750{:,3};  %Mean Cylinder Pressure @ 1750 RPM
Pressure_1500 = RPM_1500{:,3};  %Mean Cylinder Pressure @ 1500 RPM
Pressure_1000 = RPM_1000{:,3};  %Mean Cylinder Pressure @ 1000 RPM




%Build a Matrix to allow for easier implementation of all RPM Ranges into
%one Graph. Used a '...' to allow all of the array to fit into the script window

Pressure_Combined = [Pressure_1000 Pressure_1500 Pressure_1750 ...
    Pressure_2000 Pressure_2500 Pressure_3000 Pressure_3500 Pressure_4000 ...
    Pressure_4250 Pressure_4500 Pressure_5000 Pressure_5500 Pressure_6000 Pressure_6500];

Volume_Combined = [RPM_1000{:,2} RPM_1500{:,2} RPM_1500{:,2} RPM_2000{:,2} RPM_2500{:,2} ...
                   RPM_3000{:,2} RPM_3500{:,2} RPM_4000{:,2} RPM_4250{:,2} RPM_4500{:,2} ...
                   RPM_5000{:,2} RPM_5500{:,2} RPM_6000{:,2} RPM_6500{:,2}];

%Array for Each RPM Range to allow for a easier integration for the plots

RPM = [1000 1500 1750 2000 2500 3000 3500 4000 4250 4500 5000 5500 6000 6500];

%Plot Volume Against Crank Angle

[Crank_grid, RPM_grid] = meshgrid(RPM_6500{:,1}, RPM);

figure('Name','Crank Angle vs Volume vs RPM')

surf(Crank_grid, RPM_grid, Volume_Combined')

xlabel('Crank Angle (deg)')
ylabel('Engine Speed (RPM)')
zlabel('Volume (cc)')

title('3D Cylinder Volume vs Crank Angle and RPM')

colormap turbo
colorbar

grid on
view(45,30)

ylim([1000 6500])
xlim([-360 360])

shading interp

%This Graph shows the cylinder volume in cc in each degree of crank angle.
%The change in volume is not a linear process as the rate of change of
%volume slows at TDC (-360, 0 and 360 degrees) and at BDC ( -180 and 180
%degrees). This graph is more a control measure to make sure that the
%volume is constant at each of the different RPMs and that the volume is
%the least at TDC and the greatest at BDC.


%Plot Crank Angle against Pressure

[Crank_grid, RPM_grid] = meshgrid(RPM_6500{:,1}, RPM);

figure('Name','Crank Angle vs Pressure')

surf(Crank_grid, RPM_grid, Pressure_Combined')

xlabel('Crank Angle (deg)')
ylabel('Engine Speed (RPM)')
zlabel('Mean Cylinder Pressure (bar)')

title('3D Cylinder Pressure vs Crank Angle and RPM')

colormap turbo
colorbar
clim([0 80])

grid on
view(45,30)

ylim([1000 6500])
xlim([-360 360])
zlim([0 80])

shading interp
lighting phong
camlight




%The dataset for cylinder pressure at 1000 RPM appears to be incomplete
%upon analysis of the created graph. After further inspection the area of 
% incomplete data is between crank angle -29 and 71. There is also a
%potential that these values are unreadable and returning NaN to MATLAB 
%preventing them from being plot on the graph. There could be a fix 
%by editing the 'Sigma 1.6.xls' however this would be considered tampering
%with the data and could produce skewed or inaccurate results. Thankfully
%this doesn't prevent further analysis where the graph shows that peak
%pressure is at crank angle 20 degrees or within 3 degrees of rotation
%about that point which would maximise the force on the crankshaft on the
%expansion stroke. The graph also shows that when in higher bands of RPM
%there is an increase in pressure with the peak from this data being at
%4000 RPM with a pressure of 68.667 Bar, this is where peak torque from the 
%engine would be.


%Plot Pressure Against Volume

figure('Name','Pressure vs Volume')
plot(Volume,Pressure_Combined)
xlabel("Cylinder Volume (cc)")
ylabel("Mean Cylinder Pressure (Bar)")
title("Volume vs Mean Cylinder Pressure")
grid on
xlim([25 450])
ylim([0 80])

%This figure shows the pressure against volume at different stages in the
%four-stroke cycle. The line that runs along the bottom of the figure is
%the inlet stroke with the lower of the two curves being the compression
%stroke with where the base of the line that goes almost straight upwards
%being the start of the pressure increase caused by the expansion of buring
%fuel. The peak of the pressure being 69.035 Bar

%Calculate the Log of Volume and Pressures 
logV = log(Volume); 
logP_1000 = log(Pressure_1000); 
logP_1500 = log(Pressure_1500); 
logP_1750 = log(Pressure_1750); 
logP_2000 = log(Pressure_2000); 
logP_2500 = log(Pressure_2500); 
logP_3000 = log(Pressure_3000); 
logP_3500 = log(Pressure_3500); 
logP_4000 = log(Pressure_4000); 
logP_4250 = log(Pressure_4250); 
logP_4500 = log(Pressure_4500); 
logP_5000 = log(Pressure_5000); 
logP_5500 = log(Pressure_5500); 
logP_6000 = log(Pressure_6000); 
logP_6500 = log(Pressure_6500); 

%Calculate the Gradients using the Values at 6500RPM 
comp_index = Crank >= -180 & Crank <= -20; 
exp_index = Crank >= 20 & Crank <= 180; 

logP = log(Pressure_6500); 

% Compression fit 

p_comp = polyfit(logV(comp_index),logP_6500(comp_index),1); 

% Expansion fit 
p_exp = polyfit(logV(exp_index),logP_6500(exp_index),1); 

% Calculate fitted lines 

fit_comp = polyval(p_comp,logV(comp_index)); 
fit_exp = polyval(p_exp,logV(exp_index)); 

% Polytropic indices 

n_comp = -p_comp(1); n_exp = -p_exp(1); 

%Plot Log of Volume vs Log of Pressure 
figure('Name','Log of Volume vs Log of Pressure') 
hold on 
plot3(logV,1000*ones(size(logV)),logP_1000) 
plot3(logV,1500*ones(size(logV)),logP_1500) 
plot3(logV,1750*ones(size(logV)),logP_1750) 
plot3(logV,2000*ones(size(logV)),logP_2000)
plot3(logV,2500*ones(size(logV)),logP_2500) 
plot3(logV,3000*ones(size(logV)),logP_3000) 
plot3(logV,3500*ones(size(logV)),logP_3500) 
plot3(logV,4000*ones(size(logV)),logP_4000) 
plot3(logV,4250*ones(size(logV)),logP_4250)
plot3(logV,4500*ones(size(logV)),logP_4500)
plot3(logV,5000*ones(size(logV)),logP_5000)
plot3(logV,5500*ones(size(logV)),logP_5500)
plot3(logV,6000*ones(size(logV)),logP_6000)
plot3(logV,6500*ones(size(logV)),logP_6500) 

% Compression line (6500 RPM) represented by a Black Line Where 0 RPM would be 
plot3(logV(comp_index),0*ones(sum(comp_index),1),fit_comp,'k','LineWidth',3) 

% Expansion line (6500 RPM) represented by a Red Line Where 0 RPM would be 
plot3(logV(exp_index),0*ones(sum(exp_index),1),fit_exp,'r','LineWidth',3) 

xlabel('Log of Volume') 
ylabel('Engine Speed (RPM)') 
zlabel('Log of Pressure') 
title('Log of Pressure vs Log of Volume') 
grid on 
view(45,30) 
ylim([0 6500])
xlim([3.5 6.25]) 
zlim([-1.25 5]) 

%Display the calculated Poltropic Index in the Command Window
%Compression Polytropic Index (6500 RPM) = 1.311 
fprintf('Compression Polytropic Index (6500 RPM) = %.3f\n',n_comp) 
%Expansion Polytropic Index (6500 RPM) = 1.365 
fprintf('Expansion Polytropic Index (6500 RPM) = %.3f\n',n_exp)

%The above section of code creates a graph showing the log of pressure
%against the log of volume for 6500 RPM. There are also two trend
%lines where the gradients of which provide  the polytropic index
%for compression and expansion. The black trendline gives polytropic
%index for compression n_comp = 1.311 (which is also displayed in the
%command window). The red trendline gives the polytropic index for
%expansion wher n_exp = 1.365. The polytropic index is calculated over
%average of available RPMs however I have been unable to make my code work
%to be able to present these trendlines and values.

%Start of Rate of Change of Pressure Against Crank Angle

dP_dTheta_6500 = gradient(Pressure_6500, Crank);

%Plot Graph
figure('Name','Rate of Pressure Change at 6500 RPM')

plot(Crank, dP_dTheta_6500,'LineWidth',2)

xlabel('Crank Angle (deg)')
ylabel('dP/dθ (bar/deg)')

title('Rate of Change of Cylinder Pressure at 6500 RPM')

grid on
xlim([-180 180])

%The graph plotted above gives the rate of change of pressure against the
%Crank angle. This graph shows where combustion starts as well as its
%duration. From this graph ignition appears to start at an approximate
%crank angle of -15 with an ignition delay of approximately 4 degrees meaning
% combustion starts at crank angle -11 with the period of combustion ending 
% at crank angle 29 giving a duration of 40 degrees of crank angle. 
%These values can now be used to calculate the mass fraction burned using 
% the Wiebe function.

%Define Combustion Parameters
Crank = -180:1:180;     % Crank angle range

a = 5;                  % Wiebe constant
m = 2;                  % Shape factor
SOC = -11;              % Start of combustion (deg)
bdur = 40;              % Burn duration (deg)

%Calculate Mass Fraction Burned
for j = 1:length(RPM)

    for i = 1:length(Crank)

        if Crank(i) < SOC
            MFB_all(i,j) = 0;

        elseif Crank(i) <= SOC + bdur
            MFB_all(i,j) = 1 - exp(-a*((Crank(i)-SOC)/bdur)^(m+1));

        else
            MFB_all(i,j) = 1;
        end

    end

    % Burn rate (derivative)
    BurnRate_all(:,j) = gradient(MFB_all(:,j));

end

%Plot MFB vs Crank Angle

[Crank_grid, RPM_grid] = meshgrid(Crank, RPM);

figure('Name','Calculated MFB From Wiebe Function')

surf(Crank_grid, RPM_grid, MFB_all')

xlabel('Crank Angle (deg)')
ylabel('Engine Speed (RPM)')
zlabel('Mass Fraction Burned')

title('Calculated MFB From Wiebe Function')

colormap turbo
colorbar

view(45,30)
grid on
xlim([-40 80])
ylim([1000 6500])
zlim([0 1])

shading interp


%Plot Burn Rate vs Crank Angle

figure('Name','Burn Rate From Wiebe Function')

surf(Crank_grid, RPM_grid, BurnRate_all')

xlabel('Crank Angle (deg)')
ylabel('Engine Speed (RPM)')
zlabel('Burn Rate')

title('Burn Rate From Wiebe Function')

colormap turbo
colorbar

view(45,30)
grid on
xlim([-40 80])
ylim([1000 6500])

shading interp