function [P,S,Vol,V,CV] = Engine_Dynamics(RPM1, RPM2, RPM3, RPM4, RPM5, RPM6, Crank, Rad, Rod, CR, Bore, RPM, CAmin, CAmax)




P = [RPM1, RPM2, RPM3, RPM4, RPM5, RPM6];

%Pressure graph 
figure
plot(Crank, P)
% plot(P)
title('Pressure vs Crank Angle')
xlabel('Crank (Degrees)')
ylabel('Pressure (Bar)')
legend(RPM)

S= Rad + Rod - (Rad .* cosd(Crank)) - (Rod .* cosd(asind((Rad/Rod) .* sind(Crank))));  %s 
V = (pi() * ((Bore^2)/4) * (Rad*2));
CV = V/(CR - 1);
Vol= (pi() * ((Bore^2)/4) * S) + CV;      

figure
plot(Crank, Vol)
title('Volume vs Crank Angle')
xlabel('Crank (Degrees)')
ylabel('Volume (m^3)')

figure
plot(Vol, P)
title('Volume vs Pressure')
xlabel('Volume (m^3)')
ylabel('Pressure (Bar)')
legend(RPM)

figure
plot(log(Vol), log(P))
title('LogVolume vs LogPressure')
xlabel('Volume (m^3)')
ylabel('Pressure (Bar)')
legend(RPM)

windowN = 11;
P_smooth = movmean(P, windowN, 1);

comb = (Crank >= CAmin) & (Crank <= CAmax);

Cr = Crank(comb);
P2 = P_smooth(comb,:);




dP_dCA = zeros(size(P2));
for k = 1:size(P2,2)
    dP_dCA(:,k) = gradient(P2(:,k), Cr);
end



[maxPPRR, idxLocal] = max(dP_dCA, [], 1);
CA_maxPPRR = Cr(idxLocal);

yMax = max(dP_dCA(:));
yMin = min(dP_dCA(:));
pad = 0.10*(yMax - yMin + eps);

figure
plot(Cr, dP_dCA, 'LineWidth', 1.5)
grid on
xline(0,'-','TDC','LabelVerticalAlignment','bottom');
xlim([0 40])
ylim([yMin-pad, yMax+pad])
title('Pressure Rise Rate (dP/d\theta) vs Crank Angle')
xlabel('Crank Angle (deg)')
ylabel('dP/d\theta (bar/deg)')
legend(string(RPM), 'Location','best')


hold on
for k = 1:numel(CA_maxPPRR)
    h = plot(CA_maxPPRR(k), maxPPRR(k), 'o', 'MarkerSize', 6, 'LineWidth', 1.5);
    h.HandleVisibility = 'off';
end
hold off

T = table(string(RPM(:)), maxPPRR(:), CA_maxPPRR(:),'VariableNames', {'RPM','Max_dP_dCA_bar_per_deg','CA_at_Max_deg'});
disp("Max Pressure Rise Rate Results (Smoothed)")
disp(T)


P_Pa = P * 1e5;                    % convert bar -> Pa


dV_dCA = gradient(Vol, Crank);     % m^3/deg

Wi = zeros(1, size(P_Pa,2));      
for k = 1:size(P_Pa,2)
    Wi(k) = trapz(Crank, P_Pa(:,k) .* dV_dCA);   
end

Vd = max(Vol) - min(Vol);          % swept volume per cylinder (m^3)

IMEP_Pa  = Wi ./ Vd;
IMEP_bar = IMEP_Pa / 1e5;

% Display table 
disp("IMEP Results")
T_imep = table(string(RPM(:)), IMEP_bar(:), Wi(:), ...
    'VariableNames', {'RPM','IMEP_bar','IndicatedWork_J'});
disp(T_imep)


figure
plot(str2double(string(RPM)), IMEP_bar, 'o-', 'LineWidth', 1.5)
grid on
title('IMEP vs RPM')
xlabel('RPM')
ylabel('IMEP (bar)')


end