clc,clear,close all

Data = readmatrix('1600 sigma Pressure Data.xlsx');

RPM_Range= Data (1,2:7);
Pres_Data=Data(4:end,2:7);  %Bar
Crank=Data(4:end,1);    %Crank angle (degrees)
RPM_Range2=(1000:1000:9000);
Pres_Data=Data(4:end,2:7);  %Bar
for i =1:length(Crank)
Pres_Data2(i,:)=interp1(RPM_Range,Pres_Data(i,:),RPM_Range2,'linear','extrap');
end
Pres_Data=Pres_Data2;
RPM_Range=RPM_Range2;

Bore=79/1000;   %Bore diameter m
Crank_Radius =40.7/1000;     %Crank radius m
Rod=134/1000;   % Con rod length m
Piston_mass=258/1000;   % pisotn mass kg
Rod_mass = 372/1000;    %con rod mass with bolts kg
Ring_mass = 9/1000;     %Ring masses kg
Clips_mass = 5/1000;    %Clips mass kg
Pin_mass = 80/1000;      %Pin mass kg
Exposed_Pin_mass = 23.28/1000;   %Exposed section of the pins mass
Stroke = Crank_Radius*2;     %Stroke m

%%%% No Clearance volume %%%%%

S = (Crank_Radius+Rod) - ((Crank_Radius*cosd(Crank))+(Rod*(cosd(asin((Crank_Radius*sind(Crank))/Rod)))));       %Calculate the value of 's' - Distance the Piston has Moved
Vol= (pi*(Bore^2)/4)*S;     %Cylinder Volume (m^3)
Vol = Vol/1e-6;     %Converting m^3 to cc

Vol3d = repmat(Vol,1,9);
Crank3d = repmat(Crank,1,9);
RPM_Range3d = repmat(RPM_Range,720,1);

figure;
plot(Crank,Vol)
title('Crank angle against volume')
xlabel('Crank angle (Degrees)')
ylabel('Volume (cc)')

figure;
surf(Crank3d,RPM_Range3d,Pres_Data)
shading interp
title('Crank angle against Rpm against Pressure')
xlabel('Crank angle (Degrees)')
ylabel ('RPM')
zlabel ('Pressure (Bar)')

for i = 1:719
    vel(i) = S(i+1)-S(i);
end
vel(end+1)=0;
vel=vel';

figure;
yyaxis left
plot(Crank,S)
ylabel('Piston pisition from TDC (m)')
hold on
yyaxis right
plot(Crank,vel)
ylabel('Piston velocity (m per degree)')
legend('Piston Position','Piston Velocity')
title('Piston position against velocity')

Vel_RPM=vel.*6.*RPM_Range;        %Piston velocity in account for the RPM

figure;
surf(Crank3d,RPM_Range3d,Vel_RPM)
shading interp
title('Crank angle against Rpm against Velocity of the piston')
xlabel('Crank angle (Degrees)')
ylabel ('RPM')
zlabel ('Velocity (m/s)')

% Angular velocity of crankshaft
Ang_vel= (RPM_Range*2*pi)./60;    % Degrees/sec

% Angular velocity of the con rod
Ang_Vel_Con_Rod=(Crank_Radius.*Ang_vel.*sind(Crank))+((Crank_Radius.^2.*Ang_vel)./(2.*Rod).*sind(2*Crank));

figure;
surf(Crank3d,RPM_Range3d,Ang_Vel_Con_Rod)
shading interp
title('Crank angle against Rpm against Angular velcoity of the rod')
xlabel('Crank angle (Degrees)')
ylabel ('RPM')
zlabel ('Angular Velocity (m/s)')

% Angular Acceleration of the con rod

Ang_Accel_Con_Rod = (Crank_Radius.*(Ang_vel.^2).*cosd(Crank))+(((Crank_Radius.^2)./Rod).*(Ang_vel.^2).*cosd(2.*Crank));

figure;
surf(Crank3d,RPM_Range3d,Ang_Accel_Con_Rod)
shading interp
title('Crank angle against Rpm against Angular Accel of the rod')
xlabel('Crank angle (Degrees)')
ylabel ('RPM')
zlabel ('Angular Accel (m/s/s)')

% Calculating inertia

% Reciprocating mass (Kg)
Recip_Mass= Pin_mass+Piston_mass+Ring_mass+((1./3).*Rod_mass);

Inertia_Con_Rod = Recip_Mass.*Ang_Accel_Con_Rod;

figure;
surf(Crank3d,RPM_Range3d,Inertia_Con_Rod)
shading interp
title('Crank angle against Rpm against Inertia of the rod')
xlabel('Crank angle (Degrees)')
ylabel ('RPM')
zlabel ('Inertia force (N)')

%Piston Forces

%Force from the gas

Area = (pi.*Bore.^2)./4;

Force_Gas = Pres_Data.*1e5.*Area;       % Force of the gas

Force_Net_Compressive= Inertia_Con_Rod-Force_Gas;    % Net forces on the assembly

figure;
surf(Crank3d,RPM_Range3d,Force_Net_Compressive)
shading interp
title('Crank angle against Rpm against Net force of the assembly')
xlabel('Crank angle (Degrees)')
ylabel ('RPM')
zlabel ('Net force (N)')

% Calculating phi (Angle from the rod to the vertical position of the
% centre of the piston (OBA))

phi = asind((Crank_Radius.*sind(Crank))./Rod);

Force_Compressive_Con_Rod=Force_Net_Compressive./cosd(phi);

figure;
surf(Crank3d,RPM_Range3d,Force_Compressive_Con_Rod)
shading interp
title('Crank angle against Rpm against Net force of the rod')
xlabel('Crank angle (Degrees)')
ylabel ('RPM')
zlabel ('Net force (N)')

% Clalculating force on the side of the piston

Force_Compressive_Side=Force_Compressive_Con_Rod.*tand(phi);

figure;
surf(Crank3d,RPM_Range3d,Force_Compressive_Side)
shading interp
title('Crank angle against Rpm against the force of the side of the piston')
xlabel('Crank angle (Degrees)')
ylabel ('RPM')
zlabel ('Net force (N)')

% Friction force

Coeff_Friction = 0.1;    %Coefficient of friction for the lubricating oil

Friction_Force_Compressive= Force_Compressive_Side.*Coeff_Friction;

figure;
surf(Crank3d,RPM_Range3d,Friction_Force_Compressive)
shading interp
title('Crank angle against Rpm against the friction force of the piston')
xlabel('Crank angle (Degrees)')
ylabel ('RPM')
zlabel ('Friction force (N)')

