diff --git a/Super-Pressure.py b/Super-Pressure.py
new file mode 100644
index 0000000..198663e
--- /dev/null
+++ b/Super-Pressure.py
@@ -0,0 +1,123 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+# VARIABLES:
+
+v0    = 0.00     # Initial Balloon Velocity in [m/s]
+s0    = 0.00     # Initial Balloon Altitude in [m]
+p0    = 101325   # (Initial) Air Pressure in [Pa]
+rho_a = 1.2250   # (Initial) Air Density in [kg/m^3]
+
+tmax = 50000
+
+T_a   = 288.15   # (Initial) Air Temperature in [K]
+g     = 9.80665  # local gravitation in [m/s^2]
+dt    = 0.01     # time step in [s]
+cD    = 0.47     # drag coefficient balloon (spherical) [-]
+
+R_a = 287.1      # Specific Gas Constant Dry Air in [J/(kg * K)]
+R_He = 2077.1    # Specific Gas Constant Helium in [J/(kg * K)]
+
+
+
+m_B = 50         # Balloon (Start) Mass in [kg]
+m_PL = 10        # Balloon Payload (Start) Mass in [kg]
+
+
+
+
+pi = np.pi
+
+def T(h):
+    if h >= 0 and h <= 11000:
+        res = 288.15 - 0.0065 * h
+        return res
+    elif h > 11000 and h <= 20000:
+        res = 216.65
+        return res
+    elif h >= 20000:
+        res = 216.65 + 0.0010 * (h - 20000)
+        return res
+
+def p(h):
+    if h >= 0 and h <= 11000:
+        res = 101325 * ((288.15 - 0.0065 * h)/288.15) ** 5.25577
+        return res
+    elif h > 11000 and h <= 20000:
+        res = 22632 * np.exp(-(h - 11000)/6341.62)
+        return res
+    elif h > 20000:
+        res = 5474.87 * ((216.65 + 0.0010 * (z - 20000))/216.65) ** (-34.163)
+        return res
+
+def rho(h):
+    res = p(h)/(R_a * T(h))
+    return res
+
+
+
+
+# p_a = p0         # air pressure = constant
+p_He = 101325       # initial gas pressure = air pressure
+T_g = 288.15        # gas temperature = air temperature = constant
+
+
+# GAS DENSITY
+rho_g = p_He / (R_He * T_g)
+
+V_B = m_B / rho_g
+D_B = (6 * V_B / pi) ** (1/3)
+
+
+v = []
+time = []
+alt = []
+
+
+v_z = v0
+z = s0
+t = 0
+D = 0
+
+rho_plot = []
+T_plot = []
+p_plot = []
+D_plot = []
+
+while t < tmax:
+    v.append(v_z)
+    alt.append(z)
+    time.append(t)
+
+    rho_plot.append(rho(z))
+    T_plot.append(T(z))
+    p_plot.append(p(z))
+    D_plot.append(D)
+
+    D = 0.5 * cD * 0.25 * pi * rho(z) * v_z ** 2 * D_B ** 2
+    I = g * V_B * (rho(z) - rho_g)
+    G = g * m_PL
+
+
+    dv_z = (I - np.sign(v_z) * D - G) / m_B * dt
+
+    v_z += dv_z * dt
+    z += v_z * dt
+
+    t += dt
+
+
+
+#plt.plot(time, rho_plot)
+#plt.plot(time, T_plot)
+###plt.plot(time, v)
+###plt.show()
+#plt.plot(time, ind)
+#plt.plot(time, D_plot)
+
+plt.plot(time, alt)
+plt.show()
+
+
+
+
diff --git a/input/__pycache__/natural_constants.cpython-38.pyc b/input/__pycache__/natural_constants.cpython-38.pyc
new file mode 100644
index 0000000..deed592
Binary files /dev/null and b/input/__pycache__/natural_constants.cpython-38.pyc differ
diff --git a/input/__pycache__/user_input.cpython-38.pyc b/input/__pycache__/user_input.cpython-38.pyc
new file mode 100644
index 0000000..a6b29c8
Binary files /dev/null and b/input/__pycache__/user_input.cpython-38.pyc differ
diff --git a/input/atmosphere.nc b/input/atmosphere.nc
index 06d7405..e69de29 100644
Binary files a/input/atmosphere.nc and b/input/atmosphere.nc differ
diff --git a/models/__init__.py b/models/__init__.py
index 06d7405..e69de29 100644
Binary files a/models/__init__.py and b/models/__init__.py differ
diff --git a/models/__pycache__/__init__.cpython-38.pyc b/models/__pycache__/__init__.cpython-38.pyc
new file mode 100644
index 0000000..cecabe0
Binary files /dev/null and b/models/__pycache__/__init__.cpython-38.pyc differ
diff --git a/models/__pycache__/drag.cpython-38.pyc b/models/__pycache__/drag.cpython-38.pyc
new file mode 100644
index 0000000..342001e
Binary files /dev/null and b/models/__pycache__/drag.cpython-38.pyc differ
diff --git a/models/__pycache__/simple_atmosphere.cpython-38.pyc b/models/__pycache__/simple_atmosphere.cpython-38.pyc
new file mode 100644
index 0000000..747855b
Binary files /dev/null and b/models/__pycache__/simple_atmosphere.cpython-38.pyc differ
diff --git a/models/__pycache__/sun.cpython-38.pyc b/models/__pycache__/sun.cpython-38.pyc
new file mode 100644
index 0000000..802efe8
Binary files /dev/null and b/models/__pycache__/sun.cpython-38.pyc differ
diff --git a/models/__pycache__/test1.cpython-38.pyc b/models/__pycache__/test1.cpython-38.pyc
new file mode 100644
index 0000000..a2fc0c6
Binary files /dev/null and b/models/__pycache__/test1.cpython-38.pyc differ
diff --git a/models/__pycache__/thermal.cpython-38.pyc b/models/__pycache__/thermal.cpython-38.pyc
new file mode 100644
index 0000000..7609bd3
Binary files /dev/null and b/models/__pycache__/thermal.cpython-38.pyc differ
diff --git a/models/test2.py b/models/test2.py
new file mode 100644
index 0000000..eda5730
--- /dev/null
+++ b/models/test2.py
@@ -0,0 +1,3 @@
+from test1 import f
+
+print(f(1))
diff --git a/models/thermal.py b/models/thermal.py
new file mode 100644
index 0000000..b96ec77
--- /dev/null
+++ b/models/thermal.py
@@ -0,0 +1,111 @@
+from input.natural_constants import *
+import numpy as np
+from datetime import datetime
+import time
+from astropy.time import Time
+
+from astropy.time import TimeISO
+class TimeYearDayTimeCustom(TimeISO):
+   """
+   day-of-year as "<DOY>".
+   The day-of-year (DOY) goes from 001 to 365 (366 in leap years).
+   The allowed subformat is:
+   - 'doy': day of year
+   """
+   name = 'doy'  # Unique format name
+   subfmts = (('doy',
+               '%j',
+               '{yday:03d}'),
+              ('doy',
+               '%j',
+               '{yday:03d}'),
+              ('doy',
+               '%j',
+               '{yday:03d}'))
+
+def AirMass(x, p_0, ELV, h):
+    p_air = x
+    ELV_rad = np.deg2rad(ELV)  # convert ELV from degree to radian
+    Dip = np.arccos(R_E / (R_E + h))  # Dip in radian
+
+    if ELV_rad >= -Dip and ELV_rad < 0:
+        res = p_air/p_0 * (1 + ELV_rad/Dip) - 70 * ELV_rad/Dip
+    else:
+        res = (p_air/p_0) * ((1229 + (614 * np.sin(ELV_rad)) ** 2) ** (1/2) - 614 * np.sin(ELV_rad))
+
+    return res
+
+#if ELV >= -(180/np.pi * np.arccos(R_e / (R_e + h))):
+#    tau_atm = 0.5 * (np.exp(-0.65 * AirMass(p_air, p_0, ELV, h)) + np.exp(-0.095 * AirMass(p_air, p_0, ELV, h)))
+#    tau_atmIR = 1.716 - 0.5 * (np.exp(-0.65 * p_air/p_0) + np.exp(-0.095 * p_air/p_0))
+#else:
+#    tau_atm = 0
+#    tau_atmIR = 0
+#
+#doy = int(utc.doy)
+#
+#MA = (357.52911 + 0.98560028 * (utc.jd - 2451545)) % 360  # in degree, reference: see folder "literature"
+#TA = MA + 2 * e * np.sin(np.deg2rad(MA)) + 5/4 * e ** 2 * np.sin(np.deg2rad(2 * MA))
+#I_Sun = 1367.5 * ((1 + e * np.cos(np.deg2rad(TA)))/(1 - e ** 2)) ** 2
+#I_SunZ = I_Sun * tau_atm
+#q_sun = I_SunZ
+#
+#q_IRground = epsilon_ground * sigma * T_ground ** 4
+#
+#q_IREarth = q_IRground * tau_atmIR
+#
+#if ELV <= 0:
+#    q_Albedo = 0
+#else:
+#    q_Albedo = Albedo * I_Sun * np.sin(np.deg2rad(ELV))
+#
+### NEEDED:
+# Diameter = ...
+# V_relative = ...
+# rho_air = ...
+# rho_gas = ...
+# T_film = ...
+# T_air = ...
+# T_gas = ...
+#
+#my_air = (1.458 * 10 ** -6 * T_air ** 1.5) / (T_air + 110.4)
+#my_gas = 1.895 * 10 ** -5 * (T_gas/273.15) ** 0.647
+#k_air = 0.0241 * (T_air/273.15) ** 0.9
+#k_gas = 0.144 * (T_gas/273.15) ** 0.7
+#Pr_air = 0.804 - 3.25 * 10 ** (-4) * T_air
+#Pr_gas = 0.729 - 1.6 * 10 ** (-4) * T_gas
+#
+#Gr_air = (rho_air ** 2 * g * np.abs(T_film - T_air) * Diameter ** 3) / (T_air * my_air ** 2)
+#Nu_air = 2 + 0.45 * (Gr_air * Pr_air) ** 0.25
+#HC_free = Nu_air * k_air / Diameter
+#Re = V_relative * Diameter * rho_air / my_air
+#HC_forced = k_air / Diameter * (2 + 0.41 * Re ** 0.55)
+#HC_internal = 0.13 * k_gas * ((rho_gas ** 2 * g * np.abs(T_film - T_gas) * Pr_gas) / (T_gas * my_air ** 2)) ** (1/3)
+#HC_external = np.maximum(HC_free, HC_forced)
+#
+#Q_Sun = alpha_VIS * A_proj * q_sun * (1 + tau_VIS/(1 - r_VIS))
+#Q_Albedo = alpha_VIS * A_surf * q_Albedo * ViewFactor * (1 + tau_VIS/(1 - r_VIS))
+#Q_IREarth = alpha_IR * A_surf * q_IREarth * ViewFactor * (1 + tau_IR/(1 - r_IR))
+#Q_IRFilm = sigma * epsilon * alpha_IR * A_surf * T_film ** 4 * 1/(1 - r_IR)
+#Q_IROut = sigma * epsilon * A_surf * T_film ** 4 * (1 * tau_IR/(1 - r_IR))
+#Q_ConvExt = HC_external * A_effective * (T_air - T_film)
+#Q_ConvInt = HC_internal * A_effective * (T_film - T_gas)
+#
+#HalfConeAngle = np.arcsin(R_E/(R_E + h))
+#ViewFactor = (1 - np.cos(HalfConeAngle))/2
+#
+#
+#
+#RoC = -V_z
+#
+#dT_gas = (Q_ConvInt/(gamma * c_v * M_gas) - (gamma - 1)/gamma * rho_air * g / (rho_gas * R_gas) * RoC) * dt
+#dT_film = ((Q_Sun + Q_Albedo + Q_IREarth + Q_IRfilm + Q_ConvExt - Q_ConvInt - Q_IRout)/(c_f * M_film)) * dt
+
+
+
+
+
+
+
+
+
diff --git a/output/placeholder.txt b/output/placeholder.txt
index 06d7405..e69de29 100644
Binary files a/output/placeholder.txt and b/output/placeholder.txt differ