Air standard Cycle

import math
import matplotlib.pyplot as plt

def engine_kinematics(bore,stroke,cr_length,start_cranck,end_cranck):
a=stroke/2
R=cr_length/a
sc=math.radians(start_cranck)
ec=math.radians(end_cranck)

n_value=25
dtheta=(ec-sc)/20
v=[]

for i in range(0,n_value):
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(vc*(1+term1*(term2-term3)))
return v


bore=0.1
stroke=0.1
cr=12
vs=math.pi/4*pow(bore,2)*stroke
vc=vs/(cr-1)
gamma=1.4
cr_length=0.1
t3=2500

t1=500
p1=101325
v1=vs+vc

v2=vc
#p1*v1^gamma=p2*v2^gamma
lhs=p1*pow(v1,gamma)
p2=lhs/(pow(v2,gamma))
#p1v1/t1=p2v2/t2
lhs=p1*v1/t1
t2=p2*v2/lhs

v_compression=engine_kinematics(bore,stroke,cr_length,180,0)

constant=p1*pow(v1,gamma)
p_compression=[]
for v in v_compression:
p_compression.append(constant/(pow(v,gamma)))

v3=v2
#p2v2/t2=p3v3/t3
lhs=p2*v2/t2
p3=lhs*t3/v3

v_expansion=engine_kinematics(bore,stroke,cr_length,0,180)
constant=p3*pow(v3,gamma)
p_expansion=[]
for v in v_expansion:
p_expansion.append(constant/(pow(v,gamma)))

v4=v1
#p3v3^gamma=p4v4^gamma
lhs=p3*pow(v3,gamma)
p4=lhs/(pow(v4,gamma))
#p3v3/t3=p4v4/t4
lhs=p3*v3/t3
t4=p4*v4/lhs
#thermal efficiency=1-(1/(cr^gamma)-1)
ef1=1/(pow(cr,gamma)-1)
efficiency=1-ef1
print(efficiency)

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.xlabel(\'volume\')
plt.ylabel(\'pressure\')
plt.grid(True)
plt.show()

efficiency=0.9681762181323296