# (c) 2015, Joachim Schwender
# Der cw Wert ist als konstant angenommen, was unterhalb der Schallgeschwindigkeit eine gute Näherung ist.
# Die berechnung stimmt nur solange die Raketenbewegung auch eindimensional bleibt, 
# also eine translatorische Bewegung auf einer Geraden ausführt.
# Sobald sich die Rakete über eine Hozizontalachse dreht werden die Abweichungen schnell sehr gross.
# Unstabilisiert ist das nicht machbar. Eine aerodynamische Stabilisierung ist unabdingbar.
#
dt = 0.001;    # Die Schrittweiter für die Berechnung [s]
#------- Raketenabmessungen --------
D_d = 0.021;   # Austrittsdurchmesser der Düse [m]
D_g = 0.08;    # Gehäusedurchmesser der Rakete [m]
m_geh = 0.049; # Gehäusemasse [kg]
vol_rak = 0.00155;  # Füllvolumen der Rakete [m³]
cw = 0.40;     # Luftwiederstandsbeiwert cw, woher man den auch immer bekommt...
#-------- Treibladung -----------
fuellgrad = 0.35;    # Füllgrad der Flasche, also Volumenanteil mit Wasser am Flaschenvolumen (%)
Druck = 800000.0;     # Druck der Gasfüllung beim Start der Rakte [Pa], Hinweis: 100000 Pa = 1 bar
#
# ------ Konstanten ------
dichte_luft = 1.234;   # Dichte der Luft [kg/m³]
dichte_w = 1000.0;      # Dichte des Wassers [kg/m³]
g = 9.81;
k_luft = 1.4;          # isentropenexponent der Luft (trocken, 0°C), zur Berechnung der adiabatischen Expanion des Treibgases
#------------------------
A_d = D_d^2.0 * pi/4.0;    # Fläche der Düse [m²]
A_g = D_g^2.0 * pi/4.0;    # Fläche des Raketengehäuses [m²]
vol_fl = fuellgrad * vol_rak;         # Füllvolumen der Flüssigkeit
vol_luft = (1-fuellgrad)*vol_rak;     # Gasvolumen in der Rakete
Energie = vol_rak * (1-fuellgrad)*Druck;  # Energie des Raketentreibsatzes der Rakete [J]
masse = m_geh + vol_fl*dichte_w + vol_luft*dichte_luft*Druck/100000.0;  # Startmasse der Rakete [kg], incl. Gasfüllung bei Anfangsdruck
Text=sprintf("Kaltwasserrakete mit p=%.2f bar, Füllgrad %.2f", Druck/100000, fuellgrad);
#
#
# Wir beginnen die Differentialgleichungen in kleinen Schritten bei Flughöhe=0 Geschwindigkeit=0 und Zeit=0
t = 0.0;
h = 0.0;
v = 0.0;
i = 1;
while h>=0.0 
    v_aus = sqrt(2.0*Druck/dichte_w);      #  Die Strahlaustrittsgeschwindigkeit, einfache Näherung ohne Berücksichtigung des Strömungsprofils in der Düse
    F_luft = A_g * v^2.0 * dichte_luft/2*cw;   #  Die Luftwiederstandskraft der Rakete
    dStrahlVol = A_d * v_aus*dt * 0.96;              #  Volumeneinheit die in der Zeit dt ausgestossen wird, und damit das Füllvolumen mindert, in simpler Korrekturfaktor soll die Reibung in der Düse berücksichtigen.
    if (vol_fl > 0.0) 
	dmasse = dStrahlVol*dichte_w;        #  Masseneinheit die in der Zeit dt ausgestossen wird
	F_schub = v_aus*dmasse/dt - g*masse - F_luft;   # Die effektive Schubkraft ist die Rückstosskraft, vermindert um die Schwerkraft und die Luftwiederstandskraft
	beschl = F_schub/masse;                   # Beschleunigung der Rakete
	masse = masse - dmasse;                   # Die Masse verringert sich um die Ausgestsenen Wassermasse
	Druck = Druck * (1.0/(1.0+dStrahlVol/vol_luft))^k_luft;      # Der Druck ändert sich mit der Volumenzunahme des Treibgases adiabatisch
	vol_luft = vol_luft - dStrahlVol;         # Das Druckluftvolumen nimmt um die ausgestossenen Wassermenge zu
	vol_fl = vol_fl - dStrahlVol;
    elseif
	F_schub =  -g * masse -F_luft;   # Die effektive Schubkraft ist die Rückstosskraft, vermindert um die Schwerkraft und die Luftwiederstandskraft
	beschl = F_schub/masse;                   # Beschleunigung der Rakete
    endif
    vv(i)=v;
    hh(i)=h;
    tt(i)=t;
    v = v + beschl*dt;
    h = h + v*dt;              # Die zurückgelegte Flugstrecke oder bei senkrechtem Flug die Steighöhe
    t = t + dt;
    i = i + 1;
#    printf("%f %f %f %f\n", t, h, v, Druck);
end
plot(tt, hh, ";h;",tt, vv,";v;")
title(Text)
xlabel("Zeit [s]")
ylabel("Hoehe [m], V [m/s]")
grid on
grid "minor" on
lim= ylim;
ylim ([0, lim(2)])
grid on
pause