Air standard cycle

%Air standard cycle
clear all
close all
clc

bore=0.1
stroke=0.1
cond_rod=0.15
cr=12

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

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

p2=p1*cr^gamma;

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

constant_c=p1*v1^gamma;
V_compression= otto_cycle(bore,stroke,cond_rod,cr,180,0);
P_compression= constant_c./V_compression.^gamma;

v3=v2;
p3=p2*t3/t2;

constant_c=p3*v3^gamma;
V_expansion= otto_cycle(bore,stroke,cond_rod,cr,0,180);
P_expansion= constant_c./V_expansion.^gamma;

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

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

xlabel(\'Volume\');
ylabel(\'Pressure\');

grid on

%Efficiency
efficiency=(1-(1/cr^(gamma-1)))*100
disp(efficiency)

%figure(2)
%plot([v1 v2 v3 v4],[p1 p2 p3 p4]);

function V=otto_cycle(bore,stroke,cond_rod,cr,start_crank,end_crank);

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

theta=linspace(start_crank,end_crank,180);
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