%Luca Gray
%Thermodynamics and Combustion
%ACHA 5011
%University of Wales Trinity Saint David
%March 2026

%load data sheets (rpms.mat) before running

%TASK2
%Base Engine Values
Bore=0.079; %m
Stroke=0.0814; %m
Rod=0.137; %Conrod Length m

%Values from data sheets
%Crank Angles (deg)
Ca1000=rpm1000{10:end,1};
Ca1500=rpm1500{10:end,1};
Ca2000=rpm2000{10:end,1};
Ca2500=rpm2500{10:end,1};
Ca3000=rpm3000{10:end,1};
Ca3500=rpm3500{10:end,1};
Ca4000=rpm4000{10:end,1};
Ca4500=rpm4500{10:end,1};
Ca5000=rpm5000{10:end,1};
Ca5500=rpm5500{10:end,1};
Ca6000=rpm6000{10:end,1};
Ca6500=rpm6000{10:end,1};

%Volumes (cc)
V1000=rpm1000{10:end,2};
V1500=rpm1500{10:end,2};
V2000=rpm2000{10:end,2};
V2500=rpm2500{10:end,2};
V3000=rpm3000{10:end,2};
V3500=rpm3500{10:end,2};
V4000=rpm4000{10:end,2};
V4500=rpm4500{10:end,2};
V5000=rpm5000{10:end,2};
V5500=rpm5500{10:end,2};
V6000=rpm6000{10:end,2};
V6500=rpm6500{10:end,2};

%logV
logV1000 = log(V1000);
logV1500 = log(V1500);
logV2000 = log(V2000);
logV2500 = log(V2500);
logV3000 = log(V3000);
logV3500 = log(V3500);
logV4000 = log(V4000);
logV4500 = log(V4500);
logV5000 = log(V5000);
logV5500 = log(V5500);
logV6000 = log(V6000);
logV6500 = log(V6500);

%Mean Pressures (bar)
P1000=rpm1000{10:end,3};
P1500=rpm1500{10:end,3};
P2000=rpm2000{10:end,3};
P2500=rpm2500{10:end,3};
P3000=rpm3000{10:end,3};
P3500=rpm3500{10:end,3};
P4000=rpm4000{10:end,3};
P4500=rpm4500{10:end,3};
P5000=rpm5000{10:end,3};
P5500=rpm5500{10:end,3};
P6000=rpm6000{10:end,3};
P6500=rpm6500{10:end,3};

%logP
logP1000 = log(P1000);
logP1500 = log(P1500);
logP2000 = log(P2000);
logP2500 = log(P2500);
logP3000 = log(P3000);
logP3500 = log(P3500);
logP4000 = log(P4000);
logP4500 = log(P4500);
logP5000 = log(P5000);
logP5500 = log(P5500);
logP6000 = log(P6000);
logP6500 = log(P6500);

%Polytropic Index
n1000=log(P1000(1:end-1)./P1000(2:end))./log(V1000(2:end)./V1000(1:end-1));
n1500=log(P1500(1:end-1)./P1500(2:end))./log(V1500(2:end)./V1500(1:end-1));
n2000=log(P2000(1:end-1)./P2000(2:end))./log(V2000(2:end)./V2000(1:end-1));
n2500=log(P2500(1:end-1)./P2500(2:end))./log(V2500(2:end)./V2500(1:end-1));
n3000=log(P3000(1:end-1)./P3000(2:end))./log(V3000(2:end)./V3000(1:end-1));
n3500=log(P3500(1:end-1)./P3500(2:end))./log(V3500(2:end)./V3500(1:end-1));
n4000=log(P4000(1:end-1)./P4000(2:end))./log(V4000(2:end)./V4000(1:end-1));
n4500=log(P4500(1:end-1)./P4500(2:end))./log(V4500(2:end)./V4500(1:end-1));
n5000=log(P5000(1:end-1)./P5000(2:end))./log(V5000(2:end)./V5000(1:end-1));
n5500=log(P5500(1:end-1)./P5500(2:end))./log(V5500(2:end)./V5500(1:end-1));
n6000=log(P6000(1:end-1)./P6000(2:end))./log(V6000(2:end)./V6000(1:end-1));
n6500=log(P6500(1:end-1)./P6500(2:end))./log(V6500(2:end)./V6500(1:end-1));

