diff --git a/models/__init__.py b/models/__init__.py
new file mode 100644
index 0000000..06d7405
Binary files /dev/null and b/models/__init__.py differ
diff --git a/models/drag.py b/models/drag.py
new file mode 100644
index 0000000..cac6e10
--- /dev/null
+++ b/models/drag.py
@@ -0,0 +1,6 @@
+import numpy as np
+
+c_d = 0.47       # drag coefficient balloon (spherical) [-]
+
+def drag(c_d, rho_air, d_b, v_z):
+    return 0.125 * np.pi * c_d * rho_air * (d_b * v_z) ** 2
\ No newline at end of file
diff --git a/models/simple_atmosphere.py b/models/simple_atmosphere.py
new file mode 100644
index 0000000..39cf97f
--- /dev/null
+++ b/models/simple_atmosphere.py
@@ -0,0 +1,28 @@
+from input.natural_constants import *
+import numpy as np
+
+# ATMOSPHERE MODEL:
+
+def T_air(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_air(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 * (h - 20000))/216.65) ** (-34.163)
+        return res
+def rho_air(h):
+    res = p_air(h)/(R_air * T_air(h))
+    return res
\ No newline at end of file
diff --git a/models/sun.py b/models/sun.py
new file mode 100644
index 0000000..85f3418
--- /dev/null
+++ b/models/sun.py
@@ -0,0 +1,58 @@
+import astropy.units as u
+import numpy as np
+from astropy.coordinates import EarthLocation, AltAz
+from astropy.coordinates import get_sun
+
+def sun_angles_astropy(lat, lon, h, utc):
+    loc = EarthLocation(lat=lat*u.deg, lon=lon*u.deg, height=h*u.m)
+    ref = AltAz(obstime=utc, location=loc)
+
+    sun_pos = get_sun(utc).transform_to(ref)
+
+    AZ = sun_pos.az.degree
+    ELV = sun_pos.alt.degree
+
+    return AZ, ELV
+
+def sun_angles_analytical(lat, lon, utc):
+    JD = utc.jd
+    JC = (JD - 2451545) / 36525
+    GML = (280.46646 + JC * (36000.76983 + JC * 0.0003032)) % 360
+    GMA = 357.52911 + JC * (35999.05029 - 0.0001537 * JC)
+    EEO = 0.016708634 - JC * (0.000042037 + 0.0000001267 * JC)
+    SEC = np.sin(np.deg2rad(GMA)) * (1.914602 - JC * (0.004817 + 0.000014 * JC)) + np.sin(np.deg2rad(2 * GMA)) * (
+                0.019993 - 0.000101 * JC) + np.sin(np.deg2rad(3 * GMA)) * 0.000289
+    STL = GML + SEC
+    # STA = GMA + SEC
+    # SRV = (1.000001018 * (1 - EEO ** 2)) / (1 + EEO * np.cos(np.deg2rad(STA)))
+    SAL = STL - 0.00569 - 0.00478 * np.sin(np.deg2rad(125.04 - 1934.136 * JC))
+    MOE = 23 + (26 + (21.448 - JC * (46.815 + JC * (0.00059 - JC * 0.001813))) / 60) / 60
+    OC = MOE + 0.00256 * np.cos(np.deg2rad(125.04 - 1934.136 * JC))
+    # SRA = np.rad2deg(np.arctan2(np.cos(np.deg2rad(OC)) * np.sin(np.deg2rad(SAL)), np.cos(np.deg2rad(SAL))))  # radian
+    SD = np.rad2deg(np.arcsin(np.sin(np.deg2rad(OC)) * np.sin(np.deg2rad(SAL))))  # radian
+    var_y = np.tan(np.deg2rad(OC / 2)) ** 2
+    EOT = 4 * np.rad2deg(
+        var_y * np.sin(2 * np.deg2rad(GML)) - 2 * EEO * np.sin(np.deg2rad(GMA)) + 4 * EEO * var_y * np.sin(
+            np.deg2rad(GMA)) * np.cos(2 * np.deg2rad(GML)) - 0.5 * var_y ** 2 * np.sin(
+            4 * np.deg2rad(GML)) - 1.25 * EEO ** 2 * np.sin(2 * np.deg2rad(GMA)))
+    TST = (((JD - 0.5) % 1) * 1440 + EOT + 4 * lon) % 1440
+
+    if TST / 4 < 0:
+        HA = TST / 4 + 180
+    else:
+        HA = TST / 4 - 180
+
+    SZA = np.rad2deg(np.arccos(
+        np.sin(np.deg2rad(lat)) * np.sin(np.deg2rad(SD)) + np.cos(np.deg2rad(lat)) * np.cos(np.deg2rad(SD)) * np.cos(
+            np.deg2rad(HA))))
+    SEA = 90 - SZA
+
+    if HA > 0:
+        SAA = (np.rad2deg(np.arccos(((np.sin(np.deg2rad(lat)) * np.cos(np.deg2rad(SZA))) - np.sin(np.deg2rad(SD))) / (
+                    np.cos(np.deg2rad(lat)) * np.sin(np.deg2rad(SZA))))) + 180) % 360
+    else:
+        SAA = (540 - np.rad2deg(np.arccos(
+            ((np.sin(np.deg2rad(lat)) * np.cos(np.deg2rad(SZA))) - np.sin(np.deg2rad(SD))) / (
+                        np.cos(np.deg2rad(lat)) * np.sin(np.deg2rad(SZA)))))) % 360
+
+    return SAA, SEA
\ No newline at end of file
diff --git a/models/test1.py b/models/test1.py
new file mode 100644
index 0000000..75275a1
--- /dev/null
+++ b/models/test1.py
@@ -0,0 +1,5 @@
+import numpy as np
+
+def f(x):
+    res = np.sin(x)
+    return res
\ No newline at end of file