%% Stress of the rod

CSA_Con_Rod =139.655.*1e-6;   %m^2
UTS = 1000*1e6;  %MPa
End_Lim=300*1e6;   %Endurance limit
Surf_Finish_Press=-100*1e6;  %%%%%%%%%%%%%%%%%%%%%%% -100 is a place holder %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Stress_Compressive=Force_Compressive_Con_Rod./CSA_Con_Rod;

Stress_Max=max(Stress_Compressive);
Stress_Min=min(Stress_Compressive);

Stress_Range=Stress_Max-Stress_Min;

Stress_Mean=(Stress_Max+Stress_Min)./2;

Stress_Allowable=2.*End_Lim.*(1-(Stress_Mean./UTS));

FSF=Stress_Allowable./Stress_Range;     %Fatigue safety factor

figure;
surf(Crank3d,RPM_Range3d,Stress_Compressive)
shading interp
title('Crank angle against Rpm against the stress (compressive) of the piston')
xlabel('Crank angle (Degrees)')
ylabel ('RPM')
zlabel ('Stress (N/m^2)')

%GOODMANS Diagram

% Diagram for no surface finish

%Upper Endurance limit plot
Y=[End_Lim-(UTS-End_Lim),End_Lim,UTS];
X=[-UTS,0,UTS];

%Lower Endurance limit plot
X1=[-UTS,0,UTS];
Y1=[-End_Lim-(UTS+End_Lim),-End_Lim,UTS];

subplot(1,2,1)
hold on
plot(X,Y)
plot(X1,Y1)

for i=1:length(RPM_Range)
% Actual range plot
X2(1,i)=Stress_Mean(i);
X2(2,i)=Stress_Mean(i);
Y2(1,i)=Stress_Min(i);
Y2(2,i)=Stress_Max(i);

%Allowable range plot
X3(1,i)=Stress_Mean(i);
X3(2,i)=Stress_Mean(i);
Y3(1,i)=Stress_Mean(i)-(Stress_Range(i)./2.*FSF(i));
Y3(2,i)=Stress_Mean(i)+(Stress_Range(i)./2.*FSF(i));
end

plot(X3,Y3, 'y','LineWidth', 5)
plot(X2,Y2, 'LineWidth',2.5)
%Line going through 0,0
Y4=[-UTS,0,UTS];
X4=[-UTS,0,UTS];
plot(Y4,X4,'g--')
%UTS line
Y5=[UTS,-End_Lim-(UTS+End_Lim)];
X5=[UTS,UTS];
Y6=[UTS,UTS];
X6=[UTS,-UTS];
plot(X5,Y5,'k--')
plot(X6,Y6,'k--')

title('Goodmans diagram for no surface finish')
legend('Upper endurance limit','Lower endurance limit', 'Allowable range @1000 rpm','Allowable range @2000 rpm','Allowable range @3500 rpm','Allowable range @4000 rpm','Allowable range @5000 rpm','Allowable range @6000 rpm','Actual range @1000 rpm','Actual range @2000 rpm','Actual range @3500 rpm','Actual range @4000 rpm','Actual range @5000 rpm','Actual range @6000 rpm','0 point','UTS')
xlabel('Mean Stress (N/m)')
ylabel('Stress (N/m)')
a=(min(Stress_Mean).*3);
b=(UTS.*1.1);
xlim([a b])
a=(min(Stress_Min).*3);
b=(UTS.*1.1);
ylim([a b])
hold off

Stress_Min_New=Stress_Min+Surf_Finish_Press;
Stress_Max_New=Stress_Max+Surf_Finish_Press;

Stress_Mean_New=(Stress_Max_New+Stress_Min_New)./2;

Stress_Range_New=Stress_Max_New-Stress_Min_New;

Stress_Allowable_New=2.*End_Lim.*(1-(Stress_Mean_New./UTS));

FSF_New = Stress_Allowable_New./Stress_Range_New;

%GOODMANS Diagram

% Diagram for a surface finish

%Upper Endurance limit plot
Y=[End_Lim-(UTS-End_Lim),End_Lim,UTS];
X=[-UTS,0,UTS];

%Lower Endurance limit plot
X1=[-UTS,0,UTS];
Y1=[-End_Lim-(UTS+End_Lim),-End_Lim,UTS];

subplot(1,2,2)
hold on
plot(X,Y)
plot(X1,Y1)


for i=1:length(RPM_Range)
% Actual range plot
X2(1,i)=Stress_Mean_New(i);
X2(2,i)=Stress_Mean_New(i);
Y2(1,i)=Stress_Min_New(i);
Y2(2,i)=Stress_Max_New(i);

%Allowable range plot
X3(1,i)=Stress_Mean_New(i);
X3(2,i)=Stress_Mean_New(i);
Y3(1,i)=Stress_Mean_New(i)-(Stress_Range_New(i)./2.*FSF_New(i));
Y3(2,i)=Stress_Mean_New(i)+(Stress_Range_New(i)./2.*FSF_New(i));
end

plot(X3,Y3, 'y','LineWidth', 5)
plot(X2,Y2, 'LineWidth',2.5)
%Line going through 0,0
Y4=[-UTS,0,UTS];
X4=[-UTS,0,UTS];
plot(Y4,X4,'g--')
%UTS line
Y5=[UTS,-End_Lim-(UTS+End_Lim)];
X5=[UTS,UTS];
Y6=[UTS,UTS];
X6=[UTS,-UTS];
plot(X5,Y5,'k--')
plot(X6,Y6,'k--')

title('Goodmans diagram with a surface finish')
legend('Upper endurance limit','Lower endurance limit', 'Allowable range @1000 rpm','Allowable range @2000 rpm','Allowable range @3500 rpm','Allowable range @4000 rpm','Allowable range @5000 rpm','Allowable range @6000 rpm','Actual range @1000 rpm','Actual range @2000 rpm','Actual range @3500 rpm','Actual range @4000 rpm','Actual range @5000 rpm','Actual range @6000 rpm','0 point','UTS')
xlabel('Mean Stress (N/m)')
ylabel('Stress (N/m)')
a=(min(Stress_Mean_New).*3);
b=(UTS.*1.1);
xlim([a b])
a=(min(Stress_Min_New).*3);
b=(UTS.*1.1);
ylim([a b])
hold off

%%  PIN and Torque

L1=0.00775;
L2=0.007;
D2=0.021;
D1=0.0123;


Force_Pin_Inertia=Exposed_Pin_mass.*Ang_Accel_Con_Rod;        %N

