PV diagram of otto cycle

The code is expalined line by line below:

clear all
close all
clc
%Pv diagram of an ottocycle
%Ottocycle consists of adiabatic compression, constant volume heat
%addition, adiabatic expansion, and constant volume heat rejection
%processes.
%imputs
% gamma for air is 1.4
gamma=1.4
% T3 is the temperature after the heat addition process.
T3=2300
%The bore diameter of the cyclinder
bore=0.1
%The stroke length of the cyclinder
stroke=0.1
%The length of the connecting rod
con_rod=0.15
%Compression ratio
cr=12
%Swept volume
vs= (pi/4)*bore^2*stroke;
%Clearance volume
vc=vs/(cr-1);
% Pressure before compression
P1=101325
%Temperature before compression
T1=300
%Volume before the compression process
v1=vs+vc;
%Volume after the compression process
v2=vc;

%Pressure after the compression process
P2=P1*cr^gamma;
%Temperature after the compression process
T2=P2*v2*T1/(P1*v1);
%Calculation of constant for adiabatic compression process
constant_c=P1*v1^gamma;
%Calling of the function engine_kinematics. This wil calculate the volume
% with varying theta value from 180 to 0 degrees
v_compression=engine_kinematics(0.1,stroke,con_rod,cr,180,0);
%Calculation of pressure during the compression
P_compression=constant_c./v_compression.^gamma;
% Pressure at the end of heat addition
P3=P2*T3/T2
%Volume at the end of heat addition is equal to volume before heat addition
v3=v2;
%Volume before the heat rejection is equal to volume after heat rejection
%process
v4=v1
%%Calling of the function engine_kinematics. This wil calculate the volume
% with varying theta value from 0 to 180 degrees
v_expansion=engine_kinematics(bore,stroke,con_rod,cr,0,180);
%Calculation of constant for adiabatic expansion process
constant_c=P3*v3^gamma;
%Calculation of pressure during the expansion
P_expansion=constant_c./v_expansion.^gamma;
%Pressure at the end of expansion process
P4=P3*(v3/v4)^gamma
%Volume at the end of expansion process is equal to volume at the end of
%heat rejection process
v4=v1

%To draw figure
figure(1)
% First all the x axis values, followed by y axis values
plot([v_compression v3 v_expansion v1],[P_compression P3 P_expansion P1])
%Labelling of x and y axis
xlabel({'Volume','(m^3)'})
ylabel({'Pressure','(N/m^2)'})

%Calculation of thermal efficiency of otto cycle
%Thermal efficiency is given by the below formula
Thermalefficiency = 1-1/(cr)^(gamma-1)

% v is the output, engine_kinematics is the function name, bore, stroke etc
% are input variables
function [v] = engine_kinematics(bore,stroke,con_rod,cr,start_crank,end_crank)

a=stroke/2
R=con_rod/a
%calculation of swept volume and the clearance volume
vs= (pi/4)*bore^2*stroke;
vc=vs/(cr-1);
%crank angle ranging from start to end with 100 values in steps
theta=linspace(start_crank,end_crank,100)
% the bigger formula for v is divided in to 3 terms and then written
% together with term1 , term2 and term3
term1=0.5*(cr-1);
term2=R+1-cosd(theta);
term3=(R^2-sind(theta).^2).^0.5;

v=(1+term1*(term2-term3)).*vc;
%the function should have a end
end