%% TASK 2 - REAL CYLINDER-PRESSURE DATA ANALYSIS
% Ford Sigma 1.6 L spark-ignition engine


%% 1. USER SETTINGS
filename = 'Sigma_1.6L.xls';
sheetName = '6500';       % Engine-speed worksheet: 1000 to 6500 rpm
pressureColumn = 3;       % 3 = MEAN(B); 4:7 = individual cylinders 1:4
exportFigures = true;

% Thermodynamic assumptions used only for heat-release and temperature estimates
gamma = 1.25;             % ratio of specific heats
Rgas = 287;               % J/(kg K), air approximation
Tintake = 300;             % K, assumed trapped-charge temperature at -180 deg
smoothWindow = 5;         % samples used for derivative calculations

%% 2. IMPORT THE MEASURED DATA
% Workbook rows 13:732 contain 720 crank-angle samples.
% Columns: A crank angle, B volume, C mean pressure, D:G cylinders 1:4.
raw = readmatrix(filename, 'Sheet', sheetName, 'Range', 'A13:G732');
raw = raw(all(isfinite(raw(:,1:3)),2),:);

CA720 = raw(:,1);                        % crank angle, deg ATDC
V720 = raw(:,2)*1e-6;                    % cylinder volume, m^3
P720 = raw(:,pressureColumn);            % measured pressure, bar

% Extract the closed compression-combustion-expansion cycle.
closed = CA720 >= -180 & CA720 <= 180;
CA = CA720(closed);
V = V720(closed);
Pbar = P720(closed);

% Sort, remove duplicated angles and place data on a one-degree grid.
[CA, order] = sort(CA(:));
V = V(order); Pbar = Pbar(order);
[CA, uniqueIndex] = unique(CA, 'stable');
V = V(uniqueIndex); Pbar = Pbar(uniqueIndex);

CAgrid = (-180:1:180)';
V = interp1(CA, V, CAgrid, 'linear', 'extrap');
Pbar = interp1(CA, Pbar, CAgrid, 'linear', 'extrap');
CA = CAgrid;
Ppa = Pbar*1e5;

if any(Pbar <= 0)
    error('Pressure must be positive. Check the worksheet and pressureColumn setting.');
end

%% 3. ENGINE GEOMETRY AND PISTON POSITION
% Geometry is read from the workbook specification for the Sigma engine.
Bore = 0.079;             % m
Stroke = 0.0814;          % m
Rod = 0.137;              % m
CR = 11.16;
A = pi*Bore^2/4;
Vd = A*Stroke;
Vc = Vd/(CR-1);
r = Stroke/2;
theta = deg2rad(CA);

% Distance travelled downwards from TDC.
pistonDisplacement = r*(1-cos(theta)) + ...
    (Rod - sqrt(Rod^2 - (r*sin(theta)).^2));
pistonPositionMM = pistonDisplacement*1000;
Vgeometry = Vc + A*pistonDisplacement;
maxVolumeDifferenceCC = max(abs(Vgeometry-V))*1e6;

%% 4. PRESSURE AGAINST CRANK ANGLE AND PISTON POSITION
figure('Name','Pressure and piston position','Color','w');
yyaxis left
plot(CA, Pbar, 'LineWidth',1.6);
ylabel('Cylinder pressure (bar)');
yyaxis right
plot(CA, pistonPositionMM, '--', 'LineWidth',1.3);
ylabel('Piston displacement from TDC (mm)');
xline(0, ':', 'TDC');
xlabel('Crank angle (deg ATDC)');
title(sprintf('Cylinder pressure and piston position at %s rpm',sheetName));
legend('Measured pressure','Piston displacement','Location','best');
grid on; xlim([-180 180]);
if exportFigures, exportgraphics(gcf,'Task2_pressure_piston_position.png','Resolution',300); end

% Interpretation: pressure rises as the piston moves towards TDC and the
% cylinder volume falls. Combustion produces a rapid pressure rise around
% TDC, after which pressure falls as the piston moves down during expansion.