Bend_Mom_Pin_Inertia=Force_Pin_Inertia.*L1;       %Nm

Piston_Comp_Mass=Clips_mass+Piston_mass+Ring_mass;    %Mass of components of the piston

Force_Pin_Piston=Piston_Comp_Mass.*Ang_Accel_Con_Rod;       %N

Force_Net=Force_Pin_Piston-Force_Gas;

Bend_Mom_Net_Force=(Force_Net./2).*L2;       %Nm

Bend_Mom=(Bend_Mom_Net_Force+Bend_Mom_Pin_Inertia);

Tau=((Force_Net./2)+Force_Pin_Inertia)./((pi./4).*((D2.^2)-(D1.^2)));

I=(pi./64).*(D2.^4-D1.^4);

y=D2/2;      %Check this      %Distnace neutral axis

Bend_Stress=(Bend_Mom.*y)./I;

Stress_Result=sqrt(((Bend_Stress).^2)+((2.*Tau).^2));

figure;
surf(Crank3d,RPM_Range3d,Tau)
shading interp
title('Crank angle against Rpm against Tau of the Piston')
xlabel('Crank angle (Degrees)')
ylabel ('RPM')
zlabel ('Tau (N/m^2)')

figure;
surf(Crank3d,RPM_Range3d,Bend_Mom)
shading interp
title('Crank angle against Rpm against Bending moment of the Piston')
xlabel('Crank angle (Degrees)')
ylabel ('RPM')
zlabel ('Bending moment (Nm)')

figure;
surf(Crank3d,RPM_Range3d,Stress_Result)
shading interp
title('Crank angle against Rpm against Resultant stress of the Piston')
xlabel('Crank angle (Degrees)')
ylabel ('RPM')
zlabel ('Resultant stress (Nm)')

figure;
hold on
plot(Crank,Bend_Stress(:,6))
plot(Crank,Tau(:,6))
plot(Crank,Stress_Result(:,6))
legend('Bending moment','Tau','Resultant stress')

%% Torque

Torque_Inertia_Cyl_1=(Inertia_Con_Rod.*tand(phi).*Crank_Radius.*cosd(Crank))-(Inertia_Con_Rod.*Crank_Radius.*sind(Crank));
Torque_Inertia_Cyl_2=circshift(Torque_Inertia_Cyl_1,-180);
Torque_Inertia_Cyl_3=circshift(Torque_Inertia_Cyl_1,-360);
Torque_Inertia_Cyl_4=circshift(Torque_Inertia_Cyl_1,180);

figure;
plot(Crank,Torque_Inertia_Cyl_1(:,6))
hold on
plot(Crank,Torque_Inertia_Cyl_2(:,6))
plot(Crank,Torque_Inertia_Cyl_3(:,6),'--')
plot(Crank,Torque_Inertia_Cyl_4(:,6),'--')
xlabel('Crank angle (degrees)')
ylabel('Torque (Nm)')
title('Torque inertia for the individual cylinders @6000rpm')
legend('Cyl 1','Cyl 3','Cyl 2','Cyl 4')

Torque_Inertia = Torque_Inertia_Cyl_1+Torque_Inertia_Cyl_2+Torque_Inertia_Cyl_3+Torque_Inertia_Cyl_4;

figure;
plot(Crank,Torque_Inertia_Cyl_1(:,6))
hold on
plot(Crank,Torque_Inertia_Cyl_2(:,6))
plot(Crank,Torque_Inertia_Cyl_3(:,6),'--')
plot(Crank,Torque_Inertia_Cyl_4(:,6),'--')
plot(Crank,Torque_Inertia(:,6))
xlabel('Crank angle (degrees)')
ylabel('Torque (Nm)')
title('Torque inertia for the individual cylinders and the crank @6000rpm')
legend('Cyl 1','Cyl 3','Cyl 2','Cyl 4','Crank')

figure;
surf(Crank3d,RPM_Range3d,Torque_Inertia)
shading interp
title('Crank angle against Rpm against inertia torque of the Piston')
xlabel('Crank angle (Degrees)')
ylabel ('RPM')
zlabel ('Inertia torque (Nm)')

% Cross plane inertia Torque
Cross_Plane_Torque_Inertia_Cyl_2=circshift(Torque_Inertia_Cyl_1,-180);
Cross_Plane_Torque_Inertia_Cyl_3=circshift(Torque_Inertia_Cyl_1,-270);
Cross_Plane_Torque_Inertia_Cyl_4=circshift(Torque_Inertia_Cyl_1,270);
Cross_Plane_Torque_Inertia = Torque_Inertia_Cyl_1+Cross_Plane_Torque_Inertia_Cyl_2+Cross_Plane_Torque_Inertia_Cyl_3+Cross_Plane_Torque_Inertia_Cyl_4;

figure;
plot(Crank,Torque_Inertia_Cyl_1(:,6))
hold on
plot(Crank,Cross_Plane_Torque_Inertia_Cyl_2(:,6))
plot(Crank,Cross_Plane_Torque_Inertia_Cyl_3(:,6))
plot(Crank,Cross_Plane_Torque_Inertia_Cyl_4(:,6))
plot(Crank,Cross_Plane_Torque_Inertia(:,6))
xlabel('Crank angle (degrees)')
ylabel('Torque (Nm)')
title('Torque inertia for the individual cylinders and the crank @6000rpm for a cross plane engine')
legend('Cyl 1','Cyl 3','Cyl 2','Cyl 4','Crank')

%% Inertia Torque for a V8 engine

V8_RPM_Range=1000:500:5000;
V8_Ang_vel=(V8_RPM_Range*2*pi)./60;
V8_Rod=145./1000; % looked up V8 CONROD LENGTH usually between 5.0 to 6.5 inches https://www.lunatipower.com/rod-length/
V8_Crank_Radius=Crank_Radius.*(V8_Rod./V8_Rod);
V8_phi=asind((V8_Crank_Radius.*sind(Crank))./V8_Rod);
V8_Ang_Accel_Con_Rod = ((V8_Crank_Radius).*(V8_Ang_vel.^2).*cosd(Crank))+(((V8_Crank_Radius.^2)./V8_Rod).*(V8_Ang_vel.^2).*cosd(2.*Crank));

V8_Inertia_Con_Rod = Recip_Mass.*V8_Ang_Accel_Con_Rod;
V8_Torque_Inertia_Cyl_1=(V8_Inertia_Con_Rod.*tand(V8_phi).*(V8_Crank_Radius).*cosd(Crank))-(V8_Inertia_Con_Rod.*(V8_Crank_Radius).*sind(Crank));

