clc; clear; close all;

%% === Engine parameters (reference) ===
Bore     = 0.104;        % m
Stroke   = 0.132;        % m
Rod      = 0.195;        % m
CrankRad = 0.066;        % m
Vs       = 0.06435505;   % m^3 (swept volume)
Vc       = 0.003900306;  % m^3 (clearance volume)

[fn,fp] = uigetfile({'*.xlsx;*.xls'},'Select claasdata Excel file');
if isequal(fn,0)
    error('No file selected.');
end
filename = fullfile(fp,fn);

raw = readmatrix(filename,'NumHeaderLines',1);
%% === Extract columns ===
crank = raw(:,1);          % Crank angle (deg)
vol   = raw(:,2);          % Cylinder volume (m^3)
Pbar  = raw(:,3:end);      % Pressure traces (bar)

% RPM values from headers: 1000 to 2200 in steps of 100
rpmVals = 1000:100:2200;
nCols   = size(Pbar,2);

% Convert pressure to Pa
pressures = Pbar * 1e5;

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

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

%% === Log–Log Plot ===
LogP = log(pressures);
LogV = log(vol);

subplot(2,2,3)
plot(LogV, LogP)
xlabel('log(Volume)')
ylabel('log(Pressure)')
title('Claas log(P) vs log(V)')

%% === Cylinder Temperature ===
R   = 287.05;       % J/(kg·K)
rho = 1.167;        % kg/m^3 (air)
mass = rho * Vs;    % trapped mass (kg)

T = (pressures .* vol) ./ (mass * R);

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

%% === Rate of Pressure Change (dP/dθ) ===
dP = diff(pressures);
dP(end+1,:) = dP(end,:);   % pad length

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

%% === Combustion timing (approximate) ===
idxSOC = find(crank >= -20 & crank <= 0);
[~,iSOC] = min(dP(idxSOC,:));
SOC = crank(idxSOC(iSOC));

[~,iEOC] = min(dP);
EOC = crank(iEOC);

BurnDuration = EOC - SOC;

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

for i = 1:nCols
    pC = polyfit(LogV(idxComp), LogP(idxComp,i), 1);
    pE = polyfit(LogV(idxExp),  LogP(idxExp,i),  1);
    nC(i) = -pC(1);
    nE(i) = -pE(1);
end

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

%% === Results ===
fprintf('\nStart of Combustion (deg):\n'); disp(SOC)
fprintf('End of Combustion (deg):\n'); disp(EOC)
fprintf('Burn Duration (deg):\n'); disp(BurnDuration)