PV diagram visualizer

https://drive.google.com/open?id=1FtsVWX4Hw7EEj1qSUK8crOt7iGjiPEz3

clear all
close all
clc

% state variables at point-1

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

%engine geometric parameters
bore=0.1;
stroke=0.1;
con_rod=0.15;
cr=12;

swept_v=(pi/4)*bore^2*stroke;
clear_v= swept_v/(cr-1);

v1=swept_v+clear_v
v2=clear_v

%state variables at point-2

 %p2*v2^gamma=p1v1^gamma
 %p2=p1*(v1/v2)^gamma
 
 p2=p1*(v1/v2)^gamma
 t2=p2*v2*t1/(p1*v1)
 
 constant_c=p1*v1^gamma;
  V_comp=engine_kinematics(bore,stroke,con_rod,cr,180,0);
  p_comp=constant_c./V_comp.^gamma;

% state variables at point 3
 v3=v2
 p3=p2*t3/t2
 constant_c=p3*v3^gamma;
  V_exp=engine_kinematics(bore,stroke,con_rod,cr,0,180);
  p_exp=constant_c./V_exp.^gamma;

% state variables at point 4
 v4=v1
 p4=p3*(v3/v4)^gamma

 figure(1)
 hold on
 plot(v1,p1,\'*\',\'color\',\'r\')
 plot(V_comp,p_comp)
 plot(v2,p2,\'*\',\'color\',\'r\')
 plot([v2 v3],[p2 p3])
 plot(v3,p3,\'*\',\'color\',\'r\')
 plot(V_exp,p_exp)
 plot(v4,p4,\'*\',\'color\',\'r\')
plot([v4 v1],[p4 p1])
 xlabel(\'Volume\')
 ylabel(\'Pressure\')
 title(\' Otto Cycle P-V dig\')



%Piston_kinematics
 
function [V]= engine_kinematics(bore, stroke, con_rod, cr, start_crank, end_crank)
%engine geometric parameters

a=stroke/2;
R=con_rod/a;

swept_v=(pi/4)*bore^2*stroke;
clear_v= swept_v/(cr-1);

theta=linspace(start_crank,end_crank,100);

%(v/clear_v)=1+0.5(cr-1)[R+1-cosd(theta)-R^2-sin(theta)^2])^0.5

term1=0.5*(cr-1);
term2=R+1-cosd(theta);
term3=(R^2-sind(theta).^2).^0.5;

V=(1+term1*(term2-term3)).*clear_v;
 


end