%% Import Engine Data

clear
clc
close all

filename = 'Sigma 1.6L.xls';

% Read the 6500 rpm worksheet
opts = detectImportOptions(filename,'Sheet','6500');
opts.DataRange = 'A12';

T = readtable(filename,opts);

% Extract variables
theta = T{:,1};          % Crank angle (deg)
V = T{:,2};              % Cylinder volume (cc)
P_bar = T{:,3};          % Mean cylinder pressure (bar)

% Convert to SI units
V = V*1e-6;              % cc to m^3
P = P_bar*1e5;           % bar to Pa

disp('Engine data successfully imported.')
%% Display Basic Information

fprintf('Number of Data Points = %d\n',length(theta));

fprintf('Minimum Pressure = %.2f bar\n',min(P_bar));
fprintf('Maximum Pressure = %.2f bar\n',max(P_bar));

fprintf('Minimum Volume = %.2f cc\n',min(V)*1e6);
fprintf('Maximum Volume = %.2f cc\n',max(V)*1e6);
%% Pressure vs Crank Angle

figure

plot(theta,P_bar,'LineWidth',2)

grid on
box on

xlabel('Crank Angle (°)')
ylabel('Cylinder Pressure (bar)')
title('Cylinder Pressure Variation During One Engine Cycle')

xlim([-360 360])

set(gca,'FontSize',12)
%% Pressure-Volume Diagram

figure

plot(V*1e6,P_bar,'b','LineWidth',2)

grid on
box on

xlabel('Cylinder Volume (cm^3)')
ylabel('Cylinder Pressure (bar)')
title('Pressure–Volume Diagram')

set(gca,'FontSize',12)
%% Pressure Rise Rate

dP = gradient(P_bar)./gradient(theta);

figure

plot(theta,dP,'LineWidth',2)

grid on
box on

xlabel('Crank Angle (°)')
ylabel('dP/d\theta (bar/deg)')
title('Cylinder Pressure Rise Rate')

set(gca,'FontSize',12)

% Locate maximum pressure rise
[maxRise,index] = max(dP);

hold on
plot(theta(index),maxRise,'ro','MarkerFaceColor','r')

text(theta(index)+5,maxRise,...
    sprintf('Maximum = %.2f bar/deg',maxRise))
%% Compression Polytropic Index

compression = theta > -150 & theta < -30;

lnP = log(P(compression));
lnV = log(V(compression));

n_comp = -gradient(lnP)./gradient(lnV);

figure

plot(theta(compression),n_comp,'LineWidth',2)

grid on
box on

xlabel('Crank Angle (°)')
ylabel('Polytropic Index')

title('Compression Polytropic Index')

set(gca,'FontSize',12)

fprintf('Average Compression n = %.3f\n',...
    mean(n_comp,'omitnan'));
%% Expansion Polytropic Index

expansion = theta > 40 & theta < 160;

lnP = log(P(expansion));
lnV = log(V(expansion));

n_exp = -gradient(lnP)./gradient(lnV);

figure

plot(theta(expansion),n_exp,'LineWidth',2)

grid on
box on

xlabel('Crank Angle (°)')
ylabel('Polytropic Index')

title('Expansion Polytropic Index')

set(gca,'FontSize',12)

fprintf('Average Expansion n = %.3f\n',...
    mean(n_exp,'omitnan'));
%% Cylinder Gas Temperature

% Gas constant for air
R = 287;                 % J/kg.K

% Assume intake conditions
P0 = 101325;             % Pa
T0 = 300;                % K

% Estimate trapped mass at BDC
VBDC = max(V);

m = (P0*VBDC)/(R*T0);

% Calculate temperature using Ideal Gas Law
T = (P.*V)/(m*R);

figure

plot(theta,T,'LineWidth',2)

grid on
box on

xlabel('Crank Angle (°)')
ylabel('Temperature (K)')
title('Estimated Cylinder Gas Temperature')

set(gca,'FontSize',12)

fprintf('Maximum Temperature = %.0f K\n',max(T));
fprintf('Minimum Temperature = %.0f K\n',min(T));
%% Indicated Work

% Calculate work using numerical integration
W = trapz(V,P);

fprintf('\n');
fprintf('Net Indicated Work = %.2f J\n',W);

% Calculate IMEP
Vd = max(V)-min(V);

IMEP = W/Vd;

fprintf('Indicated Mean Effective Pressure = %.2f bar\n',IMEP/1e5);