%Standard V8

V8_Torque_Inertia_Cyl_2=circshift(V8_Torque_Inertia_Cyl_1,-180);
V8_Torque_Inertia_Cyl_3=circshift(V8_Torque_Inertia_Cyl_1,-360);
V8_Torque_Inertia_Cyl_4=circshift(V8_Torque_Inertia_Cyl_1,180);
V8_Torque_Inertia_Cyl_5=circshift(V8_Torque_Inertia_Cyl_1,(0));
V8_Torque_Inertia_Cyl_6=circshift(V8_Torque_Inertia_Cyl_1,(360));
V8_Torque_Inertia_Cyl_7=circshift(V8_Torque_Inertia_Cyl_1,(180));
V8_Torque_Inertia_Cyl_8=circshift(V8_Torque_Inertia_Cyl_1,(-180));
V8_Torque_Inertia = V8_Torque_Inertia_Cyl_1+V8_Torque_Inertia_Cyl_2+V8_Torque_Inertia_Cyl_3+V8_Torque_Inertia_Cyl_4+V8_Torque_Inertia_Cyl_5+V8_Torque_Inertia_Cyl_6+V8_Torque_Inertia_Cyl_7+V8_Torque_Inertia_Cyl_8;
figure;
plot(Crank,V8_Torque_Inertia_Cyl_1(:,6))
hold on
plot(Crank,V8_Torque_Inertia_Cyl_2(:,6))
plot(Crank,V8_Torque_Inertia_Cyl_3(:,6))
plot(Crank,V8_Torque_Inertia_Cyl_4(:,6))
plot(Crank,V8_Torque_Inertia_Cyl_5(:,6),'--')
plot(Crank,V8_Torque_Inertia_Cyl_6(:,6),'--')
plot(Crank,V8_Torque_Inertia_Cyl_7(:,6),'--')
plot(Crank,V8_Torque_Inertia_Cyl_8(:,6),'--')
plot(Crank,V8_Torque_Inertia(:,6),'k--')
xlabel('Crank angle (degrees)')
ylabel('Torque (Nm)')
title('Torque inertia for the individual cylinders and the crank @3500 rpm for a V8 standard engine')
legend('Cyl 1','Cyl 5','Cyl 3','Cyl 7','Cyl 4','Cyl 8','Cyl 2','Cyl 6','Crank')

%Cross Plane V8

V8_Cross_Plane_Torque_Inertia_Cyl_2=circshift(V8_Torque_Inertia_Cyl_1,-180);
V8_Cross_Plane_Torque_Inertia_Cyl_3=circshift(V8_Torque_Inertia_Cyl_1,-270);
V8_Cross_Plane_Torque_Inertia_Cyl_4=circshift(V8_Torque_Inertia_Cyl_1,270);
V8_Cross_Plane_Torque_Inertia_Cyl_5=circshift(V8_Torque_Inertia_Cyl_1,-90);
V8_Cross_Plane_Torque_Inertia_Cyl_6=circshift(V8_Torque_Inertia_Cyl_2,90);
V8_Cross_Plane_Torque_Inertia_Cyl_7=circshift(V8_Torque_Inertia_Cyl_3,360);
V8_Cross_Plane_Torque_Inertia_Cyl_8=circshift(V8_Torque_Inertia_Cyl_4,180);
V8_Cross_Plane_Torque_Inertia = V8_Torque_Inertia_Cyl_1+V8_Cross_Plane_Torque_Inertia_Cyl_2+V8_Cross_Plane_Torque_Inertia_Cyl_3+V8_Cross_Plane_Torque_Inertia_Cyl_4+V8_Cross_Plane_Torque_Inertia_Cyl_5+V8_Cross_Plane_Torque_Inertia_Cyl_6+V8_Cross_Plane_Torque_Inertia_Cyl_7+V8_Cross_Plane_Torque_Inertia_Cyl_8;
figure;
plot(Crank,V8_Torque_Inertia_Cyl_1(:,6))
hold on
plot(Crank,V8_Cross_Plane_Torque_Inertia_Cyl_2(:,6))
plot(Crank,V8_Cross_Plane_Torque_Inertia_Cyl_3(:,6))
plot(Crank,V8_Cross_Plane_Torque_Inertia_Cyl_4(:,6))
plot(Crank,V8_Cross_Plane_Torque_Inertia_Cyl_5(:,6))
plot(Crank,V8_Cross_Plane_Torque_Inertia_Cyl_6(:,6))
plot(Crank,V8_Cross_Plane_Torque_Inertia_Cyl_7(:,6))
plot(Crank,V8_Cross_Plane_Torque_Inertia_Cyl_8(:,6))
plot(Crank,V8_Cross_Plane_Torque_Inertia(:,6),'k--')
xlabel('Crank angle (degrees)')
ylabel('Torque (Nm)')
title('Torque inertia for the individual cylinders and the crank @3500 rpm for a V8 cross plane engine')
legend('Cyl 1','Cyl 5','Cyl 3','Cyl 7','Cyl 4','Cyl 8','Cyl 2','Cyl 6','Crank')

%% Cross plane against Standard for inline 4 and V8

figure;
plot(Crank,Torque_Inertia(:,6))
hold on
plot(Crank,Cross_Plane_Torque_Inertia(:,6))
plot(Crank,V8_Torque_Inertia(:,6))
plot(Crank,V8_Cross_Plane_Torque_Inertia(:,6))
legend('Standard engine','Cross plane engine','Standard V8 engine','Cross plane V8 engine')
title('Standard vs Cross plane engine setup Torque inertia @6000rpm of the crank for the inline 4 engines and 3000 rpm for the V8 engines')
xlabel('Crank angle (degrees)')
ylabel('Torque Inertia (Nm)')
hold off
%% Torque of the gas

Torque_Gas_Cyl_1 = (Force_Gas.*tand(phi).*Crank_Radius.*cosd(Crank))+(Force_Gas.*Crank_Radius.*sind(Crank));
Torque_Gas_Cyl_2=circshift(Torque_Gas_Cyl_1,-180);
Torque_Gas_Cyl_3=circshift(Torque_Gas_Cyl_1,-360);
Torque_Gas_Cyl_4=circshift(Torque_Gas_Cyl_1,180);

Torque_Gas=Torque_Gas_Cyl_1+Torque_Gas_Cyl_2+Torque_Gas_Cyl_3+Torque_Gas_Cyl_4;

