Otto Cycle

import math
import matplotlib.pyplot as plt
import numpy
def eng_kin(bore, stroke, con_rod_length, CR, start_angle, end_angle):
	
	a = stroke/2
	R = con_rod_length/a
	
	# volume calculation
	
	# swept volume
	V_s = (math.pi/4)*(bore**2)*stroke
	
	# clearance volume
	V_c = V_s/(CR-1)
	
	sc = math.radians(start_angle)
	ec = math.radians(end_angle)
	
	intervals = 30
	
	theta_int = (ec - sc)/(intervals-1)
	
	volume = [] 
	
	for i in range(0, intervals):
		crank_angle = sc + (i*theta_int)
		t1 = 0.5*(CR-1)
		t2 = R + 1 - math.cos(crank_angle)
		t3 = (R**2) - ((math.sin(crank_angle)**2)**0.5)
		
		volume.append((1 + t1*(t2-t3))*V_c)
		
	return volume
	
# inputs

P1 = 101325
T1 = 270
gamma = 1.4
T3 = 1900

bore = 0.1
stroke = 0.1
con_rod_length = 0.15
CR  = 15 


V_s = (math.pi/4)*(bore**2)*stroke
V_c = V_s/(CR-1)

V1 = V_s + V_c

V2 = V_c

P2 = P1*(V1/V2)**gamma

T2 = ((P2*V2)/(P1*V1))*T1

v_comp = eng_kin(bore, stroke, con_rod_length, CR, 180, 0)

const_1 = P1*V1**gamma

p_comp = []

for vol in v_comp:
	p_comp.append(const_1/vol**gamma)
	
V3 = V2 

P3 = (P2*V2/T2)*(T3/V3)

v_exp = eng_kin(bore, stroke, con_rod_length, CR, 0, 180)

const_2 = P3*V3**gamma

p_exp = []

for v in v_exp:
	p_exp.append(const_2/v**gamma)
	
V4 = V1

P4 = P1*V1**gamma/V4**gamma

T4 = (P4*V4*T1)/(P1*V1)

plt.plot([V2,V3],[P2,P3])
plt.plot(v_comp,p_comp)
plt.plot(v_exp,p_exp)
plt.plot([V4,V1],[P4,P1])

plt.xlabel(\'Volume\')
plt.ylabel(\'Pressure\')
plt.show()