PV diagram for air standard otto cycle using matlab

The aim of the project isto create a Pressure-Volume plot for an air standard OTTO cycle using MATLAB.

An otto cycle consists of 4 process (isentropic compression,constant volume heat addition,isentropic expansion ans constant volume heat rejection )which occurs in a cycle in one rotation of crankshaft. TO plot the pv diagram initially the state points are stated for the given bore and stroke with the known equations.But when plotting using plot they are connected as straight lines as shown. But the isentropic plots are curves.

then a function is created for fnding the volume and it is imported to our function by specifiying function name and by stating pv^gamma is cionstant and then applying plot command,we get the resultant curves.

The main error which occurs during the program is when specifying the function, we must strictly follow the constraint order.

The code for function is as follows

function[V]=Engine_kinematics(bore,stroke,con_rod,c_r,start_crank,end_crank)
a=stroke/2
R=con_rod/a
v_s=(pi*bore^2*stroke)/4
v_c=v_s/(c_r-1)

theeta=linspace(start_crank,end_crank,180);
term1=0.5*(c_r-1);
term2=R+1-cosd(theeta);
term3=(R^2-sind(theeta).^2).^0.5;
V=(1+term1*(term2-term3)).*v_c

end

And the code for program is

clear all
close all
clc
%Inputs
p1=101325
t1=500
gamma=1.4
t3=2300

%Engine geometric parameters
bore=0.1;
stroke=0.1;
con_rod=0.15;
c_r=12;

%calculating swept and clearence volumes
v_swept=(pi*bore^2*stroke)/4
v_clc=v_swept/(c_r-1)


%state variabless 
v1=v_swept+v_clc
v2=v_clc

%state variables at point 2
%p1.v1^gamma=p2.v2^gamma
p2=p1*(c_r^gamma)
%p1v1/t1=p2v2/t2 ideal gas equation
t2=p2*v2*t1/(p1*v1)

constant_c=p1*v1^gamma;
V_expansion=Engine_kinematics(bore,stroke,con_rod,c_r,180,0);
P_expansion=constant_c./V_expansion.^gamma;

%state variables at point 3
v3=v2
p3=p2*t3/t2

constant_c2=p3*v3^gamma;
V_compression=Engine_kinematics(bore,stroke,con_rod,c_r,0,180);
P_compression=constant_c2./(V_compression.^gamma);


%state varibles ath point 4
v4=v1
p4=p3*(v3/v4)^gamma

figure(1) 
hold on
plot(v1,p1,'*','color','r')
plot([v1 v4],[p1 p4])
plot(V_expansion,P_expansion)
plot(v2,p2,'*','color','r')
plot(v3,p3,'*','color','r')
plot(V_compression,P_compression)
plot([v2 v3],[p2 p3])

plot(v4,p4,'*','color','r')

the result is