P-V diagram visualizer of Otto cycle in PYTHON

import math
import matplotlib.pyplot as plt 


cr = 8
bore_radius = 1.2
stroke =6
con_rod = 4
a = stroke/2
R = con_rod /a
theta = math.radians(0)


v_s = (math.pi/4)*bore_radius*bore_radius*stroke
v_c = v_s/(cr - 1)

def theta_range(t_s,t_e,n):
	theta_start = math.radians(t_s)
	theta_end = math.radians(t_e)
	num = n

	theta = []
	for i in range(0,num):
		d = (theta_end - theta_start)/(num - 1)
		theta.append(theta_start + d*i )
	return theta


def vol_range(bore_radius,cr,stroke,con_rod,theta):
	a = stroke/2
	R = con_rod /a
	v_s = (math.pi/4)*bore_radius*bore_radius*stroke
	v_c = v_s/(cr - 1)


	term1 = 0.5*(cr-1)
	term2 = R + 1 - math.cos(theta)
	term3 =( R**2 - pow(( math.sin(theta)),2) )**.5

	vol = ((1 + term1*( term2 - term3))*v_c)
	return vol
#print(vol,v_c,v_s)

#expansion

theta_1 = theta_range(0,180,80)
vol_exp = []
for q in theta_1:
	vol = vol_range(bore_radius,cr,stroke,con_rod,q)
	vol_exp.append(vol)

#compression

theta_2 = theta_range(180,0,80)
vol_com = []
for q in theta_2:
	vol = vol_range(bore_radius,cr,stroke,con_rod,q)
	vol_com.append(vol)
#print(vol_com)
#print(v_s + v_c)

#initial conditions

P1 =1
T1 = 300
gamma = 1.4
const = P1*(pow((v_s + v_c),gamma))
p_com = []

for v in vol_com:
	p_com.append(const/(pow(v,gamma)))
#print(p_com) 

#max condition

p_max = 50
P3 = p_max
V1 = v_c + v_s
V2 = v_c
V3 = V2
const = P3*pow(V3,gamma)
p_exp = []
for v in vol_exp:
	p_exp.append(const/(pow(v,gamma)))



P2 = p_com[len(vol_com)-1]
P4 = p_exp[len(vol_exp)-1]

# THERMAL EFFICIENCY

TE = 1 - 1/cr**(gamma - 1)
print('Thermal efficieny is',TE)


# plooting the graph

plt.plot(vol_com,p_com)
plt.plot([V2,V2],[P2,P3])
plt.plot(vol_exp,p_exp)
plt.plot([V1,V1],[P4,P1])
plt.xlabel('volume')
plt.ylabel('pressure')
plt.title('P-V diagram of otto cycle')
plt.show()