clc; clear; close all;

%% === Select Excel File ===
[filename, pathname] = uigetfile({'*.xlsx;*.xls','Excel Files (*.xlsx, *.xls)'}, ...
                                 'Select the Excel file');
if isequal(filename,0)
    error('No file selected. Exiting.');
end
fullFile = fullfile(pathname, filename);

%% === Read Data (skip header row) ===
data = readmatrix(fullFile, 'Range', 'A2:H10000');

%% === Extract Columns ===
crank   = data(:,1);        % Crank Angle (deg)
vol     = data(:,2);        % Cylinder Volume (m^3)
rpmCols = 3:8;              % Pressure columns (Pa)
rpmVals = [1000 2000 3000 4000 5000 6000];
nCols   = length(rpmCols);

%% === Convert Pressure to bar ===
P_bar = data(:,rpmCols) / 1e5;     % bar
P_Pa  = data(:,rpmCols);           % keep Pa for temperature calc

%% === Pressure vs Crank Angle ===
figure(1)
subplot(2,2,1)
plot(crank, P_bar, 'LineWidth',1.2)
title('Pressure vs Crank Angle')
xlabel('Crank Angle (°)')
ylabel('Pressure (bar)')
legend(string(rpmVals)+" rpm",'Location','best')

%% === Pressure vs Volume ===
subplot(2,2,2)
plot(vol, P_bar, 'LineWidth',1.2)
title('Pressure vs Volume')
xlabel('Volume (m^3)')
ylabel('Pressure (bar)')

%% === Log–Log Pressure vs Volume ===
LogP = log(P_bar);
LogV = repmat(log(vol), 1, nCols);

subplot(2,2,3)
plot(LogV, LogP, 'LineWidth',1.2)
title('Log(Pressure) vs Log(Volume)')
xlabel('log(Volume)')
ylabel('log(Pressure [bar])')

%% === Cylinder Temperature Estimate (uses Pa) ===
R = 287.052;                 % J/(kg·K)
density = 1.167;             % kg/m^3
mass = ((max(vol)-min(vol))/1e6) * density * 0.9;

T = (P_Pa .* vol) / (mass * R) * 0.1;

subplot(2,2,4)
plot(crank, T, 'LineWidth',1.2)
axis([-180 180 0 3000])
title('Estimated Cylinder Temperature')
xlabel('Crank Angle (°)')
ylabel('Temperature (K)')
legend(string(rpmVals)+" rpm",'Location','best')

%% === Rate of Pressure Change dP/dθ (bar/deg) ===
dP = diff(P_bar);
dP(end+1,:) = NaN;

figure(2)
plot(crank, dP, 'LineWidth',1.2)
title('Rate of Pressure Change (dP/dθ)')
xlabel('Crank Angle (°)')
ylabel('dP/dθ (bar/deg)')
legend(string(rpmVals)+" rpm",'Location','best')

%% === Polytropic Index (Fixed Crank Windows) ===
idxComp = find(crank >= -140 & crank <= -30);
idxExp  = find(crank >=   30 & crank <= 140);

nC = zeros(1,nCols);
nE = zeros(1,nCols);

for i = 1:nCols
    % Compression
    pC = polyfit(LogV(idxComp,i), LogP(idxComp,i), 1);
    nC(i) = -pC(1);

    % Expansion
    pE = polyfit(LogV(idxExp,i), LogP(idxExp,i), 1);
    nE(i) = -pE(1);
end

%% === Polytropic Index vs RPM ===
figure(3)
plot(rpmVals, nC, '-o', rpmVals, nE, '-s', 'LineWidth',1.2)
grid on
xlabel('Engine Speed (rpm)')
ylabel('Polytropic Index (n)')
title('Polytropic Index vs RPM')
legend('Compression','Expansion','Location','best')

%% === Print Results ===
fprintf('\nPolytropic Index Results:\n');
for i = 1:nCols
    fprintf('%4d rpm  |  n_c = %.3f   n_e = %.3f\n', ...
            rpmVals(i), nC(i), nE(i));
end