figure;
plot(Crank,Torque_Gas_Cyl_1(:,6))
hold on
plot(Crank,Torque_Gas_Cyl_2(:,6))
plot(Crank,Torque_Gas_Cyl_3(:,6))
plot(Crank,Torque_Gas_Cyl_4(:,6))
plot(Crank,Torque_Gas(:,6))
xlabel('Crank angle (degrees)')
ylabel('Torque (Nm)')
title('Torque of the gas for the individual cylinders @6000rpm')
legend('Cyl 1','Cyl 3','Cyl 2','Cyl 4','Crank')

% Cross Plane engine

Cross_Plane_Torque_Gas_Cyl_2=circshift(Torque_Gas_Cyl_1,-180);
Cross_Plane_Torque_Gas_Cyl_3=circshift(Torque_Gas_Cyl_1,-270);
Cross_Plane_Torque_Gas_Cyl_4=circshift(Torque_Gas_Cyl_1,270);

Cross_Plane_Torque_Gas=Torque_Gas_Cyl_1+Cross_Plane_Torque_Gas_Cyl_2+Cross_Plane_Torque_Gas_Cyl_3+Cross_Plane_Torque_Gas_Cyl_4;

figure;
plot(Crank,Torque_Gas_Cyl_1(:,6))
hold on
plot(Crank,Cross_Plane_Torque_Gas_Cyl_2(:,6))
plot(Crank,Cross_Plane_Torque_Gas_Cyl_3(:,6))
plot(Crank,Cross_Plane_Torque_Gas_Cyl_4(:,6))
plot(Crank,Cross_Plane_Torque_Gas(:,6))
xlabel('Crank angle (degrees)')
ylabel('Torque (Nm)')
title('Torque of the gas for the individual cylinders @6000rpm for a Cross plane engine')
legend('Cyl 1','Cyl 3','Cyl 2','Cyl 4','Crank')
hold off

%% Gas torque for V8 engine
for i =1:length(Crank)
V8_P1500(i)=interp1(RPM_Range,Pres_Data(i,:),1500);
V8_P2500(i)=interp1(RPM_Range,Pres_Data(i,:),2500);
V8_P3500(i)=interp1(RPM_Range,Pres_Data(i,:),3500);
V8_P4500(i)=interp1(RPM_Range,Pres_Data(i,:),4500);
end
V8_P1500=V8_P1500';
V8_P2500=V8_P2500';
V8_P3500=V8_P3500';
V8_P4500=V8_P4500';

V8_Pres_Data=[Pres_Data(:,1),V8_P1500,Pres_Data(:,2),V8_P2500,Pres_Data(:,3),V8_P3500,Pres_Data(:,4),V8_P4500,Pres_Data(:,5)];

V8_Force_Gas = V8_Pres_Data.*1e5.*Area;       % Force of the gas

V8_Torque_Gas_Cyl_1 = (V8_Force_Gas.*tand(V8_phi).*V8_Crank_Radius.*cosd(Crank))+(V8_Force_Gas.*V8_Crank_Radius.*sind(Crank));
V8_Torque_Gas_Cyl_2=circshift(V8_Torque_Gas_Cyl_1,-180);
V8_Torque_Gas_Cyl_3=circshift(V8_Torque_Gas_Cyl_1,-360);
V8_Torque_Gas_Cyl_4=circshift(V8_Torque_Gas_Cyl_1,180);
V8_Torque_Gas_Cyl_5=circshift(V8_Torque_Gas_Cyl_1,0);
V8_Torque_Gas_Cyl_6=circshift(V8_Torque_Gas_Cyl_1,360);
V8_Torque_Gas_Cyl_7=circshift(V8_Torque_Gas_Cyl_1,180);
V8_Torque_Gas_Cyl_8=circshift(V8_Torque_Gas_Cyl_1,-180);

V8_Torque_Gas=V8_Torque_Gas_Cyl_1+V8_Torque_Gas_Cyl_2+V8_Torque_Gas_Cyl_3+V8_Torque_Gas_Cyl_4+V8_Torque_Gas_Cyl_5+V8_Torque_Gas_Cyl_6+V8_Torque_Gas_Cyl_7+V8_Torque_Gas_Cyl_8;

figure;
plot(Crank,V8_Torque_Gas_Cyl_1(:,6))
hold on
plot(Crank,V8_Torque_Gas_Cyl_2(:,6))
plot(Crank,V8_Torque_Gas_Cyl_3(:,6))
plot(Crank,V8_Torque_Gas_Cyl_4(:,6))
plot(Crank,V8_Torque_Gas_Cyl_5(:,6))
plot(Crank,V8_Torque_Gas_Cyl_6(:,6))
plot(Crank,V8_Torque_Gas_Cyl_7(:,6))
plot(Crank,V8_Torque_Gas_Cyl_8(:,6))
plot(Crank,V8_Torque_Gas(:,6),'k--')
xlabel('Crank angle (degrees)')
ylabel('Torque (Nm)')
title('Torque of the gas for the individual cylinders @3500 rpm')
legend('Cyl 1','Cyl 5','Cyl 3','Cyl 7','Cyl 4','Cyl 8','Cyl 2','Cyl 6','Crank')

% Cross Plane engine

V8_Cross_Plane_Torque_Gas_Cyl_2=circshift(V8_Torque_Gas_Cyl_1,-180);
V8_Cross_Plane_Torque_Gas_Cyl_3=circshift(V8_Torque_Gas_Cyl_1,-270);
V8_Cross_Plane_Torque_Gas_Cyl_4=circshift(V8_Torque_Gas_Cyl_1,270);
V8_Cross_Plane_Torque_Gas_Cyl_5=circshift(V8_Torque_Gas_Cyl_1,90);
V8_Cross_Plane_Torque_Gas_Cyl_6=circshift(V8_Torque_Gas_Cyl_1,-90);
V8_Cross_Plane_Torque_Gas_Cyl_7=circshift(V8_Torque_Gas_Cyl_1,180);
V8_Cross_Plane_Torque_Gas_Cyl_8=circshift(V8_Torque_Gas_Cyl_1,360);

V8_Cross_Plane_Torque_Gas=V8_Torque_Gas_Cyl_1+V8_Cross_Plane_Torque_Gas_Cyl_2+V8_Cross_Plane_Torque_Gas_Cyl_3+V8_Cross_Plane_Torque_Gas_Cyl_4+V8_Cross_Plane_Torque_Gas_Cyl_5+V8_Cross_Plane_Torque_Gas_Cyl_6+V8_Cross_Plane_Torque_Gas_Cyl_7+V8_Cross_Plane_Torque_Gas_Cyl_8;