%Rate of change of pressure
dP1000=gradient(P1000)./gradient(Ca1000);
dP1500=gradient(P1500)./gradient(Ca1500);
dP2000=gradient(P2000)./gradient(Ca2000);
dP2500=gradient(P2500)./gradient(Ca2500);
dP3000=gradient(P3000)./gradient(Ca3000);
dP3500=gradient(P3500)./gradient(Ca3500);
dP4000=gradient(P4000)./gradient(Ca4000);
dP4500=gradient(P4500)./gradient(Ca4500);
dP5000=gradient(P5000)./gradient(Ca5000);
dP5500=gradient(P5500)./gradient(Ca5500);
dP6000=gradient(P6000)./gradient(Ca6000);
dP6500=gradient(P6500)./gradient(Ca6500);

%Rate of change of volume
dV1000=gradient(V1000)./gradient(Ca1000);
dV1500=gradient(V1500)./gradient(Ca1500);
dV2000=gradient(V2000)./gradient(Ca2000);
dV2500=gradient(V2500)./gradient(Ca2500);
dV3000=gradient(V3000)./gradient(Ca3000);
dV3500=gradient(V3500)./gradient(Ca3500);
dV4000=gradient(V4000)./gradient(Ca4000);
dV4500=gradient(V4500)./gradient(Ca4500);
dV5000=gradient(V5000)./gradient(Ca5000);
dV5500=gradient(V5500)./gradient(Ca5500);
dV6000=gradient(V6000)./gradient(Ca6000);
dV6500=gradient(V6500)./gradient(Ca6500);

%Torque
T1000=P1000.*dV1000;
T1500=P1500.*dV1500;
T2000=P2000.*dV2000;
T2500=P2500.*dV2500;
T3000=P3000.*dV3000;
T3500=P3500.*dV3500;
T4000=P4000.*dV4000;
T4500=P4500.*dV4500;
T5000=P5000.*dV5000;
T5500=P5500.*dV5500;
T6000=P6000.*dV6000;
T6500=P6500.*dV6500;

%Crank Angle vs Pressure 
%At low engine speed
%Note:due to the Mean Pressure values being incomplete at 1000rpm between
%crank angles -29 and 71, I will be taking low engine speed at 1500rpm
figure
plot(Ca1500,P1500)
xlabel('Crank Angle (deg)')
ylabel('Mean Cylinder Pressure (bar)')
title('Cylinder Pressure vs Crank Angle at 1500RPM')
xlim([-360,360])
grid on

%At medium engine speed
figure
plot(Ca3500,P3500)
xlabel('Crank Angle (deg)')
ylabel('Mean Cylinder Pressure (bar)')
title('Cylinder Pressure vs Crank Angle at 3500RPM')
xlim([-360,360])
grid on

%At high engine speed
figure
plot(Ca6500,P6500)
xlabel('Crank Angle (deg)')
ylabel('Mean Cylinder Pressure (bar)')
title('Cylinder Pressure vs Crank Angle at 6500RPM')
xlim([-360,360])
grid on

%At all engine speeds 
figure
rpm=[1000 1500 2000 2500 3000 3500 4000 4500 5000 5500 6000 6500];
P=[P1000 P1500 P2000 P2500 P3000 P3500 P4000 P4500 P5000 P5500 P6000 P6500];
[RPM, Cagrid]=meshgrid(rpm,Ca1000);
surf(Cagrid,RPM,P,'EdgeColor',[0.7 0.7 0.7])
shading interp
xlabel('Crank Angle (deg)')
ylabel('Engine Speed (RPM)')
zlabel('Mean Cylinder Pressure (bar)')
title('Cylinder Pressure vs Crank Angle')
colormap turbo
colorbar
view(45,30)
xlim([-360 360])

%Pressure vs Volume
%At low engine speed
figure
plot(V1500,P1500)
xlabel('Volume (cc)')
ylabel('Mean Cylinder Pressure (bar)')
title('Pressure vs Volume at 1500RPM')
grid on

%At medium engine speed
figure
plot(V3500,P3500)
xlabel('Volume (cc)')
ylabel('Mean Cylinder Pressure (bar)')
title('Pressure vs Volume at 3500RPM')
grid on

%At high engine speed
figure
plot(V6500,P6500)
xlabel('Volume (cc)')
ylabel('Mean Cylinder Pressure (bar)')
title('Pressure vs Volume at 6500RPM')
grid on

