To calculate the efficency of an otto cycle using matlab

To calculate the efficency and plot the otto cycle:-

function v = engine_kinematics(bore,stroke,cond_rod,cr, start_crank,end_crank)



v_s = (pi/4)*bore^2*stroke;
v_c = v_s/(cr-1);
a= stroke/2;
R= cond_rod/a;


theta = linspace(start_crank,end_crank,180);

term1 = 0.5*(cr-1);
term2 = R +1- cosd(theta);
term3= (R^2- sind(theta).^2).^0.5;


v = (1+term1*(term2-term3))*v_c;




end 




  

close all
clear all
clc

%inputs 


p1= 101325;
t1= 500;
t3= 2300;
gamma= 1.4;

%engine geometric properties

bore = 0.1;
stroke= 0.1;
con_rod= 0.15;
cr = 12;

%calculating the clearnce and swept volume

v_swept = (pi/4)*(bore^2)*stroke;
v_clearance = v_swept/(cr-1);


v1= v_swept+ v_clearance;
v2= v_clearance;

c1= p1*v1^gamma;
v_compression= engine_kinematics(bore, stroke,con_rod,cr,180,0)
p_compression= c1./(v_compression.^gamma)

% state variables at point 2
% p2*(v2)^gamma = p1*(v1)^gamma
p2= p1*(cr)^gamma


%calculating t2
%p1*v1/t1=p2*v2/t2

t2 = (p2/p1)*(v2/v1)*t1


%stae variable at point 3 
v3 = v2;

%calculating tempertaure at point 3 
%p3/t3=p2/t2

p3 = p2*(t3/t2);

c2 = p3*v3^gamma;
v_expansion = engine_kinematics(bore, stroke,con_rod,cr,0,180)
p_expansion= c2./(v_expansion.^gamma)

%state variables at point 4

v4 = v1

%p3v3^gamma =p4*v4^gamma

p4 = p3*(v3/v4)^gamma

%calculating temperature at point 4
t4 = (p4*v4*t3)/(p3*v3);

hold on
plot(v1, p1, '*', 'color', 'g');
plot(v_compression, p_compression );
plot(v2,p2,'*','color','g');
plot(v3,p3,'*','color','g');
plot(v_expansion,p_expansion)
plot(v4,p4,'*','color','g');



plot([v1 v_compression v3 v_expansion],[p1 p_compression p3 p_expansion], 'color', 'r')

xlabel('volume')
ylabel('pressure')

% calculating the efficency of otto cycle in terms of percentage

eeta1= (1- ((t4-t1)/(t3-t2)))*100