engine kinematics otto cycle using Python

This program solve the otto cycle i.e., it finds the value of temperature and pressure at various stages of the cycle and then generates a plot between pressure and volume and it also calculates the efficiency of the engine (opersting on otto cycle).

#otto cycle
# proces(1-2) - isentropic compression
# process(2-3) - constant vol. heat addition
# process(3-4) - isentropic expansion
# process(4-1) - constant vol. heat rejection

# v_swept = swept volume (v1 - v2)
# v_clearance = clearance volume (v2)
# c_r = compression ratio (v1/v2) 

import math
import matplotlib.pyplot as plt

#defining a function for volume trace by the piston in cylinder w.r.t. crank angle
def engine_kinematics(bore,stroke,con_rod,c_r,start_crank,end_crank):
    
    v_s = math.pi/4 * bore**2 * stroke #swpet volume    
    v_c = v_swept/(c_r-1)  #clearance volume
    
    a = stroke/2  #crank pin radius
    R = con_rod/a
    
    
    n_theta = 50 #no. of values for crank angle
    theta = []
    theta_start = math.radians(start_crank)
    theta_end = math.radians(end_crank)
    
    V = []
    for i in range(0,n_theta):
        theta = theta_start + i*(theta_end - theta_start)/(n_theta-1) 
        term1 = 0.5*(c_r-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 
#state point 1
p1 = 101325  # 1 atm
t1 = 500  # Kelvin
gamma = 1.4
t3 = 2300  # Kelvin

#engine geometric parameters
bore = 0.1
stroke = 0.1
con_rod = 0.15
c_r = 12

#calculating v1 ans v2
v_swept = math.pi/4 * pow(bore,2) * stroke
v_clearance = v_swept/(c_r-1)
v1 = v_clearance + v_swept
v2 = v_clearance

#state point 2
#p1*v1^gamma = p2*v2^gamma | p1*v1/t1 = p2*v2/t2
p2 = p1 * pow((v1/v2),gamma)
t2 = (p2*v2*t1)/(p1*v1)

#for compression stroke starting crank angle is 180 and end crank angle is 0
v_compression = engine_kinematics(stroke,bore,con_rod,c_r,180,0)

constant1 = p1*pow(v1,gamma)

p_compression = []
for v in v_compression:
    p_compression.append(constant1/pow(v,gamma))

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

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

#for expansion stroke starting crank angle is 0 and end crank is 180
v_expansion = engine_kinematics(stroke,bore,con_rod,c_r,0,180)

constant2 = p3*pow(v3,gamma)

p_expansion = []
for v in v_expansion:
    p_expansion.append(constant2/pow(v,gamma))
    

\'\'\'
   thermal efficiency n_th of the otto cycle (engine)
   efficiency = fraction of heat converted to work 
   n_th = W/Qh
   n_th = 1 - (Ql/Qh)
\'\'\'
n_th = 1 - 1/(pow(c_r,(gamma-1)))
print(\'thermal efficiency of the engine is: \' + str(n_th*100) + \' %\')

plt.plot(v1,p1,\'*\')
plt.plot(v2,p2,\'*\')
plt.plot(v3,p3,\'*\')
plt.plot(v4,p4,\'*\')
plt.plot([v2,v3],[p2,p3])
plt.plot([v4,v1],[p4,p1])
plt.plot(v_compression,p_compression)
plt.plot(v_expansion,p_expansion)
plt.xlabel(\'Volume\')
plt.ylabel(\'Pressure\')
plt.show()

output:-

thermal efficiency of the engine is: 62.98928275128466 %