Otto cycle : PV diagram visualizer in MATLAB

%PV Diagram visualizer for otto cycle

clear all
close all
clc

%input
gamma=1.4;
t3 = 2500;

%state variables
p1=101300;
t1=600;

%engine geometric parameters
bore = 0.1;
stroke = 0.1;
con_rod = 0.15;
cr = 12;

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

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

%p1v1/t1=p2v2/t2 | t2= p2*v2*t1/(p1*v1)
t2= p2*v2*t1/(p1*v1);

constant_c = p1*v1^gamma;

V_compression= engine_kinematics (bore,stroke,con_rod,cr,180,0);

P_compression = constant_c ./ V_compression.^gamma;

%state varriable at point 3
v3=v2;

%p3v3/t3=p2v2/t2 | p3 = p2*v2*t3/(v3*t2)
p3 = p2*t3/t2;

constant_c = p3*v3^gamma;

V_expansion= engine_kinematics (bore,stroke,con_rod,cr,0,180);

P_expansion = constant_c ./ V_expansion.^gamma;

% state varriable at point 4
v4=v1;

%p3v3^gamma=p4v4^gamma
p4= p3*(v1/v4);

%p4v4^gamma = p3*v3^gamma;
p4= p3*(v3/v4)^gamma;

%plotting
figure(1)
hold on
plot (v1,p1,'*','color','r');
plot (V_compression , P_compression );
plot (v2,p2,'*','color','r');
plot (v3,p3,'*','color','r');
plot (V_expansion , P_expansion);
plot (v4,p4,'*','color','r');
plot ([v2 v3],[p2 p3])
plot ([v4 v1],[p4 p1])
xlabel('Volume')
ylabel('Pressure')
grid on

%efficiency of otto cycle
eff = 1-(1/(cr^(gamma-1)))
txt1 = strcat ('Efficiency =',num2str(eff));
txt2 = strcat ('Compression ratio = ',num2str(cr));
text(v1/2 , p3/1.3,txt2)
text(v1/2,p3/1.5,txt1)

%Function for volume trace in process 1-2 and process3-4

function [V] = engine_kinematics (bore,stroke,con_rod,cr,start_crank,end_crank)

a=stroke/2;
R=con_rod/a;

v_s = pi/4*bore^2*stroke;

v_c = v_s/(cr-1);

theta = linspace(start_crank,end_crank,100);

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