%At all engine speeds
figure
hold on
c = turbo(length(rpm));
plot(V1000,P1000,'Color',c(1,:))
plot(V1500,P1500,'Color',c(2,:))
plot(V2000,P2000,'Color',c(3,:))
plot(V2500,P2500,'Color',c(4,:))
plot(V3000,P3000,'Color',c(5,:))
plot(V3500,P3500,'Color',c(6,:))
plot(V4000,P4000,'Color',c(7,:))
plot(V4500,P4500,'Color',c(8,:))
plot(V5000,P5000,'Color',c(9,:))
plot(V5500,P5500,'Color',c(10,:))
plot(V6000,P6000,'Color',c(11,:))
plot(V6500,P6500,'Color',c(12,:))
xlabel('Volume (cc)')
ylabel('Mean Cylinder Pressure (bar)')
title('Pressure vs Volume')
legend('1000 RPM','1500 RPM','2000 RPM','2500 RPM','3000 RPM',...
    '3500 RPM','4000 RPM','4500 RPM','5000 RPM','5500 RPM',...
    '6000 RPM','6500 RPM')
grid on

%Log Pressure vs Log Volume
%At low engine speed
figure
hold on
plot(logV1500,logP1500)
xlabel('Log Volume')
ylabel('Log Pressure')
title('Log of Pressure vs Log of Volume at 1500RPM')
grid on

%At medium engine speed
figure
hold on
plot(logV3500,logP3500)
xlabel('Log Volume')
ylabel('Log Pressure')
title('Log of Pressure vs Log of Volume at 3500RPM')
grid on
%to find n (polytropic index)
V_points = [3.8 5.8];
P_points = [3.2 0.6];
coeff = polyfit(V_points,P_points,1);
xfit = linspace(3.8,5.8,100);
yfit = polyval(coeff,xfit);
plot(xfit,yfit,'k')
legend('Data','Expansion Fit')
n = -coeff(1);
m = coeff(1);
c = coeff(2);
eqn = sprintf('y = %.2fx + %.2f\nn = %.2f', m, c, n);
text(4.4,2.6, eqn, 'FontSize',12, 'Color','k')
%n=1.3000

%At all engine speeds
figure
rpm=[1000 1500 2000 2500 3000 3500 4000 4500 5000 5500 6000 6500];
logP=[logP1000 logP1500 logP2000 logP2500 logP3000 logP3500 ...
      logP4000 logP4500 logP5000 logP5500 logP6000 logP6500];
[RPM, Vgrid]=meshgrid(rpm,logV1000);
surf(Vgrid,RPM,logP,'EdgeColor',[0.7 0.7 0.7])
shading interp
xlabel('Log Volume')
ylabel('Engine Speed (RPM)')
zlabel('Log Pressure')
title('Log Pressure vs Log Volume')
colormap turbo
colorbar
view(45,30)

%Polytropic Index vs Crank Angle
%At low engine speed
theta1500_n = Ca1500(1:end-1);
figure
plot(theta1500_n,n1500)
xlabel('Crank Angle (deg)')
ylabel('Polytropic Index (n)')
title('Polytropic Index vs Crank Angle at 1500RPM')
grid on
xlim([-360 360])

%At medium engine speed
theta3500_n = Ca3500(1:end-1);
figure
plot(theta3500_n,n3500)
xlabel('Crank Angle (deg)')
ylabel('Polytropic Index (n)')
title('Polytropic Index vs Crank Angle at 3500RPM')
grid on
xlim([-360 360])

%At high engine speed
theta6500_n = Ca6500(1:end-1);
figure
plot(theta6500_n,n6500)
xlabel('Crank Angle (deg)')
ylabel('Polytropic Index (n)')
title('Polytropic Index vs Crank Angle at 6500RPM')
grid on
xlim([-360 360])

%At all engine speeds
figure
rpm=[1000 1500 2000 2500 3000 3500 4000 4500 5000 5500 6000 6500];
n=[n1000 n1500 n2000 n2500 n3000 n3500 n4000 n4500 n5000 n5500 n6000 n6500];
theta=Ca1000(1:end-1);
[RPM, THETA]=meshgrid(rpm,theta);
surf(THETA,RPM,n,'EdgeColor',[0.7 0.7 0.7])
shading interp
xlabel('Crank Angle (deg)')
ylabel('Engine Speed (RPM)')
zlabel('Polytropic Index (n)')
title('Polytropic Index vs Crank Angle')
colormap turbo
colorbar
view(45,30)
xlim([-360 360])

%Rate of Change of Pressure against Crank Angle 
%At low engine speed
figure
plot(Ca1500,dP1500)
xlabel('Crank Angle (deg)')
ylabel('dP/d\theta')
title('Rate of Change of Pressure at 1500RPM')
grid on
xlim([-360 360])

%At medium engine speed 
figure
plot(Ca3500,dP3500)
xlabel('Crank Angle (deg)')
ylabel('dP/d\theta')
title('Rate of Change of Pressure at 3500RPM')
grid on
xlim([-360 360])

