Air Standard Cycle Using Python

"""
Air Standard Cycle
A code that creates a PV diagram and outputs the thermal efficiency of the engine
"""

import matplotlib.pyplot as plt
import math

#inputs
p1=101325
t1=500
t3=2300

gamma=1.4

# Geometric paprameters
bore = 0.1
stroke = 0.1
con_rod = 0.15
cr = 12

a = stroke/2
R = con_rod/a

# Volume
v_s=(math.pi/4)*pow(bore,2)*stroke
v_c=v_s/(cr-1)
v1=v_s+v_c

def engine_kinematics(bore, stroke, con_rod, cr, start_crank, end_crank):
	sc=math.radians(start_crank)
	ec=math.radians(end_crank)

	n=100

	dtheta=(ec-sc)/(n-1)
	v=[]

	for i in range(0,n):
		theta=sc+i*dtheta
		term1=0.5*(cr-1)
		term2=R+1-math.cos(theta)
		term3=pow(pow(R,2)-pow(math.sin(theta),2),0.5)

		v.append((1+term1*(term2-term3))*v_c)
	return v	
		

# State point 2
v2= v_c
p2= p1*pow((v1/v2),gamma)
t2= t1*p2*v2/(p1*v1)

V_compression=engine_kinematics(bore,stroke,con_rod,cr,180,0)
P_compression=[]

for a in V_compression:
	P_compression.append(p1*pow((v1/a),gamma))

#state point 3
v3=v2
p3=p2*v2*t3/(t2*v3)

V_expansion=engine_kinematics(bore,stroke,con_rod,cr,0,180)
P_expansion=[]

for a in V_expansion:
	P_expansion.append(p3*pow((v3/a),gamma))

#state point 4
v4=v1
p4=p3*pow((v3/v4),gamma)
t4=p4*v4*t3/(p3*v3)

# Thermal Efficiency
Therm = 1 - ((t4-t1)/(t3-t2))

plt.plot([v2,v3],[p2,p3])
plt.plot(V_compression,P_compression)
plt.plot(V_expansion,P_expansion)
plt.plot([v4,v1],[p4,p1])
plt.show()
print(Therm)