figure;
plot(Crank,V8_Torque_Gas_Cyl_1(:,6))
hold on
plot(Crank,V8_Cross_Plane_Torque_Gas_Cyl_2(:,6))
plot(Crank,V8_Cross_Plane_Torque_Gas_Cyl_3(:,6))
plot(Crank,V8_Cross_Plane_Torque_Gas_Cyl_4(:,6))
plot(Crank,V8_Cross_Plane_Torque_Gas_Cyl_5(:,6))
plot(Crank,V8_Cross_Plane_Torque_Gas_Cyl_6(:,6))
plot(Crank,V8_Cross_Plane_Torque_Gas_Cyl_7(:,6))
plot(Crank,V8_Cross_Plane_Torque_Gas_Cyl_8(:,6))
plot(Crank,V8_Cross_Plane_Torque_Gas(:,6),'k--')
xlabel('Crank angle (degrees)')
ylabel('Torque (Nm)')
title('Torque of the gas for the individual cylinders @3500 rpm for a Cross plane engine')
legend('Cyl 1','Cyl 5','Cyl 3','Cyl 7','Cyl 4','Cyl 8','Cyl 2','Cyl 6','Crank')

%% Cross plane against standard for inline 4 and V8

figure;
plot(Crank,Torque_Gas(:,6))
hold on
plot(Crank,Cross_Plane_Torque_Gas(:,6))
plot(Crank,V8_Torque_Gas(:,6))
plot(Crank,V8_Cross_Plane_Torque_Gas(:,6),'--')
legend('Standard engine','Cross plane engine','V8 Standard engine','V8 Cross plane engine')
title('Standard vs Cross plane engine setup Torque gas @6000rpm of the crank for the inline 4 engines and 3000 rpm for the V8 engines')
xlabel('Crank angle (degrees)')
ylabel('Torque gas (Nm)')



%% Engine Torque

Crank_Torque= Torque_Gas+Torque_Inertia;
Cross_Plane_Crank_Torque= Cross_Plane_Torque_Gas+Cross_Plane_Torque_Inertia;
V8_Crank_Torque=V8_Torque_Gas+V8_Torque_Inertia;
V8_Cross_Plane_Crank_Torque= V8_Cross_Plane_Torque_Gas+V8_Cross_Plane_Torque_Inertia;

figure;
plot(Crank,Crank_Torque(:,6))
hold on
plot(Crank,Cross_Plane_Crank_Torque(:,6))
xlabel('Crank angle (degrees)')
ylabel('Torque (Nm)')
title('Torque of the crank @6000rpm')
legend('Standard engine','Cross plane engine')
hold off

figure;
plot(Crank,Crank_Torque(:,6),'--')
hold on
plot(Crank,Cross_Plane_Crank_Torque(:,6),'--')
plot(Crank,V8_Crank_Torque(:,6))
plot(Crank,V8_Cross_Plane_Crank_Torque(:,6))
xlabel('Crank angle (degrees)')
ylabel('Torque (Nm)')
title('Torque of the crank')
legend('Standard engine','Cross plane engine','Standard engine V8','Cross plane engine V8')
hold off

Torque_average=sum(Crank_Torque)./720;
Cross_Plane_Torque_average=sum(Cross_Plane_Crank_Torque)./720;
V8_Torque_average=sum(V8_Crank_Torque)./720;
V8_Cross_Plane_Torque_average=sum(V8_Cross_Plane_Crank_Torque)./720;

Num_of_Cyl=4;      %Number of cylinders
V8_Num_of_Cyl=8;

Engine_Torque = Torque_average.*Num_of_Cyl;
Cross_Plane_Engine_Torque = Cross_Plane_Torque_average.*Num_of_Cyl;
V8_Engine_Torque = V8_Torque_average.*V8_Num_of_Cyl;
V8_Cross_Plane_Engine_Torque = V8_Torque_average.*V8_Num_of_Cyl;

for i=1:length(RPM_Range)
Max_Torque(i)=max(Crank_Torque(:,i));
Cross_Plane_Max_Torque(i)=max(Cross_Plane_Crank_Torque(:,i));
V8_Max_Torque(i)=max(V8_Crank_Torque(:,i));
V8_Cross_Plane_Max_Torque(i)=max(V8_Cross_Plane_Crank_Torque(:,i));
end

figure;
plot(RPM_Range,Engine_Torque)
hold on
plot(RPM_Range,Torque_average)
plot(RPM_Range,Max_Torque)
plot(RPM_Range,Cross_Plane_Max_Torque)
% Engine torque and Torque average are the same for Cross plane and standard

ylabel('Torque (Nm)')
xlabel('RPM')
legend('engine torque','average torque','Max torque Standard engine at the crank','Max torque Cross plane engine at the crank')
title('Torque for the engine at varying RPMs')
hold off

figure;
hold on
plot(V8_RPM_Range,V8_Engine_Torque)
plot(V8_RPM_Range,V8_Torque_average)
plot(V8_RPM_Range,V8_Max_Torque)
plot(V8_RPM_Range,V8_Cross_Plane_Max_Torque)

ylabel('Torque (Nm)')
xlabel('RPM')
legend('engine torque','average torque','Max torque Standard V8 engine at the crank','Max torque Cross plane V8 engine at the crank')
title('Torque for the engine at varying RPMs on V8 engines')
hold off

%% Work loss due to piston ring

Piston_Ring_Thick=1./1000;      %m Pisotn ring thickness

Ring_Force=Pres_Data.*1e5.*pi*Bore.*Piston_Ring_Thick;

Ring_Friction=Ring_Force.*Coeff_Friction;

for i=1:719
    Work_Ring_Friction(i,1:9)=((Ring_Friction(i,1:9)+Ring_Friction(i+1,1:9))./2).*(S(i+1)-S(i));
end

Ring_Friction_Total=sum(Work_Ring_Friction);

for i=1:719
    Work_Side_Friction(i,1:9)=((Friction_Force_Compressive(i,1:9)+Friction_Force_Compressive(i+1,1:9))./2).*(S(i+1)-S(i));
end

Work_Side_Friction_Total=sum(Work_Side_Friction);

Work_Loss_Total=Ring_Friction_Total+Work_Side_Friction_Total;

%% Power

% Work_Torque=sum(Crank_Torque);
Work_Torque=Torque_average;

Brake_Torque=Work_Torque-Work_Loss_Total;

Indicated_Power=Work_Torque.*Ang_vel;

Friction_Power=Work_Loss_Total.*Ang_vel;

Brake_Power=Indicated_Power-Friction_Power;

