%% ---------------- Part 3 -----------------%%

%% Simulated pressure (with combustion) v crank angle %%
Bore = 0.079;   
Stroke = 0.0814; 
Rod = 0.137; 
Clearance = 39.27;
gamma = 1.25; 
R = 287; 
rho = 1.167; 
LHV = 44e6; 
AF = 14.7;   
Voleff = 0.9;
Combeff = 0.9; 

clearanceV = Clearance / 1e6;

Crank = -180:1:180;
angle = deg2rad(Crank);

Rad = Stroke/2;
s = Rad*(1 - cos(angle)) + Rod - sqrt(Rod^2 - (Rad*sin(angle)).^2);
pistonA = pi/4 * Bore^2;
Vol = clearanceV + pistonA .* s;

displacementV = pi*(Bore^2)/4 * Stroke;
Air = rho * displacementV * Voleff;
Fuel = Air / AF;
Qtot = Fuel * LHV * Combeff;

a = 5;
m = 2;
Ignition = -12;
IgnDel = 5;
bdur = 30;

SoC = Ignition + IgnDel;
burn_start = find(Crank >= SoC, 1);
burn_end = burn_start + bdur;

x = linspace(0,1,bdur);
Wiebe = 1 - exp(-a * x.^(m+1));

dQr = zeros(size(Crank));
dQr(burn_start:burn_end-1) = Qtot * diff([0 Wiebe]);

P_sim = zeros(size(Crank));
P_sim(1) = 1e5;

n_comp = 1.33;
n_exp = 1.27;

for i = 1:length(Crank)-1

    if Crank(i) <= 0
        n = n_comp;
    else
        n = n_exp;
    end

    P_adiab = P_sim(i) * (Vol(i)/Vol(i+1))^n;
    P_sim(i+1) = P_adiab + (gamma-1)*dQr(i)/Vol(i+1);
end

P_sim = P_sim / 1e5;

figure;
plot(Crank, P_sim, 'r', 'LineWidth', 1.5);
xlabel('Crank Angle (°)');
ylabel('Pressure (bar)');
title('Simulated Cylinder Pressure vs Crank Angle');
grid on;

% We can see the effects of combustion on the cylinder pressure. That is, 
% the pressure increases during compression as the piston approaches TDC, 
% then rises sharply when combustion begins, reaching a peak close to TDC 
% where the burn rate is highest. After the peak, pressure decreases through 
% the expansion stroke as the gases do work on the piston. The smooth shape
%  reflects the idealised nature of the Wiebe-based heat release model. %


%% Real v Simulated pressure %%
data = readmatrix('Book5.xlsx', 'Range', 'B5:P724');
angle = data(:,1);  
pressure = data(:,2:end);
P_real = pressure(:,1); 
P_sim;  
Crank; 
angle; 
P_sim_interp = interp1(Crank, P_sim, angle, 'linear', 'extrap');
figure;
hold on;

plot(angle, P_real, 'k', 'LineWidth', 1.4); 
plot(angle, P_sim_interp, 'r', 'LineWidth', 1.6);  

xlabel('Crank Angle (°)');
ylabel('Pressure (bar)');
title('Real vs Simulated Cylinder Pressure');
legend('Real Pressure', 'Simulated Pressure (Wiebe Model)', ...
       'Location', 'best');
grid on;
hold off;

% We can see the comparison of the real engine pressure against the 
% simulated pressure. Both curves show a rapid pressure rise near TDC with 
% a peak, indicating that the model captures the main features of the 
% combustion process. The real pressure line has a less symmetrical peak 
% due to the flame development, heat transfer, and gas dynamics that are 
% not included in the simplified model. The simulated curve is smoother and
%  produces a sharper peak because it assums ideal polytropic behaviour and
%  a prescribed heat release profile. This comparison highlights the 
% strengths and limitations of the model. % 

%% Simulated pressure (with and no combustion) v crank angle %%
Bore = 0.079;            
Stroke = 0.0814;         
Rod = 0.137;             
Clearance = 39.27;       
gamma = 1.25;            
R = 287;                 
rho = 1.167;             
LHV = 44e6;              
AF = 14.7;               
Voleff = 0.9;            
Combeff = 0.9;           

clearanceV = Clearance / 1e6;

Crank = -180:1:180;
angle = deg2rad(Crank);

Rad = Stroke/2;
s = Rad*(1 - cos(angle)) + Rod - sqrt(Rod^2 - (Rad*sin(angle)).^2);
pistonA = pi/4 * Bore^2;
Vol = clearanceV + pistonA .* s;

displacementV = pi*(Bore^2)/4 * Stroke;
Air = rho * displacementV * Voleff;
Fuel = Air / AF;
Qtot = Fuel * LHV * Combeff;

nComp = 1.33;
nExp  = 1.27;

P_nc = zeros(size(Crank));
P_nc(1) = 1e5;   % Pa

for i = 1:length(Crank)-1
    if Crank(i) <= 0
        n = nComp;
    else
        n = nExp;
    end
    P_nc(i+1) = P_nc(i) * (Vol(i)/Vol(i+1))^n;
end

P_nc_bar = P_nc / 1e5; 

a = 5;
m = 2;
Ignition = -12;
IgnDel = 5;
bdur = 30;

