Air standard Cycle Challenge

%program to solve an otto cycle, to create a PV diagram & to calculate thermal efficiency of the engine
%inputs given
gamma=1.4;
t3=2000;

%input state variables at point 1
%p1=1atm
p1=101325;
t1=350;

%engine geometric parameters in meter
bore=0.1;
stroke=0.15;
c_rod=0.2;
cr=12;

%calculating clearance and swept volume
v_swept=(pi/4)*bore^2*stroke;
v_clear=v_swept/(cr-1);
v1=v_swept+v_clear;
v2=v_clear;

%state variables at point 2
%by isentropic process equation and ideal gas law equation
p2=p1*cr^gamma;
t2=p2*v2*t1/(p1*v1);

%to incoporate the engine kinematics and get a smooth plot during compression
constant_c=p1*v1^gamma;
%calling the function
v_compression=kinematics(bore,stroke,c_rod,cr,180,0);
p_compression=constant_c./v_compression.^gamma;

%state variables at point 3
%v2 and v3 are equal (constant volume heat addition)
v3=v2;
%by ideal gas law equation
p3=p2*t3/t2;

%to incoporate the engine kinematics and get a smooth plot during expansion
constant_c=p3*v3^gamma;
%calling the function
v_expansion=kinematics(bore,stroke,c_rod,cr,0,180);
p_expansion=constant_c./v_expansion.^gamma;

%state variables at point 4
%v4 and v1 are equal (constant volume heat rejection)
v4=v1;
%by ideal gas law equation
p4=p3*(v3/v4)^gamma;

%to calculate the thermal efficiency of the engine
eta=1-(1/cr^(gamma-1));
eff=eta*100;
efficiency=['Efficiency=',num2str(eff),'%'];
disp(efficiency)

%plotting
figure(1)
hold on
plot(v1,p1,'*','Color','r')
plot(v2,p2,'*','Color','r')
plot(v_compression,p_compression,'Color','r')
plot(v3,p3,'*','Color','r')
plot([v2,v3],[p2,p3],'Color','b')
plot(v4,p4,'*','Color','r')
plot(v_expansion,p_expansion,'Color','r')
plot([v4,v1],[p4,p1],'Color','b')
xlabel('volume')
ylabel('pressure')

%function to incoporate the engine kinematics and get a smooth plot during
%compression and expansion
function[vol]=kinematics(bore,stroke,c_rod,cr,start_crank,end_crank)

%engine geometric calculations
a=stroke/2;
R=c_rod/a;
v_s=(pi/4)*bore^2*stroke;
v_c=v_s/(cr-1);

%angle of rotation of crank
theta=linspace(start_crank,end_crank,100);

%to calculate volume
term1=0.5*(cr-1);
term2=R+1-cosd(theta);
term3=(R^2-sind(theta).^2).^0.5;
vol=(1+term1*(term2-term3)).*v_c;
end