%At high engine speed
figure
plot(Ca6500,dP6500)
xlabel('Crank Angle (deg)')
ylabel('dP/d\theta')
title('Rate of Change of Pressure at 6500RPM')
grid on
xlim([-360 360])

% for all RPMs
figure
hold on
surface([Ca1000 Ca1000],[dP1000 dP1000],[zeros(size(dP1000)) zeros(size(dP1000))],[dP1000 dP1000],'EdgeColor','interp','FaceColor','none','LineWidth',2)
surface([Ca1500 Ca1500],[dP1500 dP1500],[zeros(size(dP1500)) zeros(size(dP1500))],[dP1500 dP1500],'EdgeColor','interp','FaceColor','none','LineWidth',2)
surface([Ca2000 Ca2000],[dP2000 dP2000],[zeros(size(dP2000)) zeros(size(dP2000))],[dP2000 dP2000],'EdgeColor','interp','FaceColor','none','LineWidth',2)
surface([Ca2500 Ca2500],[dP2500 dP2500],[zeros(size(dP2500)) zeros(size(dP2500))],[dP2500 dP2500],'EdgeColor','interp','FaceColor','none','LineWidth',2)
surface([Ca3000 Ca3000],[dP3000 dP3000],[zeros(size(dP3000)) zeros(size(dP3000))],[dP3000 dP3000],'EdgeColor','interp','FaceColor','none','LineWidth',2)
surface([Ca3500 Ca3500],[dP3500 dP3500],[zeros(size(dP3500)) zeros(size(dP3500))],[dP3500 dP3500],'EdgeColor','interp','FaceColor','none','LineWidth',2)
surface([Ca4000 Ca4000],[dP4000 dP4000],[zeros(size(dP4000)) zeros(size(dP4000))],[dP4000 dP4000],'EdgeColor','interp','FaceColor','none','LineWidth',2)
surface([Ca4500 Ca4500],[dP4500 dP4500],[zeros(size(dP4500)) zeros(size(dP4500))],[dP4500 dP4500],'EdgeColor','interp','FaceColor','none','LineWidth',2)
surface([Ca5000 Ca5000],[dP5000 dP5000],[zeros(size(dP5000)) zeros(size(dP5000))],[dP5000 dP5000],'EdgeColor','interp','FaceColor','none','LineWidth',2)
surface([Ca5500 Ca5500],[dP5500 dP5500],[zeros(size(dP5500)) zeros(size(dP5500))],[dP5500 dP5500],'EdgeColor','interp','FaceColor','none','LineWidth',2)
surface([Ca6000 Ca6000],[dP6000 dP6000],[zeros(size(dP6000)) zeros(size(dP6000))],[dP6000 dP6000],'EdgeColor','interp','FaceColor','none','LineWidth',2)
surface([Ca6500 Ca6500],[dP6500 dP6500],[zeros(size(dP6500)) zeros(size(dP6500))],[dP6500 dP6500],'EdgeColor','interp','FaceColor','none','LineWidth',2)
colormap turbo
colorbar
xlabel('Crank Angle (deg)')
ylabel('dP/d\theta')
title('Rate of Change of Pressure vs Crank Angle')
grid on
xlim([-360 360])

%Torque against crank angle
%At low engine speed
figure
plot(Ca1500,T1500)
xlabel('Crank Angle (deg)')
ylabel('Torque (Nm)')
title('Torque vs Crank Angle at 1500RPM')
grid on
xlim([-360 360])

%At medium engine speed 
figure
plot(Ca3500,T3500)
xlabel('Crank Angle (deg)')
ylabel('Torque (Nm)')
title('Torque vs Crank Angle at 3500RPM')
grid on
xlim([-360 360])

%At high engine speed
figure
plot(Ca6500,T6500)
xlabel('Crank Angle (deg)')
ylabel('Torque (Nm)')
title('Torque vs Crank Angle at 6500RPM')
grid on
xlim([-360 360])

%At all engine speeds
figure
rpm=[1000 1500 2000 2500 3000 3500 4000 4500 5000 5500 6000 6500];
T=[T1000(:) T1500(:) T2000(:) T2500(:) T3000(:) T3500(:) ...
   T4000(:) T4500(:) T5000(:) T5500(:) T6000(:) T6500(:)];
theta=Ca1000;
[RPM, THETA]=meshgrid(rpm,theta);
surf(THETA,RPM,T,'EdgeColor',[0.7 0.7 0.7])
shading interp
xlabel('Crank Angle (deg)')
ylabel('Engine Speed (RPM)')
zlabel('Torque (Nm)')
title('Torque vs Crank Angle')
colormap turbo
colorbar
view(45,30)
xlim([-360 360])