Week 2 Air standard Cycle

""" 
Otto Cycle Simulator
 """
import math 
import matplotlib.pyplot as plt
def engine_kinematics (bore,stroke,conn_rod,cr,start_crank,end_crank):
 """
 Engine Kinematics
 """
 # Geometric paramters
 a=stroke/2
 R=conn_rod/a
 # Volume parameters
 V_s = math.pi*(1/4)*pow(bore,2)*stroke
 V_c = V_s/(cr-1)
 sc=math.radians(start_crank)
 ec=math.radians(end_crank)
 num_values=50
 dtheta=(ec-sc)/(num_values-1)
 V=[]
 for i in range(0,num_values):
  theta=sc+i*dtheta
  term1=0.5*(cr-1)
  term2=R+1-math.cos(theta)
  term3=pow(R,2)-pow(math.sin(theta),2)
  term3=pow(term3,0.5)
  V.append((1+term1*(term2-term3))*V_c)
 return V
  # inputs
p1=101325
t1=500
gamma=1.4
t3=2300
# volume parametrs
bore=0.1
stroke=0.1
conn_rod=0.15
cr=8
#volume calculation
V_s=(math.pi/4)*pow(bore,2)*stroke
V_c=V_s/(cr-1)
v1=V_s+V_c
# state 2
v2=V_c
#p2v2^gamma = p1v1^gamma
p2=p1*pow(v1,gamma)/pow(v2,gamma)
# p2v2/t2 = p1v1/t1 | x=p1v1/t1 | x=p2v2/t2 | t2=p2v2/x
V_compression =engine_kinematics(bore,stroke,conn_rod,cr,180,0)
constant=p1*pow(v1,gamma)
P_compression=[]
for v in V_compression:
	P_compression.append(constant/pow(v,gamma))
	print(P_compression)
x=p1*v1/t1
t2=p2*v2/x
# state 3
v3=v2
# p3v3/t3=p2v2/t2 | p2v2/t2=x | p3=xt3/v3
x=p2*v2/t2
p3=x*t3/v3
V_expansion =engine_kinematics(bore,stroke,conn_rod,cr,0,180)
constant=p3*pow(v3,gamma)
P_expansion=[]
for v in V_expansion:
	P_expansion.append(constant/pow(v,gamma))
	print(P_expansion)

# state 4
v4=v1
#p4v4^gamma =p3v3^gamma
p4=p3*pow(v3,gamma)/pow(v4,gamma)
#p4v4/t4=x
t4=p4*v4/x
plt.plot(V_compression,P_compression,label="Adiabatic Compression")
plt.plot([v2,v3],[p2,p3],label="Constant Volume heat addition")
plt.plot(V_expansion,P_expansion,label="Adiabatic Expansion")
plt.plot([v4,v1],[p4,p1],label="Constant Volume heat rejection")
plt.xlabel("Volume")
plt.ylabel("Pressure")
plt.legend()
plt.show()
Thermal_efficiency= 1-pow(cr,(1-gamma))
print(Thermal_efficiency)

Thermal efficincy 56 % and cr =8