Otto Cycle - PV Diagram Efficiency

######comments#######

You havent attached the function file

#####################################

% Otto Cycle

clear all

close all

clc

% Inputs

gamma = 1.4;

t3 = 2300;

% State Variables

p1 = 101325;

t1 = 500;

% Engine Geometric Parameters

bore = 0.1;

stroke = 0.1;

con_rod = 0.15;

cr = 12;

% Calculating the 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 Variables at State Point 2

% p2v2^gamma = p1v1^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 Variables at State Point 3

v3 = v2;

% p3v3/t3 = p2v2/t2 | p3 = p2*t3/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 Variables at State Point 4

v4 = v1;

% p3v3^gamma = p4v4^gamma | p4 = p3(v3/v4)^gamma

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

figure(1)

hold on

plot(v1,p1,\'*\',\'color\',\'r\')

plot(v_compression,p_compression)

plot(v2,p2,\'*\',\'color\',\'r\')

plot([v2,v3],[p2,p3],\'color\',\'b\')

plot(v3,p3,\'*\',\'color\',\'r\')

plot(v_expansion,p_expansion)

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

plot([v4,v1],[p4,p1],\'color\',\'b\')

xlabel(\'Volume\')

ylabel(\'Pressure\')

grid(\'on\')

% Thermal Efficiency Calculation

t_eff = (1 - (1/(cr^(gamma-1))))*100;

%t_eff = (1-(v2/v1)^(gamma-1))*100; (Tried using this equation too)

fprintf(\'Thermal Efficiency of the given Otto Cycle is %f percentage\\n\',t_eff);

% Engine Kinematics Function

function [v] = Engine_Kinematics(bore, stroke, con_rod, cr, start_crank, end_crank)

a = stroke/2;
r = con_rod/a;

v_swept = pi/4*bore^2*stroke;
v_compression = v_swept/(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_compression;

end