% Ford data https://www.fordownersclub.com/forums/topic/97057-are-modern-cars-really-that-fast/
% Actual_Brake_Torque=[108,138,]

figure;
hold on
title('Power and torque at different RPMs')
yyaxis left
plot(RPM_Range,Brake_Torque)
ylabel('Torque (Nm)')
yyaxis right
plot(RPM_Range,Brake_Power./1000)
ylabel('Power (Kw)')
xlabel('RPM')
legend('Torque','Power')

%Cross plane

% Work_Torque=sum(Crank_Torque);
Cross_Plane_Work_Torque=Cross_Plane_Torque_average;

Cross_Plane_Brake_Torque=Cross_Plane_Work_Torque-Work_Loss_Total;

Cross_Plane_Indicated_Power=Cross_Plane_Work_Torque.*Ang_vel;

Cross_Plane_Friction_Power=Work_Loss_Total.*Ang_vel;

Cross_Plane_Brake_Power=Cross_Plane_Indicated_Power-Cross_Plane_Friction_Power;

% Ford data https://www.fordownersclub.com/forums/topic/97057-are-modern-cars-really-that-fast/
% Actual_Brake_Torque=[108,138,]

figure;
hold on
title('Power and torque at different RPMs for cross plane engine')
yyaxis left
plot(RPM_Range,Cross_Plane_Brake_Torque)
ylabel('Torque (Nm)')
yyaxis right
plot(RPM_Range,Cross_Plane_Brake_Power./1000)
ylabel('Power (Kw)')
xlabel('RPM')
legend('Torque','Power')


%%  Effective pressure

Work_Gas=sum(Torque_Gas);

IMEP=Work_Gas./(Crank_Radius.*2.*Area);     % indicated mean effective pressure

FMEP=Work_Loss_Total./(Crank_Radius.*2.*Area);      %Frictional mean effective pressure

BMEP = IMEP-FMEP;       % bRAKE MEAN EFFECTIVE PRESSURE