SoC = Ignition + IgnDel;
burn_start = find(Crank >= SoC, 1);
burn_end = burn_start + bdur;

x = linspace(0,1,bdur);
Wiebe = 1 - exp(-a * x.^(m+1));

dQr = zeros(size(Crank));
dQr(burn_start:burn_end-1) = Qtot * diff([0 Wiebe]);

P_sim = zeros(size(Crank));
P_sim(1) = 1e5;   % Pa

for i = 1:length(Crank)-1
    if Crank(i) <= 0
        n = nComp;
    else
        n = nExp;
    end

    P_adiab = P_sim(i) * (Vol(i)/Vol(i+1))^n;
    P_sim(i+1) = P_adiab + (gamma-1)*dQr(i)/Vol(i+1);
end

P_sim_bar = P_sim / 1e5;

figure;
plot(Crank, P_nc_bar, 'b', 'LineWidth', 1.5); hold on;
plot(Crank, P_sim_bar, 'r', 'LineWidth', 1.5);

xlabel('Crank Angle (°)');
ylabel('Pressure (bar)');
title('Simulated Cylinder Pressure with and no Combustion vs Crank Angle');
legend('No Combustion','With Combustion');
grid on;

% We can see the ideal polytropic compression and expansion process without
%  any combustion. The pressure rises smoothly as the piston compresses the 
% trapped air, and falls during expansion as the volume increases. Because 
% no heat is added, the peak pressure is much lower than in the with 
% combustion, and we see a much smoother and symmertrical line. This graph 
% provides a baseline for understanding how combustion alters the pressure 
% trace. % 

%% Simulated pressure (no combustion) v crank angle %%
Bore = 0.079;
Stroke = 0.0814;
Rod = 0.137; 
Clearance = 39.27;
nComp = 1.33;
nExp  = 1.27; 

clearanceV = Clearance / 1e6;
Crank = -180:1:180;
angle = deg2rad(Crank); 
Rad = Stroke/2;
s = Rad*(1 - cos(angle)) + Rod - sqrt(Rod^2 - (Rad*sin(angle)).^2);
pistonA = pi/4 * Bore^2;
Vol = clearanceV + pistonA .* s; 

P_nc = zeros(size(Crank));
P_nc(1) = 1e5;

for i = 1:length(Crank)-1
    
    if Crank(i) <= 0
        n = nComp;  
    else
        n = nExp; 
    end

    P_nc(i+1) = P_nc(i) * (Vol(i)/Vol(i+1))^n;
end

P_nc_bar = P_nc / 1e5;

figure;
plot(Crank, P_nc_bar, 'b', 'LineWidth', 1.5);
xlabel('Crank Angle (°)');
ylabel('Pressure (bar)');
title('Simulated Cylinder Pressure (No Combustion) vs Crank Angle');
grid on;

% Plotting both the with and without combustion together highlights the 
% effects of combustion on the pressure cycle. We can see that they follow 
% the same path during early compression, but once combustion begins, the 
% simulated combustion curve rises sharply above the non-combustion curve. 
% This is due to the combustion assing a large amount of heat to the 
% cylinder gases over a short crank angle interval. This rapid heat release
%  increases the gas temperature, and the pressure, near the TDC where the 
% cylinder volume is at its minimum. %


%% Burn rate from Wiebe function %%
Crank = -180:1:180;      
a = 5;         
m = 2;      
SoC = -15;  
bdur = 40;  

Wiebe = 1 - exp(-a*((Crank - SoC)./bdur).^(m+1)); 

for i = 1:1:360
    if (Crank(i) < SoC)
        stop = 0;
    else
        stop = 1;
    end

    MFB(i) = (a*(m+1)/bdur)*((Crank(i)-SoC)/bdur)^m * ...
             exp(-a*((Crank(i)-SoC)/bdur)^(m+1)) * stop;
end

MFB(361) = 0;
burn_rate = MFB;

figure;
plot(Crank, burn_rate, 'r', 'LineWidth', 1.8);
xlabel('Crank Angle (°)');
ylabel('Burn Rate (per degree)');
title('Burn Rate from Wiebe Function');
grid on;

% We can see how quickly the fuel-air mixture burns during the combustion 
% phase. The burn rate rises sharply once the crank angle reaches the start
% of combustion point. It then peaks near TDC, where heat release is most 
% intense, and then falls as combustion completes. The single peak shape 
% reflects the Wiebe model's assumption of smooth, continuous heat release 
% over the crank angle duration. %

%% Limitations and Improvements %%
% The simulation coded works well for showing the general shape of the 
% pressure cycle, but it still has a few clear limitations. The model 
% assumes ideal polytropic compression and expansion, and the Wiebe 
% function gives a smooth, simplified burn that doesn't include real‑engine
% effects like heat losses, turbulence, flame speed variation, or residual 
% gases. These missing factors explain why the simulated pressure peak is 
% sharper and more symmetrical than the real trace. The model also uses 
% fixed values for gamma and the polytropic indices, even though these 
% actually change throughout the cycle. To improve the accuracy, the 
% simulation could include temperature‑dependent gas properties, heat 
% transfer models, variable polytropic exponents, and a more detailed 
% combustion model that accounts for flame development and burn rate 
% changes. Adding these would make the simulated pressure curve match the 
% real data more closely. %