import numpy as np
import matplotlib.pyplot as plt

l_barrel = 2 #m
d_barrel = 0.2 #m

l_sabot = 0.03 #m
mu_sabot = 0.04

m_payload = 6 #kg

vol_chamber = 0.015 #m^3
p_chamber = 20e5 #Pa
t_chamber = 300 #K
gamma = 1.4

dt = 0.00003

a_b = d_barrel**2 * np.pi / 4

def step(t, f, a, v, l, vol, p, temp):
    global dt
    ft = a_b * p - mu_sabot*d_barrel*np.pi*l_sabot*p
    at = ft / m_payload
    vt = v + at*dt
    lt = l + vt*dt
    volt = vol + vt*dt*a_b
    pt = p * (vol/volt)**gamma
    tt = t+dt
    tempt = t_chamber * (vol_chamber/vol)**(gamma-1)
    return tt, ft, at, vt, lt, volt, pt, tempt

data = [(
    0, 0, 0, 0, 0, vol_chamber, p_chamber, t_chamber,
),]

i = 0
print(data[0])
while data[-1][4] < l_barrel:
    i+=1
    print(i)
    erg = step(*data[-1])
    data.append(erg)

print("done")

tt, ft, at, vt, lt, volt, pt, tempt = np.transpose(data)

plt.plot(lt, pt/1000)
plt.xlabel("Weg ([m]")
plt.ylabel("Pressure [bar]")
plt.show()