import numpy as np
from scipy import constants
import matplotlib.pyplot as plt

# Konstanten und Physik

gasses = { 
    "air": {
        "c_p": 1.005e3, # J/kgK
        "c_v": 0.718e3, # J/kgK
        "m_mol": 28.949e-3 , # kg/mol
    },
}

gas = gasses['air']
temp_r = 273.15+25
r_s = constants.R / gas['m_mol']


# Modellparameter

sample_d = 0.2 #m 
sample_l = 0.3 #m 
sample_vel = 100 #m/s
sample_m = 6 #kg

chamber_p = 20e5 # Pa 

sample_a = (sample_d/2)**2 * np.pi
sample_v = sample_a * sample_l
density = sample_m / sample_v #kg/m³
e_kin = sample_m * sample_vel**2 / 2
print(f"Probenenergie: {e_kin/1000:.0f} kJ")

chamber_v = e_kin/ chamber_p
# chamber_v = 0.05
print(f"Chamber Volume = {chamber_v*1000:.1f} l")

m_g = chamber_v*chamber_p / (r_s * temp_r)
print(f"Gas-Masse = {m_g:.3f} kg = {m_g / (sample_m+m_g) *100:.2f} % der Systemmasse")

# e_kin = (sample_m * sample_vel**2 / 2) + (m_g * sample_vel**2 /6)
# print(f"kinetische System-energie = {e_kin:.0f} J")


## isothermer Ansatz

e_pot = chamber_p * chamber_v
print(f"p V = {e_pot} J")
l_f = chamber_v / sample_a
print(f"l_f = V_c / A = {l_f:.3f} m")
barrel_l = l_f * np.exp(e_kin /  e_pot)-l_f
print(f"Lauflänge = {barrel_l:.3f} m")

# Adiabatischer Ansatz

gamma = gas['c_p']/gas['c_v']
print(f"γ = {gamma}")

l_f = chamber_v / sample_a
const = chamber_v*chamber_p *(l_f**(gamma-1)) / (1-gamma)
L = (e_kin/const + l_f**(1-gamma))**(1/(1-gamma)) - l_f

print(f"Lauflänge L = {L:.3f}")

# Grafiken

def work(x):
    return const *((l_f + x)**(1-gamma) - l_f**(1-gamma))

def vel(x):
    return np.sqrt(2*work(x) /sample_m)

def vol(x):
    return chamber_v + sample_a*x

t_0 = 273.15+25

def temp(x):
    return t_0 * (chamber_v / vol(x))**(gamma-1)

def pres(x):
    return e_pot * (chamber_v / vol(x))**(gamma-1) / chamber_v
    

print(f"L = {L:.3f}")

l_v = 3
print(f"W(x={L:.3f}, p={chamber_p:.0e}) = {work(L):.3f}")
print(f"p(x={L:.3f}, p={chamber_p:.0e}) = {pres(L):.2e}")
print(f"V(x={L:.3f}, p={chamber_p:.0e}) = {vol(L):.2e}")
print(f"ΔT(x={L:.3f}, p={chamber_p:.0e}) = { temp(L):.2f}")

print(e_kin + vol(L)*pres(L))


l_array = np.linspace(0,L,400)
plt.plot(l_array, vel(l_array))
plt.ylabel('Speed [m/s]')
plt.xlabel('Sample Travel [m]')
plt.show()