clc; clear; close all;

%% === Load Excel data (robust) ===
[fn,fp] = uigetfile('*.xlsx','Select husqdata Excel file');
if isequal(fn,0)
    error('No file selected.');
end
filename = fullfile(fp,fn);

raw = readmatrix(filename,'NumHeaderLines',2);

%% === Extract data ===
crank = raw(:,1);          % Crank angle (deg)
vol   = raw(:,2);          % Cylinder volume (cc)
Pbar  = raw(:,3:end);      % Pressure traces (bar)

rpmVals = 2000:1000:10000;
nCols = size(Pbar,2);

% Unit conversions
pressures = Pbar * 1e5;    % bar → Pa
vol_m3 = vol * 1e-6;       % cc → m³

%% === Pressure vs Crank Angle ===
figure(1)
subplot(2,2,1)
plot(crank, Pbar, 'LineWidth',1.1)
title('Husqvarna Cylinder 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, Pbar, 'LineWidth',1.1)
title('Husqvarna Pressure vs Volume')
xlabel('Volume (cc)')
ylabel('Pressure (bar)')

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

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

%% === Cylinder Temperature ===
R = 287.05;              % J/(kg·K)
rho = 1.167;             % kg/m³
Vs  = 248.870438;        % swept volume (cc)
mass = Vs*1e-6 * rho;    % kg

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

subplot(2,2,4)
plot(crank, T)
axis([-180 180 0 3000])
title(' Husqvarna Cylinder Temperature')
xlabel('Crank Angle (°)')
ylabel('Temperature (K)')
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)
title('Husqvarna Rate of Pressure Change')
xlabel('Crank Angle (°)')
ylabel('dP/dθ (Pa/deg)')
legend(string(rpmVals)+" rpm",'Location','best')

%% === Combustion Timing (approximate) ===
thetaSOC = find(crank >= -30 & crank <= 0);
[~,iSOC] = min(dP(thetaSOC,:));
SOC = crank(thetaSOC(iSOC));

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

BurnDuration = EOC - SOC;

%% === Polytropic Index ===
idxComp = find(crank >= -140 & crank <= -20);
idxExp  = find(crank >=  20 & 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
legend('Compression','Expansion','Location','best')
xlabel('Engine Speed (rpm)')
ylabel('Polytropic Index (n)')
title('Husqvarna Polytropic Index vs RPM')

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