plotting the pv diagram and obtaining the thermal efficiency

I have first written the funciton and then the main program. The thermal efficiency will be 1-(t1/t3) which comes out to be 78.261%. Code for this is in the main program.

                    % function for piston movement in the cylinder

function[v] = cylinder(bore, stroke, con_rod, cr, start,last);

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

% known parameters
a= stroke/2;
R = con_rod/a;
vs = (pi/4)*bore^2*stroke;
vc = vs/(cr-1);

% value of theta
theta = linspace(start,last,180);

% defining the path of the piston
term1 = 0.5*(cr-1);
term2 = R+1- cosd(theta);
term3 = (R^2- sind(theta).^2).^0.5;
v= (1+term1*(term2- term3))*vc;

end

 

          % plotting pv diagram and calculating thermal efficiency

clear all
close all
clc

% inputs
gamma = 1.4;

% known values
p1 = 101325;
t1 = 500;
t3 = 2300;

% geometric values of cylinder
bore = 0.1;
stroke = 0.1;
con_rod  = 0.15;
cr = 12;

% calculating swept volume of the clyinder
vs = (pi/4)*bore^2*stroke;
vc = vs/(cr-1);

v1 = vs+vc;
v2 = vc;

% calculating values at point 2
p2 = p1*(cr)^gamma;
t2 = v2*p2*t1/(p1*v1);
constant = p1*v1^gamma;
vcompression = cylinder(bore, stroke, con_rod, cr,180, 0);
pcompression = constant./vcompression.^gamma;


% calculating values at point 3
v3 = v2;
p3 = p2*(t3/t2);
constant = p3*v3^gamma;
vexpansion = cylinder(bore, stroke, con_rod, cr,0,180);
pexpansion = constant./vexpansion.^gamma;
%calculating values at point 4
v4 = v1;
p4 = p3*(v3/v4)^gamma;

% calculating thermal efficiency
Te = 1-(t1/t3); 

% plotting the pv diagram
plot(v1,p1,\'*\',\'color\',\'r\');
hold on
plot(vcompression,pcompression,\'color\', \'g\');
plot(v2,p2,\'*\',\'color\',\'r\');
hold on
plot(v3,p3,\'*\',\'color\',\'r\');
plot([v2 v3],[p2 p3], \'color\', \'b\');
hold on 
plot(vexpansion, pexpansion,\'color\', \'m\');
plot(v4, p4,\'*\',\'color\',\'r\');
plot([v4 v1], [p4 p1], \'color\', \'b\');
xlabel(\'volume\')
ylabel(\'pressure\')
legend(\'pv diagram\')