Creating a Pressure vs Volume curve of Otto Cycle

Otto Cycle programm and PV diagram

clear all
close all
clc

%otto cycle

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

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

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

v1=v_s+v_c;

constant_c=p1*v1^gamma;
v_comp=engine_kinematics(bore,stroke,cr,con_rod,180,0);
p_comp=constant_c./v_comp.^gamma;

v2=v_c;

p2=p1*cr^gamma;

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

v3=v2;

p3=p2*t3/t2;

constant_c=p3*v3^gamma;
v_exp=engine_kinematics(bore,stroke,cr,con_rod,0,180);
p_exp=constant_c./v_exp.^gamma;

v4=v1;

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

t4=p4*v4*t3/(p3*v3);

%thermal efficiency of otto cycle
%thermal efficiency=work done/heat input

%work done=Cv(t3-t2)-Cv(t4-t1)
%heat input=Cv(t3-t2)
%thermal efficiency=1-((t4-t1)/(t3-t2))

T=1-((t4-t1)/(t3-t2))

%plotting
hold on
plot(v1,p1,\'*\',\'color\',\'r\')
plot(v_comp,p_comp,\'linewidth\',1,\'color\',\'g\')
plot([v2 v3], [p2 p3],\'linewidth\',1,\'color\',\'b\')
plot(v2,p2,\'*\',\'color\',\'r\')
plot(v3,p3,\'*\',\'color\',\'r\')
plot(v_exp,p_exp,\'linewidth\',1,\'color\',\'g\')
plot([v4 v1], [p4 p1],\'linewidth\',1,\'color\',\'b\')
plot(v4,p4,\'*\',\'color\',\'r\')
xlabel(\'VOLUME\')
ylabel(\'PRESSURE\')

Defining Function:-

function [V]=engine_kinematics(bore,stroke,cr,con_rod,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);

%formula of volume change for each angle is V=1+0.5(cr-1)[R+1-cosd(theta)-(R^2-sind(theta)^2)^0.5]

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