%% 5. PRESSURE-VOLUME DIAGRAM
figure('Name','PV diagram','Color','w');
plot(V*1e6, Pbar, 'LineWidth',1.6);
xlabel('Cylinder volume (cm^3)');
ylabel('Cylinder pressure (bar)');
title(sprintf('Measured pressure-volume diagram at %s rpm',sheetName));
grid on;
if exportFigures, exportgraphics(gcf,'Task2_PV_diagram.png','Resolution',300); end

% Interpretation: motion from large to small volume is the compression
% process. The high-pressure loop around minimum volume is produced by
% combustion. Expansion returns the piston towards maximum volume and the
% enclosed P-V area represents net indicated work from the closed cycle.

%% 6. POLYTROPIC INDICES
% For a polytropic process, P*V^n = constant, so:
% log(P) = constant - n*log(V).
% Windows avoid valve effects and the main combustion period.
compressionWindow = CA >= -140 & CA <= -35;
expansionWindow = CA >= 50 & CA <= 140;

fitComp = polyfit(log(V(compressionWindow)), log(Ppa(compressionWindow)), 1);
fitExp = polyfit(log(V(expansionWindow)), log(Ppa(expansionWindow)), 1);
nCompression = -fitComp(1);
nExpansion = -fitExp(1);

PcompFit = exp(polyval(fitComp,log(V(compressionWindow))))/1e5;
PexpFit = exp(polyval(fitExp,log(V(expansionWindow))))/1e5;

figure('Name','Polytropic fits','Color','w');
plot(V(compressionWindow)*1e6,Pbar(compressionWindow),'LineWidth',1.5); hold on;
plot(V(compressionWindow)*1e6,PcompFit,'--','LineWidth',1.5);
plot(V(expansionWindow)*1e6,Pbar(expansionWindow),'LineWidth',1.5);
plot(V(expansionWindow)*1e6,PexpFit,'--','LineWidth',1.5);
xlabel('Cylinder volume (cm^3)');
ylabel('Cylinder pressure (bar)');
title('Compression and expansion polytropic fits');
legend('Measured compression',sprintf('Fit, n = %.3f',nCompression), ...
       'Measured expansion',sprintf('Fit, n = %.3f',nExpansion), ...
       'Location','best');
grid on;
if exportFigures, exportgraphics(gcf,'Task2_polytropic_fits.png','Resolution',300); end

% Interpretation: n is the slope magnitude of the log(P)-log(V) relation.
% Values between approximately 1 and gamma indicate that heat transfer and
% changing gas properties make the real process different from an ideal
% isothermal or perfectly isentropic process.

%% 7. PRESSURE DERIVATIVE
PsmBar = movmean(Pbar,smoothWindow);
dPdCAbar = gradient(PsmBar,CA);          % bar/deg CA
[maxPressureRise,maxRiseIndex] = max(dPdCAbar);
maxRiseAngle = CA(maxRiseIndex);

figure('Name','Pressure derivative','Color','w');
plot(CA,dPdCAbar,'LineWidth',1.5); hold on;
xline(0,':','TDC');
xline(maxRiseAngle,'--',sprintf('Maximum rise: %.0f deg',maxRiseAngle));
yline(0,':');
xlabel('Crank angle (deg ATDC)');
ylabel('dP/d\theta (bar/deg CA)');
title('Measured cylinder-pressure derivative');
grid on; xlim([-80 120]);
if exportFigures, exportgraphics(gcf,'Task2_pressure_derivative.png','Resolution',300); end

% Interpretation: a strong positive dP/dtheta marks the rapid pressure rise
% during combustion. Negative values later in the cycle indicate pressure
% falling during expansion.

