Air Standard Otto Cycle(PV diagram)

clear all
close all
clc

% inputs
g=1.4;
t3=2300;
% state1 inputs
p1=101325;
t1=25;
% engine desc
   bore=0.1;
   conn_rod=2;
   stroke=0.1;
   c_r=12;
   % calculation of v1
 v_swept=(pi/4)*bore^2*stroke;  
v_clc=v_swept/(c_r-1);
 v1=v_swept+v_clc;

v2=v_clc;
% state2 inputs
p2=p1 *c_r^g;
t2=p2*v2*t1/(p1*v1);
c=p1*v1^g;
v_compression=nk2222(bore,stroke,conn_rod,c_r,180,0);
p_compression=c/v_compression.^g;
% state3 inputs
v3=v2;
p3=p2*t3/t2;
% state4 inputs
v4=v1;
p4=p3*(v3/v4).^g;
c_constant=p3*v3^g;
v_expansion=nk2222(bore,stroke,conn_rod,c_r,0,180);
p_expansion=c_constant./v_expansion.^g;
figure(1)
hold on
plot(v1,p1)
plot(v2,p2)
plot(v_compression,p_compression,'linewidth',3,'color','r')
plot(v3,p3)
plot([v2 v3],[p2 p3],'linewidth',3,'color','b')
plot(v4,p4)
plot(v_expansion,p_expansion,'linewidth',3,'color','r')
plot([v4 v1],[p4 p1],'linewidth',3,'color','b')





function[V]=nk2222(bore,stroke,conn_rod,c_r,start_crank,end_crank)
a=stroke/2;
R=conn_rod/a;
v_s=pi/4*bore.^2*stroke;
v_c=v_s/(c_r-1);
theta=linspace(start_crank,end_crank,100)
term1=0.5*(c_r-1);
term2=R+1-cosd(theta);
term3=(R^2-sind(theta).^2).^0.5;
V=(1+term1*(term2-term3)).*v_c;
end