%% Stresses Concentration on Con Rod
% 
% Rot_Mass=(2/3.*Rod_mass)+Pin_mass;
% Horiz_Acc_Crank_Pin=Crank_Radius.*(Ang_vel.^2).*cosd(Crank);   %Chat gpt says its the horizontal acceleration of the Crank pin
% Rod_Force_Whip=Rot_Mass.*Horiz_Acc_Crank_Pin;
% Rod_Force_Inertia=Recip_Mass.*Ang_Accel_Con_Rod;
% Rod_Force_Total=Rod_Force_Inertia+Rod_Force_Whip;
% 
% Rod_Tau=((Force_Net./2)+Rod_Force_Total)/CSA_Con_Rod;
% 
% % Rod_Stress_Result=sqrt(Stress);
% 
% %% Stress on the bolts
% 
% UTS_Bolt=1100.*1e6;      %Pa
% Yield_Stress=0.9.*UTS_Bolt;
% Axial_Stress=0.6.*Yield_Stress;
% CSA_Bolt=42.43.*1e-6;        %m^2
% CSA_Cap=210.*1e-6;      %m^2
% Bolt_Cap_Diameter=12.9/1000;     %m
% Bolt_Diameter=(sqrt(CSA_Bolt.*4./pi));     %m
% Bolt_Effective_Radius=0.5.*(Bolt_Diameter+Bolt_Cap_Diameter)./2;        %m
% Bolt_Coefficient_Force=0.12;
% Bolt_Coefficient_Torque=0.18;
% Mass_Rotational=((2/3).*Rod_mass);
% Length1=25/1000;        %Length of Rod to Cap intersection  m
% Length2=27.26/1000;        %Length end of thread to Cap intersection  m
% Bolt_Modulus_Of_Elasticity=210000.*1e6;        %N/m^2
% Bolt_Cap_Modulus_Of_Elasticity=170000.*1e6;        %N/m^2
% Bolt_Inner_Thread_Diameter=6.318/1000;     %m
% Bolt_End_Lim=0.45.*UTS_Bolt;      %Pa
% 
% Axial_Force=Axial_Stress.*CSA_Bolt;
% 
% Bolt_Torque_Friction=Bolt_Coefficient_Torque.*Axial_Force.*Bolt_Effective_Radius;
% Bolt_Force_Friction=Bolt_Coefficient_Force.*Axial_Force;
% Bolt_Torque_Friction_Thread=Bolt_Effective_Radius.*Bolt_Force_Friction;
% 
% Total_Torque=Bolt_Torque_Friction+Bolt_Torque_Friction_Thread;
% 
% Bolt_Force_Inertia=(Recip_Mass.*(Crank_Radius.*(Ang_vel.^2).*cosd(Crank))+((Crank_Radius.^2)./Rod).*(Ang_vel.^2).*cosd(2*Crank))+(Mass_Rotational.*(Crank_Radius.*(Ang_vel.^2).*cosd(2*Crank)))./4;
% 
% %Finding Force of extra bolt
% 
% Bolt_Cap_Stiffness=Bolt_Cap_Modulus_Of_Elasticity.*(CSA_Cap./Length1);
% Bolt_Stiffness=Bolt_Modulus_Of_Elasticity.*(CSA_Bolt./Length2);
% Bolt_Force_Extra=(Bolt_Force_Inertia.*Bolt_Stiffness)./(Bolt_Stiffness+Bolt_Cap_Stiffness);
% 
% Bolt_Force_New=Axial_Force+Bolt_Force_Extra;
% 
% %Stresses of bolt
% 
% Bolt_Stress_Min=min(Axial_Force./((pi.*((Bolt_Diameter-(Bolt_Diameter-Bolt_Inner_Thread_Diameter)).^2))./4));
% Bolt_Stress_Min=repmat(Bolt_Stress_Min,1,6);
% Bolt_Stress_Max=max(Bolt_Force_New./((pi.*((Bolt_Diameter-(Bolt_Diameter-Bolt_Inner_Thread_Diameter)).^2))./4));
% Bolt_Stress_Mean=(Bolt_Stress_Max+Bolt_Stress_Min)/2;
% Bolt_Stress_Range=Bolt_Stress_Max-Bolt_Stress_Min;
% Bolt_Stress_Allowable=2.*Bolt_End_Lim.*(1-(Bolt_Stress_Mean./UTS_Bolt));
% Bolt_FSF=Bolt_Stress_Allowable./Bolt_Stress_Range;
% 
% %GOODMANS Diagram
% 
% % Diagram for no surface finish
% 
% %Upper Endurance limit plot
% Y=[Bolt_End_Lim-(UTS_Bolt-Bolt_End_Lim),Bolt_End_Lim,UTS_Bolt];
% X=[-UTS_Bolt,0,UTS_Bolt];
% 
% %Lower Endurance limit plot
% X1=[-UTS_Bolt,0,UTS_Bolt];
% Y1=[-Bolt_End_Lim-(UTS_Bolt+Bolt_End_Lim),-Bolt_End_Lim,UTS_Bolt];
% 
% figure;
% hold on
% plot(X,Y)
% plot(X1,Y1)
% 
% for i=1:length(RPM_Range)
% % Actual range plot
% X2(1,i)=Bolt_Stress_Mean(i);
% X2(2,i)=Bolt_Stress_Mean(i);
% Y2(1,i)=Bolt_Stress_Min(i);
% Y2(2,i)=Bolt_Stress_Max(i);
% 
% %Allowable range plot
% X3(1,i)=Bolt_Stress_Mean(i);
% X3(2,i)=Bolt_Stress_Mean(i);
% Y3(1,i)=Bolt_Stress_Mean(i)-(Bolt_Stress_Range(i)./2.*Bolt_FSF(i));
% Y3(2,i)=Bolt_Stress_Mean(i)+(Bolt_Stress_Range(i)./2.*Bolt_FSF(i));
% end
% 
% plot(X3,Y3, 'y','LineWidth', 5)
% plot(X2,Y2, 'LineWidth',2.5)
% %UTS line
% Y5=[UTS_Bolt,-Bolt_End_Lim-(UTS_Bolt+Bolt_End_Lim)];
% X5=[UTS_Bolt,UTS_Bolt];
% Y6=[UTS_Bolt,UTS_Bolt];
% X6=[UTS_Bolt,-UTS_Bolt];
% plot(X5,Y5,'k--')
% plot(X6,Y6,'k--')
% 
% title('Goodmans diagram for no surface finish on the Bolt')
% legend('Upper endurance limit','Lower endurance limit', 'Allowable range','','','','','','Actual range @1000 rpm','Actual range @2000 rpm','Actual range @3000 rpm','Actual range @4000 rpm','Actual range @5000 rpm','Actual range @6000 rpm','UTS','location','best')
% xlabel('Mean Stress (N/m)')
% ylabel('Stress (N/m)')
% Xlowlim=(UTS_Bolt);
% xuplim=(-UTS_Bolt./2);
% xlim([xuplim Xlowlim ])
% ylowlim=(UTS_Bolt);
% Yuplim=(-UTS_Bolt);
% ylim([Yuplim ylowlim])
% hold off
% 
% % The addition of stress concentrations
% 
% h=(Bolt_Diameter-Bolt_Inner_Thread_Diameter)./2;      %m
% r=0.108/1000;           %m
% h_r=h/r;
% 
% C1=0.833+(2.069.*sqrt(h_r))-(0.009.*h_r);
% C2=2.732-(4.157.*sqrt(h_r))+(0.176.*h_r);
% C3=-8.859+(5.327*sqrt(h_r))-(0.320.*h_r);
% C4=6.294-(3.239.*sqrt(h_r))+(0.154.*h_r);
% Kt=3.065-(3.370.*(2.*h./Bolt_Diameter))+(0.647.*((2.*h./Bolt_Diameter).^2))+(0.658.*((2.*h./Bolt_Diameter).^3));
% 
% Bolt_Stress_Allowable_Kt=2.*(Bolt_End_Lim./Kt).*(1-(Bolt_Stress_Mean./UTS_Bolt));
% 
% Bolt_FSF_Kt=Bolt_Stress_Allowable_Kt./Bolt_Stress_Range;
% 
% %GOODMANS Diagram
% 
% % Diagram for no surface finish
% 
% %Upper Endurance limit plot
% Y=[Bolt_End_Lim-(UTS_Bolt-Bolt_End_Lim),Bolt_End_Lim,UTS_Bolt];
% X=[-UTS_Bolt,0,UTS_Bolt];
% 
% %Lower Endurance limit plot
% X1=[-UTS_Bolt,0,UTS_Bolt];
% Y1=[-Bolt_End_Lim-(UTS_Bolt+Bolt_End_Lim),-Bolt_End_Lim,UTS_Bolt];
% 
% figure;
% hold on
% plot(X,Y)
% plot(X1,Y1)
% 
% 
% for i=1:length(RPM_Range)
% % Actual range plot
% X2(1,i)=Bolt_Stress_Mean(i);
% X2(2,i)=Bolt_Stress_Mean(i);
% Y2(1,i)=Bolt_Stress_Min(i);
% Y2(2,i)=Bolt_Stress_Max(i);
% 
% %Allowable range plot
% X3(1,i)=Bolt_Stress_Mean(i);
% X3(2,i)=Bolt_Stress_Mean(i);
% Y3(1,i)=Bolt_Stress_Mean(i)-(Bolt_Stress_Range(i)./2.*Bolt_FSF_Kt(i));
% Y3(2,i)=Bolt_Stress_Mean(i)+(Bolt_Stress_Range(i)./2.*Bolt_FSF_Kt(i));
% end
% 
% plot(X3,Y3, 'y','LineWidth', 5)
% plot(X2,Y2, 'LineWidth',2.5)
% %UTS line
% Y5=[UTS_Bolt,-Bolt_End_Lim-(UTS_Bolt+Bolt_End_Lim)];
% X5=[UTS_Bolt,UTS_Bolt];
% Y6=[UTS_Bolt,UTS_Bolt];
% X6=[UTS_Bolt,-UTS_Bolt];
% plot(X5,Y5,'k--')
% plot(X6,Y6,'k--')
% 
% title('Goodmans diagram for no surface finish on the Bolt including stress concentration')
% legend('Upper endurance limit','Lower endurance limit', 'Allowable range','','','','','','Actual range @1000 rpm','Actual range @2000 rpm','Actual range @3000 rpm','Actual range @4000 rpm','Actual range @5000 rpm','Actual range @6000 rpm','UTS','location','best')
% xlabel('Mean Stress (N/m)')
% ylabel('Stress (N/m)')
% Xlowlim=(UTS_Bolt);
% xuplim=(-UTS_Bolt./2);
% xlim([xuplim Xlowlim ])
% ylowlim=(UTS_Bolt);
% Yuplim=(-UTS_Bolt);
% ylim([Yuplim ylowlim])
% hold off
% 
% figure;
% plot(RPM_Range,Bolt_FSF)
% hold on
% plot(RPM_Range,Bolt_FSF_Kt)