Otto Cycle Python

import math
import matplotlib.pyplot as plt

def engine_kinematics(bore,stroke,cond_rod,cr,start_crank,end_crank):

	

	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(pow(R,2)- pow(math.sin(theta),2),0.5);
		V.append((1+term1*(term2- term3))*v_c)

	return V


bore=0.1
stroke=0.1
cond_rod=0.15
cr=12

p1= 101325
t1=500
gamma=1.4;
t3=2300;
a=stroke/2;
R=cond_rod/a;

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

#point 2
v2=v_c;


#pressure
p2=p1*pow(v1,gamma)/pow(v2,gamma);

#p2v2/t2=p1v1/t1
t2=p2*v2*t1/(p1*v1);
rhs=p2*v2/t2

V_compression= engine_kinematics(bore,stroke,cond_rod,cr,180,0)

constant = p1*pow(v1,gamma)
P_compression= []

for v in V_compression:
	P_compression.append(constant/pow(v,gamma))


#point 3
v3=v2;
p3=p2*t3/t2;

V_expansion= engine_kinematics(bore,stroke,cond_rod,cr,0,180)

constant = p3*pow(v3,gamma)
P_expansion= []

for v in V_expansion:
	P_expansion.append(constant/pow(v,gamma))

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

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()

efficiency=(1-(1/pow(cr,(gamma-1))))*100
print(efficiency)