%% 8. APPARENT HEAT RELEASE AND COMBUSTION TIMING
% Single-zone, closed-system first-law approximation:
% dQ/dtheta = gamma/(gamma-1)*P*dV/dtheta
%            + 1/(gamma-1)*V*dP/dtheta.
% This is net apparent heat release because wall heat transfer, blow-by and
% crevice losses are not separately measured.
PsmPa = PsmBar*1e5;
dVdCA = gradient(V,CA);
dPdCAPa = gradient(PsmPa,CA);
dQnet = gamma/(gamma-1).*PsmPa.*dVdCA + ...
        1/(gamma-1).*V.*dPdCAPa;          % J/deg CA

combustionSearchWindow = CA >= -40 & CA <= 80;
dQpositive = zeros(size(dQnet));
dQpositive(combustionSearchWindow) = max(dQnet(combustionSearchWindow),0);
Qcumulative = cumtrapz(CA,dQpositive);
if max(Qcumulative) <= 0
    error('No positive apparent heat release detected. Check the data selection.');
end
xb = Qcumulative/max(Qcumulative);

CA05 = crossingAngle(CA,xb,0.05);
CA10 = crossingAngle(CA,xb,0.10);
CA50 = crossingAngle(CA,xb,0.50);
CA90 = crossingAngle(CA,xb,0.90);
CA95 = crossingAngle(CA,xb,0.95);
combustionDuration = CA90-CA10;

figure('Name','Apparent heat release','Color','w');
yyaxis left
plot(CA,dQnet,'LineWidth',1.3);
ylabel('Net apparent heat-release rate (J/deg CA)');
yyaxis right
plot(CA,xb,'LineWidth',1.5); hold on;
yline(0.10,':'); yline(0.50,':'); yline(0.90,':');
xline(CA10,'--','CA10'); xline(CA50,'--','CA50'); xline(CA90,'--','CA90');
ylabel('Cumulative apparent burned fraction');
xlabel('Crank angle (deg ATDC)');
title('Apparent heat release and combustion timing');
grid on; xlim([-50 100]);
if exportFigures, exportgraphics(gcf,'Task2_heat_release_combustion_timing.png','Resolution',300); end

% CA10 and CA90 are used as practical estimates of the start and end of the
% main combustion period. CA50 identifies the centre of heat release. They
% are more robust than defining combustion from a noisy derivative crossing.

%% 9. CYLINDER TEMPERATURE ESTIMATE
% Estimate trapped mass by applying the ideal-gas equation at -180 deg CA.
% The result is sensitive to the assumed intake temperature and residual gas.
referenceIndex = 1;
trappedMass = Ppa(referenceIndex)*V(referenceIndex)/(Rgas*Tintake);
Temperature = Ppa.*V/(trappedMass*Rgas);
[peakTemperature,peakTemperatureIndex] = max(Temperature);
peakTemperatureAngle = CA(peakTemperatureIndex);

figure('Name','Estimated temperature','Color','w');
plot(CA,Temperature,'LineWidth',1.5); hold on;
xline(0,':','TDC');
xline(peakTemperatureAngle,'--',sprintf('Peak: %.0f deg',peakTemperatureAngle));
xlabel('Crank angle (deg ATDC)');
ylabel('Estimated cylinder temperature (K)');
title('Estimated mean cylinder-gas temperature');
grid on; xlim([-180 180]);
if exportFigures, exportgraphics(gcf,'Task2_estimated_temperature.png','Resolution',300); end

% Interpretation: temperature rises during compression, increases sharply
% during combustion and then falls during expansion. It is a bulk-gas
% estimate, not a measured local flame temperature.

%% 10. WORK DONE DURING DIFFERENT STAGES
% Boundary work done by the gas is W = integral(P dV). Positive work is done
% by the gas; negative work is done on the gas.
preCombustion = CA >= -180 & CA <= CA10;
mainCombustion = CA >= CA10 & CA <= CA90;
postCombustion = CA >= CA90 & CA <= 180;
compressionStroke = CA >= -180 & CA <= 0;
expansionStroke = CA >= 0 & CA <= 180;

