Code to animate p-V of Air standard Otto cycle

%enginegeometry%

d=0.1;

L=0.1;

conrodlength=0.15;

crankradius=L/2;

R=conrodlength/crankradius;

rc=12;vs=(pi/4)*d*d*L;

vc=vs/(rc-1);gamma=1.4;

%state 1%

p1=101325;

T1=500;

v1=vs+vc;

%state 2%

v2=vc;

p2=p1*(rc)^gamma;

T2=T1*(rc)^(gamma-1);

theta_comp=linspace(180,360,180);

term1_comp=0.5*(rc-1);

term2_comp=R+1-cosd(theta_comp);

term3_comp=((R^2)-(sind(theta_comp)).^2).^(0.5);

instantvolume_comp=(1+term1_comp*(term2_comp-term3_comp)).*vc;

const_comp=p1*(v1^gamma);

instantpressure_comp=const_comp./(instantvolume_comp).^(gamma);

%state 3%

T3=2500;

v3=vc;p3=p2*(T3/T2);

%state 4%

T4=((T3)/((rc)^(gamma-1)));

p4=p3*(1/rc)^gamma;

v4=v1;

theta_exp=linspace(360,540,180);

term1_exp=0.5*(rc-1);

term2_exp=R+1-cosd(theta_exp);

term3_exp=((R^2)-(sind(theta_exp)).^2).^(0.5);

instantvolume_exp=(1+term1_exp*(term2_exp-term3_exp)).*vc;

const_exp=p3*(v3^gamma);

instantpressure_exp=const_exp./(instantvolume_exp).^(gamma);

i=1;

plot(instantvolume_comp,instantpressure_comp,'color','k');

pause(10);

M(i)=getframe(gcf);

line([v2 v3],[p2 p3],'color','r');

pause(10);

M(i+1)=getframe(gcf);

hold on;

plot(instantvolume_exp,instantpressure_exp,'color','b');

pause(10);

M(i+2)=getframe(gcf);

line([v4 v1],[p4 p1],'color','g');

pause(10);

M(i+3)=getframe(gcf);

legend('Isentropic compression','Isochoric Heat addition','Isentropic expansion','Isochoric Heat rejection');

xlabel('Volume in cu.m');

ylabel('Pressure in Pa');

title('p-V plot of Otto Cycle');

movie(M);

videofile=VideoWriter('ottoanimation.avi','Uncompressed AVI');

open(videofile);

writeVideo(videofile,M);

close(videofile);

thermalefficiency=1-(1/(rc^(gamma-1)))