WpreCombustion = trapz(V(preCombustion),Ppa(preCombustion));
WmainCombustion = trapz(V(mainCombustion),Ppa(mainCombustion));
WpostCombustion = trapz(V(postCombustion),Ppa(postCombustion));
Wcompression = trapz(V(compressionStroke),Ppa(compressionStroke));
Wexpansion = trapz(V(expansionStroke),Ppa(expansionStroke));
WclosedCycle = trapz(V,Ppa);
IMEPbar = WclosedCycle/Vd/1e5;

workTable = table( ...
    [WpreCombustion; WmainCombustion; WpostCombustion; ...
     Wcompression; Wexpansion; WclosedCycle], ...
    'VariableNames',{'Work_J'}, ...
    'RowNames',{'PreCombustion','MainCombustion','PostCombustion', ...
                'CompressionStroke','ExpansionStroke','ClosedCycleNet'});

figure('Name','Stage work','Color','w');
bar(categorical(workTable.Properties.RowNames),workTable.Work_J);
yline(0,':');
ylabel('Boundary work (J)');
title('Work done by or on the cylinder gas');
grid on;
if exportFigures, exportgraphics(gcf,'Task2_stage_work.png','Resolution',300); end

% Negative compression work means work is supplied to the gas. Positive
% expansion work means the gas performs work on the piston. Their difference
% gives the net closed-cycle indicated work.

%% 11. SUMMARY OUTPUT
[peakPressure,peakPressureIndex] = max(Pbar);
peakPressureAngle = CA(peakPressureIndex);

summaryTable = table( ...
    nCompression,nExpansion,peakPressure,peakPressureAngle, ...
    maxPressureRise,maxRiseAngle,CA05,CA10,CA50,CA90,CA95, ...
    combustionDuration,peakTemperature,peakTemperatureAngle, ...
    Wcompression,Wexpansion,WclosedCycle,IMEPbar,trappedMass, ...
    'VariableNames',{ ...
    'nCompression','nExpansion','PeakPressure_bar','PeakPressureAngle_deg', ...
    'MaximumPressureRise_barPerDeg','MaximumPressureRiseAngle_deg', ...
    'CA05_deg','CA10_deg','CA50_deg','CA90_deg','CA95_deg', ...
    'CA10toCA90Duration_deg','PeakTemperature_K','PeakTemperatureAngle_deg', ...
    'CompressionWork_J','ExpansionWork_J','NetClosedCycleWork_J', ...
    'IMEP_bar','EstimatedTrappedMass_kg'});

disp(' ')
disp('TASK 2 RESULTS')
disp(summaryTable)
disp(workTable)

fprintf('\nEngine speed                         : %s rpm\n',sheetName);
fprintf('Peak pressure                       : %.2f bar at %.1f deg ATDC\n',peakPressure,peakPressureAngle);
fprintf('Compression polytropic index        : %.3f\n',nCompression);
fprintf('Expansion polytropic index          : %.3f\n',nExpansion);
fprintf('CA10 / CA50 / CA90                  : %.1f / %.1f / %.1f deg ATDC\n',CA10,CA50,CA90);
fprintf('Main combustion duration (10-90%%)   : %.1f deg CA\n',combustionDuration);
fprintf('Estimated peak temperature          : %.0f K at %.1f deg ATDC\n',peakTemperature,peakTemperatureAngle);
fprintf('Compression / expansion work        : %.2f / %.2f J\n',Wcompression,Wexpansion);
fprintf('Net closed-cycle work / IMEP         : %.2f J / %.2f bar\n',WclosedCycle,IMEPbar);
fprintf('Max geometry-volume difference       : %.4f cm^3\n',maxVolumeDifferenceCC);

writetable(summaryTable,'Task2_summary_results.csv');
writetable(workTable,'Task2_stage_work_results.csv','WriteRowNames',true);

%% LOCAL FUNCTION
function angle = crossingAngle(CA,y,target)
    idx = find(y >= target,1,'first');
    if isempty(idx)
        angle = NaN;
    elseif idx == 1
        angle = CA(1);
    else
        angle = interp1(y(idx-1:idx),CA(idx-1:idx),target,'linear');
    end
end