From d1c15090dc618851432a908b0c8746d6c17f4a0d Mon Sep 17 00:00:00 2001
From: Brantegger Georg <Georg.brantegger@kelag.at>
Date: Wed, 27 Jul 2022 11:40:58 +0200
Subject: [PATCH] code cleanup: consistenly use getter and setter methods
 commenting etc

---
 .../Ausgleichsbecken_class_file.py            |  33 ++-
 .../Ausgleichsbecken_test_steady_state.ipynb  |  38 +--
 .../Druckrohrleitung_class_file.py            |  21 +-
 .../Druckrohrleitung_test_steady_state.ipynb  | 236 +++++++++------
 Kraftwerk/Kraftwerk_class_file.py             |  19 ++
 Kraftwerk/Kraftwerk_class_test.ipynb          | 105 +++++++
 ...gler_test.ipynb => Pegelregler_test.ipynb} | 275 ++++++++----------
 Regler/Regler_class_file.py                   | 117 +++++---
 Turbinen/Turbinen_class_file.py               | 121 ++++++--
 ...st.ipynb => Turbinen_test_Kennlinie.ipynb} |  78 +++--
 Turbinen/messy.ipynb                          |  90 ++++++
 Untertweng.ipynb                              | 170 +++++------
 Untertweng_mit_Pegelregler.ipynb              | 237 ++++++++-------
 13 files changed, 956 insertions(+), 584 deletions(-)
 create mode 100644 Kraftwerk/Kraftwerk_class_file.py
 create mode 100644 Kraftwerk/Kraftwerk_class_test.ipynb
 rename Regler/{regler_test.ipynb => Pegelregler_test.ipynb} (59%)
 rename Turbinen/{Turbinen_test.ipynb => Turbinen_test_Kennlinie.ipynb} (99%)
 create mode 100644 Turbinen/messy.ipynb

diff --git a/Ausgleichsbecken/Ausgleichsbecken_class_file.py b/Ausgleichsbecken/Ausgleichsbecken_class_file.py
index e8242d6..937e082 100644
--- a/Ausgleichsbecken/Ausgleichsbecken_class_file.py
+++ b/Ausgleichsbecken/Ausgleichsbecken_class_file.py
@@ -44,8 +44,8 @@ class Ausgleichsbecken_class:
     density_unit_print      = 'kg/m³'
     flux_unit_print         = 'm³/s'
     level_unit_print        = 'm'
-    time_unit_print         = 's'
     pressure_unit_print     = '--' # will be set by .set_pressure() method
+    time_unit_print         = 's'
     velocity_unit_print     = 'm/s'
     volume_unit_print       = 'm³'
 
@@ -102,8 +102,7 @@ class Ausgleichsbecken_class:
         # set the steady state (ss) condition in which the net flux is zero
             # set pressure acting on the outflux area so that the level stays constant
         ss_outflux = ss_influx
-        ss_outflux_vel = ss_outflux/self.area_outflux
-        ss_pressure = self.density*self.g*ss_level-ss_outflux_vel**2*self.density/2
+        ss_pressure = self.density*self.g*ss_level-(ss_outflux/self.area_outflux)**2*self.density/2
 
         self.set_influx(ss_influx)
         self.set_initial_level(ss_level)
@@ -169,6 +168,9 @@ class Ausgleichsbecken_class:
     def update_level(self,timestep):
         # update level based on net flux and timestep by calculating the volume change in
             # the timestep and the converting the new volume to a level by assuming a cuboid reservoir
+        
+        # cannot set new level directly in this method, because it gets called to calcuate during the Runge Kutta
+            # to calculate a ficticious level at half the timestep
         net_flux = self.influx-self.outflux
         delta_V = net_flux*timestep
         new_level = (self.volume+delta_V)/self.area
@@ -178,31 +180,32 @@ class Ausgleichsbecken_class:
         # sets volume in reservoir based on self.level
         return self.level*self.area 
 
+    def update_pressure(self):
+        p_new = self.density*self.g*self.level-(self.outflux/self.area_outflux)**2*self.density/2
+        return p_new
 
     def timestep_reservoir_evolution(self):
         # update outflux and outflux velocity based on current pipeline pressure and waterlevel in reservoir
-        yn      = self.outflux/self.area_outflux # outflux velocity
-        h       = self.level
         dt      = self.timestep
-        p       = self.pressure
-        # assume constant pipeline pressure during timestep
-            # e_RK_4 timestep is way smalle than timestep of characteristic method, so this should be a valid approx.
-            # (furthermore I have no idea how to approximate p_hs otherwise :/ )
-        p_hs    = self.pressure
-        A_a     = self.area_outflux
-        A       = self.area
-        h_hs    = self.update_level(dt/2)
         rho     = self.density
         g       = self.g
+        A       = self.area
+        A_a     = self.area_outflux
+        yn      = self.outflux/A_a # outflux velocity
+        h       = self.level
+        h_hs    = self.update_level(dt/2)
+        p       = self.pressure
+        p_hs    = self.pressure + rho*g*(h_hs-h)
         # explicit 4 step Runge Kutta
         Y1      = yn
-        Y2      = yn + dt/2*FODE_function(Y1,h,A,A_a,self.pressure,rho,g)
+        Y2      = yn + dt/2*FODE_function(Y1,h,A,A_a,p,rho,g)
         Y3      = yn + dt/2*FODE_function(Y2,h_hs,A,A_a,p_hs,rho,g)
         Y4      = yn + dt*FODE_function(Y3,h_hs,A,A_a,p_hs,rho,g)
         ynp1    = yn + dt/6*(FODE_function(Y1,h,A,A_a,p,rho,g)+2*FODE_function(Y2,h_hs,A,A_a,p_hs,rho,g)+ \
             2*FODE_function(Y3,h_hs,A,A_a,p_hs,rho,g)+ FODE_function(Y4,h,A,A_a,p,rho,g))
 
-        self.outflux    = ynp1*self.area_outflux
+        self.outflux    = ynp1*A_a
         self.level      = self.update_level(dt)
         self.volume     = self.update_volume()
+        self.pressure   = self.update_pressure()
         
diff --git a/Ausgleichsbecken/Ausgleichsbecken_test_steady_state.ipynb b/Ausgleichsbecken/Ausgleichsbecken_test_steady_state.ipynb
index 8838088..1b82115 100644
--- a/Ausgleichsbecken/Ausgleichsbecken_test_steady_state.ipynb
+++ b/Ausgleichsbecken/Ausgleichsbecken_test_steady_state.ipynb
@@ -28,9 +28,9 @@
     "# define constants\n",
     "initial_level               = 10.        # m\n",
     "initial_influx              = 5.       # m³/s\n",
-    "initial_outflux             = 1.        # m³/s\n",
-    "initial_pipeline_pressure   = 10.\n",
-    "initial_pressure_unit       = 'mWS'\n",
+    "# initial_outflux             = 1.        # m³/s\n",
+    "# initial_pipeline_pressure   = 10.\n",
+    "# initial_pressure_unit       = 'mWS'\n",
     "conversion_pressure_unit    = 'mWS'\n",
     "\n",
     "area_base                   = 1.        # m²\n",
@@ -41,12 +41,12 @@
     "\n",
     "# for while loop\n",
     "total_min_level             = 0.01      # m\n",
-    "total_max_time              = 1000       # s"
+    "total_max_time              = 1000     # s"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 3,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -72,7 +72,6 @@
     "i_max = -1\n",
     "\n",
     "for i in range(np.size(time_vec)-1):\n",
-    "    # update to include p_halfstep\n",
     "    V.set_pressure(pressure_vec[i])\n",
     "    V.timestep_reservoir_evolution()\n",
     "    outflux_vec[i+1] = V.get_current_outflux()\n",
@@ -85,12 +84,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": 7,
    "metadata": {},
    "outputs": [],
    "source": [
     "\n",
-    "fig1, (ax1, ax2, ax3, ax4) = plt.subplots(4, 1)\n",
+    "fig1, (ax1, ax2, ax3) = plt.subplots(3, 1)\n",
     "fig1.set_figheight(10)\n",
     "fig1.suptitle('Ausgleichsbecken')\n",
     "\n",
@@ -109,29 +108,6 @@
     "ax3.set_xlabel(r'$t$ ['+V.time_unit+']')\n",
     "ax3.legend()\n",
     "\n",
-    "# plt.subplots_adjust(left=0.2, bottom=0.2)\n",
-    "ax4.set_axis_off()\n",
-    "cell_text = np.array([[level_vec[0], V.level_unit], \\\n",
-    "            [initial_influx, V.flux_unit], \\\n",
-    "            [outflux_vec[0], V.flux_unit], \\\n",
-    "            [simulation_timestep, V.time_unit], \\\n",
-    "            [area_base, V.area_unit], \\\n",
-    "            [area_outflux, V.area_unit]])\n",
-    "\n",
-    "row_labels =['initial_level', \\\n",
-    "            'initial_influx', \\\n",
-    "            'initial_outflux', \\\n",
-    "            'simulation_timestep', \\\n",
-    "            'area_base', \\\n",
-    "            'area_outflux']\n",
-    "\n",
-    "plt.table(cellText=cell_text, \\\n",
-    "            cellLoc='center', \\\n",
-    "            colWidths=[0.3,0.1,0.3], \\\n",
-    "            rowLabels=row_labels, \\\n",
-    "            loc = 1, \\\n",
-    "            rowLoc='left', \\\n",
-    "            fontsize = 15.)\n",
     "\n",
     "fig1.tight_layout()            "
    ]
diff --git a/Druckrohrleitung/Druckrohrleitung_class_file.py b/Druckrohrleitung/Druckrohrleitung_class_file.py
index df9a3c1..71d553b 100644
--- a/Druckrohrleitung/Druckrohrleitung_class_file.py
+++ b/Druckrohrleitung/Druckrohrleitung_class_file.py
@@ -33,6 +33,8 @@ class Druckrohrleitung_class:
         self.density    = rho                                           # density of the liquid in the pipeline
         self.g          = g                                             # gravitational acceleration
         
+        self.A          = (diameter/2)**2*np.pi
+
         self.dx         = total_length/number_segments                  # length of each segment
         self.l_vec      = np.arange(0,(number_segments+1),1)*self.dx    # vector giving the distance from each node to the start of the pipeline
 
@@ -98,23 +100,22 @@ class Druckrohrleitung_class:
         p_old_res               = self.p_old[1]     # @ second  node (the one after the reservoir)
         v_old_res               = self.v_old[1]     # @ second  node (the one after the reservoir)
         # set the boundary conditions derived from reservoir and turbine
-        self.v_boundary_tur     = v_turbine         # at new timestep
-        self.p_boundary_res     = p_reservoir       # at new timestep
+        v_boundary_tur      = v_turbine         # at new timestep
+        p_boundary_res      = p_reservoir       # at new timestep
         # calculate the missing boundary conditions
-        self.v_boundary_res     = v_old_res+1/(rho*c)*(p_reservoir-p_old_res)+dt*g*np.sin(alpha)-f_D*dt/(2*D)*abs(v_old_res)*v_old_res
-        self.p_boundary_tur     = p_old_tur-rho*c*(v_turbine-v_old_tur)+rho*c*dt*g*np.sin(alpha)-f_D*rho*c*dt/(2*D)*abs(v_old_tur)*v_old_tur
+        v_boundary_res      = v_old_res+1/(rho*c)*(p_boundary_res-p_old_res)+dt*g*np.sin(alpha)-f_D*dt/(2*D)*abs(v_old_res)*v_old_res
+        p_boundary_tur      = p_old_tur-rho*c*(v_boundary_tur-v_old_tur)+rho*c*dt*g*np.sin(alpha)-f_D*rho*c*dt/(2*D)*abs(v_old_tur)*v_old_tur
 
         # write boundary conditions to the velocity/pressure vectors of the next timestep
-        self.v[0]               = self.v_boundary_res.copy()
-        self.v[-1]              = self.v_boundary_tur.copy()
-        self.p[0]               = self.p_boundary_res.copy()
-        self.p[-1]              = self.p_boundary_tur.copy()
-
+        self.v[0]           = v_boundary_res
+        self.v[-1]          = v_boundary_tur
+        self.p[0]           = p_boundary_res
+        self.p[-1]          = p_boundary_tur
 
     def set_steady_state(self,ss_flux,ss_level_reservoir,pl_vec,h_vec):
         # set the pressure and velocity distributions, that allow a constant flow of water from the (steady-state) reservoir to the (steady-state) turbine
             # the flow velocity is given by the constant flow through the pipe
-        ss_v0 = np.full(self.n_seg+1,ss_flux/(self.dia**2/4*np.pi))
+        ss_v0 = np.full(self.n_seg+1,ss_flux/self.A)
             # the static pressure is given by the hydrostatic pressure, corrected for friction losses and dynamic pressure
         ss_pressure = (self.density*self.g*(ss_level_reservoir+h_vec)-ss_v0**2*self.density/2)-(self.f_D*pl_vec/self.dia*self.density/2*ss_v0**2)
 
diff --git a/Druckrohrleitung/Druckrohrleitung_test_steady_state.ipynb b/Druckrohrleitung/Druckrohrleitung_test_steady_state.ipynb
index bc464e4..e92c8b6 100644
--- a/Druckrohrleitung/Druckrohrleitung_test_steady_state.ipynb
+++ b/Druckrohrleitung/Druckrohrleitung_test_steady_state.ipynb
@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -16,64 +16,98 @@
     "current = os.path.dirname(os.path.realpath('Main_Programm.ipynb'))\n",
     "parent = os.path.dirname(current)\n",
     "sys.path.append(parent)\n",
-    "from functions.pressure_conversion import pressure_conversion"
+    "from functions.pressure_conversion import pressure_conversion\n",
+    "from Ausgleichsbecken.Ausgleichsbecken_class_file import Ausgleichsbecken_class"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
     "%matplotlib qt5\n",
-    "#define constants\n",
+    "#define constants pipe\n",
     "\n",
-    "g       = 9.81              # gravitational acceleration [m/s²]\n",
-    "rho     = 1000              # density of water [kg/m³]\n",
+    "g       = 9.81                  # gravitational acceleration [m/s²]\n",
+    "rho     = 1000.                 # density of water [kg/m³]\n",
     "\n",
-    "L       = 1000              # length of pipeline [m]\n",
-    "D       = 1                 # pipe diameter [m]\n",
-    "Q0      = 2                 # initial flow in whole pipe [m³/s]\n",
-    "h_res   = 20                # water level in upstream reservoir [m]\n",
-    "n       = 10                # number of pipe segments in discretization\n",
-    "nt      = 100               # number of time steps after initial conditions\n",
-    "f_D     = 0.01              # Darcy friction factor\n",
-    "c       = 400               # propagation velocity of the pressure wave [m/s]\n",
-    "h_pipe  = 200              # hydraulic head without reservoir [m] \n",
-    "alpha   = np.arcsin(h_pipe/L)  # Höhenwinkel der Druckrohrleitung \n",
+    "L       = 1000.                 # length of pipeline [m]\n",
+    "D       = 0.9                   # pipe diameter [m]\n",
+    "h_res   = 10.                   # water level in upstream reservoir [m]\n",
+    "n       = 50                    # number of pipe segments in discretization\n",
+    "nt      = 5000                   # number of time steps after initial conditions\n",
+    "f_D     = 0.01                  # Darcy friction factor\n",
+    "c       = 400.                  # propagation velocity of the pressure wave [m/s]\n",
+    "h_pipe  = 105.                  # hydraulic head without reservoir [m] \n",
+    "alpha   = np.arcsin(h_pipe/L)   # Höhenwinkel der Druckrohrleitung \n",
     "\n",
     "\n",
     "# preparing the discretization and initial conditions\n",
-    "initial_influx  = 2.                                    # m³/s\n",
-    "initial_level   = 10.                                   # m\n",
+    "initial_flux    = 0.8                                    # m³/s\n",
+    "initial_level   = h_res                                   # m\n",
     "dx              = L/n                                   # length of each pipe segment\n",
     "dt              = dx/c                                  # timestep according to method of characterisitics\n",
     "nn              = n+1                                   # number of nodes\n",
     "pl_vec          = np.arange(0,nn*dx,dx)                 # pl = pipe-length. position of the nodes on the pipeline\n",
     "t_vec           = np.arange(0,nt*dt,dt)                 # time vector\n",
-    "h_vec           = np.arange(0,h_pipe+h_pipe/n,h_pipe/n) # hydraulic head of pipeline at each node\n"
+    "h_vec           = np.arange(0,h_pipe+h_pipe/n,h_pipe/n) # hydraulic head of pipeline at each node\n",
+    "\n",
+    "\n",
+    "# define constants reservoir\n",
+    "conversion_pressure_unit    = 'mWS'\n",
+    "\n",
+    "area_base                   = 75.                # m²\n",
+    "area_pipe                   = (D/2)**2*np.pi    # m²\n",
+    "critical_level_low          = 0.                # m\n",
+    "critical_level_high         = 100.              # m\n",
+    "\n",
+    "# make sure e-RK4 method of reservoir has a small enough timestep to avoid runaway numerical error\n",
+    "nt_eRK4                     = 1              # number of simulation steps of reservoir in between timesteps of pipeline              \n",
+    "simulation_timestep         = dt/nt_eRK4"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
+    "V = Ausgleichsbecken_class(area_base, area_pipe, critical_level_low, critical_level_high,simulation_timestep)\n",
+    "V.set_steady_state(initial_flux,initial_level,conversion_pressure_unit)\n",
+    "\n",
     "pipe = Druckrohrleitung_class(L,D,n,alpha,f_D)\n",
     "pipe.set_pressure_propagation_velocity(c)\n",
     "pipe.set_number_of_timesteps(nt)\n",
-    "pipe.set_steady_state(initial_influx,initial_level,pl_vec,h_vec)"
+    "pipe.set_steady_state(initial_flux,initial_level,pl_vec,h_vec)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(V.get_current_influx())\n",
+    "print(V.get_current_outflux())\n",
+    "print(V.get_current_level())\n",
+    "print(V.get_current_pressure())\n",
+    "print(pipe.get_current_pressure_distribution()[0])\n",
+    "print(pipe.get_current_velocity_distribution()*area_pipe)\n",
+    "print(pipe.get_current_velocity_distribution())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
    "metadata": {},
    "outputs": [],
    "source": [
     "# initialization for timeloop\n",
     "\n",
+    "level_vec       = np.zeros_like(t_vec)\n",
+    "level_vec[0]    = V.get_current_level()\n",
+    "\n",
     "# prepare the vectors in which the pressure and velocity distribution in the pipeline from the previous timestep are stored\n",
     "v_old = pipe.get_current_velocity_distribution()\n",
     "p_old = pipe.get_current_pressure_distribution()\n",
@@ -94,80 +128,110 @@
     "p_boundary_tur[0]   = p_old[-1]\n"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": 5,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "fig2,axs2 = plt.subplots(2,1)\n",
-    "axs2[0].set_title('Pressure distribution in pipeline')\n",
-    "axs2[0].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
-    "axs2[0].set_ylabel(r'$p$ [mWS]')\n",
-    "lo_00, = axs2[0].plot(pl_vec,pressure_conversion(p_old,'Pa','mWS'),marker='.')\n",
-    "axs2[0].set_ylim([0.9*np.min(pressure_conversion(p_old,'Pa','mWS')),1.1*np.max(pressure_conversion(p_old,'Pa','mWS'))])\n",
-    "\n",
-    "axs2[1].set_title('Velocity distribution in pipeline')\n",
-    "axs2[1].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
-    "axs2[1].set_ylabel(r'$p$ [mWS]')\n",
-    "lo_01, = axs2[1].plot(pl_vec,v_old,marker='.')\n",
-    "axs2[1].set_ylim([0.9*np.min(v_old),1.1*np.max(v_boundary_res)])\n",
-    "\n",
-    "fig2.tight_layout()\n",
-    "plt.pause(5)\n",
-    "\n",
-    "\n",
-    "for it in range(1,pipe.nt):\n",
-    "    pipe.set_boundary_conditions_next_timestep(p_boundary_res[0],v_boundary_tur[0])\n",
-    "    pipe.timestep_characteristic_method()\n",
-    "    lo_00.set_ydata(pressure_conversion(pipe.get_current_pressure_distribution(),'Pa','mWS'))\n",
-    "    lo_01.set_ydata(pipe.get_current_velocity_distribution())\n",
-    "\n",
-    "    v_boundary_res[it] = pipe.get_current_velocity_distribution()[0]\n",
-    "    v_boundary_tur[it] = pipe.get_current_velocity_distribution()[-1]\n",
-    "    p_boundary_res[it] = pipe.get_current_pressure_distribution()[0]\n",
-    "    p_boundary_tur[it] = pipe.get_current_pressure_distribution()[-1]\n",
-    "\n",
-    "\n",
-    "    \n",
-    "    fig2.suptitle(str(it))\n",
-    "    fig2.canvas.draw()\n",
-    "    fig2.tight_layout()\n",
-    "    plt.pause(0.2)\n"
-   ]
-  },
   {
    "cell_type": "code",
    "execution_count": 6,
    "metadata": {},
    "outputs": [],
    "source": [
-    "fig3,axs3 = plt.subplots(2,2)\n",
-    "axs3[0,0].set_title('Pressure Reservoir')\n",
-    "axs3[0,0].plot(t_vec,pressure_conversion(p_boundary_res,'Pa','mWS'))\n",
-    "axs3[0,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
-    "axs3[0,0].set_ylabel(r'$p$ [mWS]')\n",
-    "axs3[0,0].set_ylim([0.9*np.min(pressure_conversion(p_boundary_res,'Pa','mWS')),1.1*np.max(pressure_conversion(p_boundary_res,'Pa','mWS'))])\n",
+    "fig1,axs1 = plt.subplots(2,1)\n",
+    "axs1[0].set_title('Pressure distribution in pipeline')\n",
+    "axs1[0].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
+    "axs1[0].set_ylabel(r'$p$ [mWS]')\n",
+    "lo_00, = axs1[0].plot(pl_vec,pressure_conversion(p_old,'Pa',conversion_pressure_unit),marker='.')\n",
+    "axs1[0].set_ylim([0.9*np.min(pressure_conversion(p_old,'Pa',conversion_pressure_unit)),1.1*np.max(pressure_conversion(p_old,'Pa',conversion_pressure_unit))])\n",
     "\n",
-    "axs3[0,1].set_title('Velocity Reservoir')\n",
-    "axs3[0,1].plot(t_vec,v_boundary_res)\n",
-    "axs3[0,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
-    "axs3[0,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n",
-    "axs3[0,1].set_ylim([0.9*np.min(v_boundary_res),1.1*np.max(v_boundary_res)])\n",
+    "axs1[1].set_title('Velocity distribution in pipeline')\n",
+    "axs1[1].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
+    "axs1[1].set_ylabel(r'$v$ [m/s]')\n",
+    "lo_01, = axs1[1].plot(pl_vec,v_old,marker='.')\n",
+    "# axs1[1].set_ylim([0.9*np.min(v_old),1.1*np.max(v_boundary_res)])\n",
     "\n",
-    "axs3[1,0].set_title('Pressure Turbine')\n",
-    "axs3[1,0].plot(t_vec,pressure_conversion(p_boundary_tur,'Pa','mWS'))\n",
-    "axs3[1,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
-    "axs3[1,0].set_ylabel(r'$p$ [mWS]')\n",
-    "axs3[1,0].set_ylim([0.9*np.min(pressure_conversion(p_boundary_tur,'Pa','mWS')),1.1*np.max(pressure_conversion(p_boundary_tur,'Pa','mWS'))])\n",
+    "fig1.tight_layout()\n",
+    "plt.pause(1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [],
+   "source": [
     "\n",
-    "axs3[1,1].set_title('Velocity Turbine')\n",
-    "axs3[1,1].plot(t_vec,v_boundary_tur)\n",
-    "axs3[1,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
-    "axs3[1,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n",
-    "axs3[1,1].set_ylim([0.9*np.min(v_boundary_tur),1.1*np.max(v_boundary_tur)])\n",
+    "for it_pipe in range(1,nt):\n",
+    "# for each pipeline timestep, execute nt_eRK4 timesteps of the reservoir code\n",
+    "    # set initial conditions for the reservoir time evolution calculted with e-RK4\n",
+    "    V.set_pressure      = p_old[0]\n",
+    "    V.set_outflux       = v_old[0]*area_pipe\n",
+    "    # calculate the time evolution of the reservoir level within each pipeline timestep to avoid runaway numerical error\n",
+    "    for it_res in range(nt_eRK4):\n",
+    "        V.timestep_reservoir_evolution()                                                             \n",
+    "    level_vec[it_pipe] = V.get_current_level() \n",
     "\n",
-    "fig3.tight_layout()\n",
+    "    \n",
+    "    # set boundary conditions for the next timestep of the characteristic method\n",
+    "    p_boundary_res[it_pipe] = V.get_current_pressure()\n",
+    "    v_boundary_tur[it_pipe] = initial_flux/area_pipe\n",
+    "\n",
+    "    # the the boundary conditions in the pipe.object and thereby calculate boundary pressure at turbine\n",
+    "    pipe.set_boundary_conditions_next_timestep(p_boundary_res[it_pipe],v_boundary_tur[it_pipe])\n",
+    "    p_boundary_tur[it_pipe] = pipe.get_current_pressure_distribution()[-1]\n",
+    "    v_boundary_res[it_pipe] = pipe.get_current_velocity_distribution()[0]\n",
+    "\n",
+    "    # perform the next timestep via the characteristic method\n",
+    "    pipe.timestep_characteristic_method()\n",
+    "\n",
+    "    # prepare for next loop\n",
+    "    p_old = pipe.get_current_pressure_distribution()\n",
+    "    v_old = pipe.get_current_velocity_distribution()\n",
+    "\n",
+    "    # plot some stuff\n",
+    "        # remove line-objects to autoscale axes (there is definetly a better way, but this works ¯\\_(ツ)_/¯ )\n",
+    "    lo_00.remove()\n",
+    "    lo_01.remove()\n",
+    "    # lo_02.remove()\n",
+    "        # plot new pressure and velocity distribution in the pipeline\n",
+    "    lo_00, = axs1[0].plot(pl_vec,pressure_conversion(p_old,'Pa', conversion_pressure_unit),marker='.',c='blue')\n",
+    "    lo_01, = axs1[1].plot(pl_vec,v_old,marker='.',c='blue')\n",
+    "    \n",
+    "    fig1.suptitle(str(round(t_vec[it_pipe],2)) + '/' + str(t_vec[-1]))\n",
+    "    fig1.canvas.draw()\n",
+    "    fig1.tight_layout()\n",
+    "    plt.pause(0.00001)\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig2,axs2 = plt.subplots(2,2)\n",
+    "axs2[0,0].set_title('Pressure Reservoir')\n",
+    "axs2[0,0].plot(t_vec,pressure_conversion(p_boundary_res,'Pa',conversion_pressure_unit))\n",
+    "axs2[0,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
+    "axs2[0,0].set_ylabel(r'$p$ [mWS]')\n",
+    "axs2[0,0].set_ylim([0.9*np.min(pressure_conversion(p_boundary_res,'Pa',conversion_pressure_unit)),1.1*np.max(pressure_conversion(p_boundary_res,'Pa',conversion_pressure_unit))])\n",
+    "\n",
+    "axs2[0,1].set_title('Velocity Reservoir')\n",
+    "axs2[0,1].plot(t_vec,v_boundary_res)\n",
+    "axs2[0,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
+    "axs2[0,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n",
+    "axs2[0,1].set_ylim([0.9*np.min(v_boundary_res),1.1*np.max(v_boundary_res)])\n",
+    "\n",
+    "axs2[1,0].set_title('Pressure Turbine')\n",
+    "axs2[1,0].plot(t_vec,pressure_conversion(p_boundary_tur,'Pa',conversion_pressure_unit))\n",
+    "axs2[1,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
+    "axs2[1,0].set_ylabel(r'$p$ [mWS]')\n",
+    "axs2[1,0].set_ylim([0.9*np.min(pressure_conversion(p_boundary_tur,'Pa',conversion_pressure_unit)),1.1*np.max(pressure_conversion(p_boundary_tur,'Pa',conversion_pressure_unit))])\n",
+    "\n",
+    "axs2[1,1].set_title('Velocity Turbine')\n",
+    "axs2[1,1].plot(t_vec,v_boundary_tur)\n",
+    "axs2[1,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
+    "axs2[1,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n",
+    "axs2[1,1].set_ylim([0.9*np.min(v_boundary_tur),1.1*np.max(v_boundary_tur)])\n",
+    "\n",
+    "fig2.tight_layout()\n",
     "plt.show()"
    ]
   }
diff --git a/Kraftwerk/Kraftwerk_class_file.py b/Kraftwerk/Kraftwerk_class_file.py
new file mode 100644
index 0000000..efcabd4
--- /dev/null
+++ b/Kraftwerk/Kraftwerk_class_file.py
@@ -0,0 +1,19 @@
+#importing Druckrohrleitung
+import sys
+import os
+current = os.path.dirname(os.path.realpath('Main_Programm.ipynb'))
+parent = os.path.dirname(current)
+sys.path.append(parent)
+from functions.pressure_conversion import pressure_conversion
+from Turbinen.Turbinen_class_file import Francis_Turbine
+
+class Kraftwerk_class:
+    def __init__(self):
+        self.turbines = []
+
+    def add_turbine(self,turbine):
+        self.turbines.append(turbine)
+
+    def print_info(self):
+        for turbine in self.turbines:
+            turbine.get_info(full=True)
diff --git a/Kraftwerk/Kraftwerk_class_test.ipynb b/Kraftwerk/Kraftwerk_class_test.ipynb
new file mode 100644
index 0000000..a820c4c
--- /dev/null
+++ b/Kraftwerk/Kraftwerk_class_test.ipynb
@@ -0,0 +1,105 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import sys\n",
+    "import os\n",
+    "from Kraftwerk_class_file import Kraftwerk_class\n",
+    "\n",
+    "current = os.path.dirname(os.path.realpath('Main_Programm.ipynb'))\n",
+    "parent = os.path.dirname(current)\n",
+    "sys.path.append(parent)\n",
+    "from functions.pressure_conversion import pressure_conversion\n",
+    "from Turbinen.Turbinen_class_file import Francis_Turbine"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[<Turbinen.Turbinen_class_file.Francis_Turbine object at 0x0000018A94FDDE80>, <Turbinen.Turbinen_class_file.Francis_Turbine object at 0x0000018A94FDDEE0>]\n",
+      "Turbine has the following attributes: \n",
+      "----------------------------- \n",
+      "Type                  =       Francis \n",
+      "Nominal flux          =       0.85       m³/s \n",
+      "Nominal pressure      =       108.09     mWS\n",
+      "Nominal LA            =       100.0      % \n",
+      "Closing time          =       500        s \n",
+      "Current flux          =       -1.0       m³/s \n",
+      "Current pipe pressure =       -1.0       mWS \n",
+      "Current LA            =       -1.0       % \n",
+      "Simulation timestep   =       -1.0       s \n",
+      "----------------------------- \n",
+      "\n",
+      "Turbine has the following attributes: \n",
+      "----------------------------- \n",
+      "Type                  =       Francis \n",
+      "Nominal flux          =       0.85       m³/s \n",
+      "Nominal pressure      =       108.09     mWS\n",
+      "Nominal LA            =       100.0      % \n",
+      "Closing time          =       500        s \n",
+      "Current flux          =       -1.0       m³/s \n",
+      "Current pipe pressure =       -1.0       mWS \n",
+      "Current LA            =       -1.0       % \n",
+      "Simulation timestep   =       -1.0       s \n",
+      "----------------------------- \n",
+      "\n"
+     ]
+    }
+   ],
+   "source": [
+    "#Turbine\n",
+    "Q_nenn          = 0.85  # m³/s\n",
+    "p_nenn          = pressure_conversion(10.6,'bar','Pa')\n",
+    "closing_time    = 500   #s\n",
+    "\n",
+    "T1 = Francis_Turbine(Q_nenn,p_nenn,closing_time)\n",
+    "T2 = Francis_Turbine(Q_nenn,p_nenn,closing_time)\n",
+    "\n",
+    "KW = Kraftwerk_class()\n",
+    "KW.add_turbine(T1)\n",
+    "KW.add_turbine(T2)\n",
+    "\n",
+    "print(KW.turbines)\n",
+    "\n",
+    "KW.print_info()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3.8.13 ('Georg_DT_Slot3')",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.13"
+  },
+  "orig_nbformat": 4,
+  "vscode": {
+   "interpreter": {
+    "hash": "84fb123bdc47ab647d3782661abcbe80fbb79236dd2f8adf4cef30e8755eb2cd"
+   }
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/Regler/regler_test.ipynb b/Regler/Pegelregler_test.ipynb
similarity index 59%
rename from Regler/regler_test.ipynb
rename to Regler/Pegelregler_test.ipynb
index bad0931..081cd94 100644
--- a/Regler/regler_test.ipynb
+++ b/Regler/Pegelregler_test.ipynb
@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 34,
+   "execution_count": 1,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -23,15 +23,16 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 35,
+   "execution_count": 2,
    "metadata": {},
    "outputs": [],
    "source": [
     "#define constants\n",
     "\n",
     "#Turbine\n",
-    "Q_nenn = 0.85\n",
-    "p_nenn,_ = pressure_conversion(10.6,'bar','Pa')\n",
+    "Q_nenn          = 0.85  # m³/s\n",
+    "p_nenn          = pressure_conversion(10.6,'bar','Pa')\n",
+    "closing_time    = 480.   #s\n",
     "\n",
     "# physics\n",
     "g               = 9.81                          # gravitational acceleration [m/s²]\n",
@@ -39,8 +40,8 @@
     "\n",
     "# define controller constants\n",
     "target_level = 8.  # m\n",
-    "Kp = 0.1\n",
-    "Ti = 100.\n",
+    "Kp = 0.01\n",
+    "Ti = 3600.\n",
     "deadband_range = 0.05 # m\n",
     "\n",
     "# reservoir\n",
@@ -59,9 +60,9 @@
     "h_fict = 100\n",
     "offset_pressure = rho*g*h_fict\n",
     "\n",
-    "t_max                       = 1e3 #s\n",
-    "nt                          = int(1e6)      # number of simulation steps of reservoir in between timesteps of pipeline              \n",
-    "dt                          = t_max/nt\n",
+    "t_max                       = 1e4 #s\n",
+    "dt                          = 1e-2      # simulation timestep\n",
+    "nt                          = int(t_max//dt)    # number of simulation steps of reservoir in between timesteps of pipeline              \n",
     "\n",
     "t_vec = np.arange(0,nt+1,1)*dt\n",
     "\n"
@@ -69,25 +70,24 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 36,
+   "execution_count": 3,
    "metadata": {},
    "outputs": [],
    "source": [
     "# create objects\n",
     "\n",
     "V = Ausgleichsbecken_class(area_base,area_outflux,critical_level_low,critical_level_high,dt)\n",
-    "V.set_steady_state(initial_influx,initial_level,initial_pressure_unit,conversion_pressure_unit)\n",
+    "V.set_steady_state(initial_influx,initial_level,conversion_pressure_unit)\n",
     "\n",
-    "T1 = Francis_Turbine(Q_nenn,p_nenn)\n",
+    "T1 = Francis_Turbine(Q_nenn,p_nenn,closing_time,dt)\n",
     "T1.set_steady_state(initial_influx,p0+offset_pressure)\n",
-    "T1.set_closing_time(500)\n",
     "\n",
     "Pegelregler = PI_controller_class(target_level,deadband_range,Kp,Ti,dt)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 37,
+   "execution_count": 4,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -96,12 +96,12 @@
     "LA_soll_vec     = np.full(nt+1,T1.LA)\n",
     "Q_vec           = np.full(nt+1,initial_influx)\n",
     "\n",
-    "Pegelregler.control_variable = T1.LA"
+    "Pegelregler.control_variable = T1.get_current_LA()"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 38,
+   "execution_count": 5,
    "metadata": {},
    "outputs": [
     {
@@ -109,106 +109,105 @@
      "output_type": "stream",
      "text": [
       "0.0\n",
-      "10.0\n",
-      "20.0\n",
-      "30.0\n",
-      "40.0\n",
-      "50.0\n",
-      "60.0\n",
-      "70.0\n",
-      "80.0\n",
-      "90.0\n",
       "100.0\n",
-      "110.0\n",
-      "120.0\n",
-      "130.0\n",
-      "140.0\n",
-      "150.0\n",
-      "160.0\n",
-      "170.0\n",
-      "180.0\n",
-      "190.0\n",
       "200.0\n",
-      "210.0\n",
-      "220.0\n",
-      "230.0\n",
-      "240.0\n",
-      "250.0\n",
-      "260.0\n",
-      "270.0\n",
-      "280.0\n",
-      "290.0\n",
       "300.0\n",
-      "310.0\n",
-      "320.0\n",
-      "330.0\n",
-      "340.0\n",
-      "350.0\n",
-      "360.0\n",
-      "370.0\n",
-      "380.0\n",
-      "390.0\n",
       "400.0\n",
-      "410.0\n",
-      "420.0\n",
-      "430.0\n",
-      "440.0\n",
-      "450.0\n",
-      "460.0\n",
-      "470.0\n",
-      "480.0\n",
-      "490.0\n",
       "500.0\n",
-      "510.0\n",
-      "520.0\n",
-      "530.0\n",
-      "540.0\n",
-      "550.0\n",
-      "560.0\n",
-      "570.0\n",
-      "580.0\n",
-      "590.0\n",
       "600.0\n",
-      "610.0\n",
-      "620.0\n",
-      "630.0\n",
-      "640.0\n",
-      "650.0\n",
-      "660.0\n",
-      "670.0\n",
-      "680.0\n",
-      "690.0\n",
       "700.0\n",
-      "710.0\n",
-      "720.0\n",
-      "730.0\n",
-      "740.0\n",
-      "750.0\n",
-      "760.0\n",
-      "770.0\n",
-      "780.0\n",
-      "790.0\n",
       "800.0\n",
-      "810.0\n",
-      "820.0\n",
-      "830.0\n",
-      "840.0\n",
-      "850.0\n",
-      "860.0\n",
-      "870.0\n",
-      "880.0\n",
-      "890.0\n",
       "900.0\n",
-      "910.0\n",
-      "920.0\n",
-      "930.0\n",
-      "940.0\n",
-      "950.0\n",
-      "960.0\n",
-      "970.0\n",
-      "980.0\n",
-      "990.0\n",
-      "1000.0\n"
+      "1000.0\n",
+      "1100.0\n",
+      "1200.0\n",
+      "1300.0\n",
+      "1400.0\n",
+      "1500.0\n",
+      "1600.0\n",
+      "1700.0\n",
+      "1800.0\n",
+      "1900.0\n",
+      "2000.0\n",
+      "2100.0\n",
+      "2200.0\n",
+      "2300.0\n",
+      "2400.0\n",
+      "2500.0\n",
+      "2600.0\n",
+      "2700.0\n",
+      "2800.0\n",
+      "2900.0\n",
+      "3000.0\n",
+      "3100.0\n",
+      "3200.0\n",
+      "3300.0\n",
+      "3400.0\n",
+      "3500.0\n",
+      "3600.0\n",
+      "3700.0\n",
+      "3800.0\n",
+      "3900.0\n",
+      "4000.0\n",
+      "4100.0\n",
+      "4200.0\n",
+      "4300.0\n",
+      "4400.0\n",
+      "4500.0\n",
+      "4600.0\n",
+      "4700.0\n",
+      "4800.0\n",
+      "4900.0\n",
+      "5000.0\n",
+      "5100.0\n",
+      "5200.0\n",
+      "5300.0\n",
+      "5400.0\n",
+      "5500.0\n",
+      "5600.0\n",
+      "5700.0\n",
+      "5800.0\n",
+      "5900.0\n",
+      "6000.0\n",
+      "6100.0\n",
+      "6200.0\n",
+      "6300.0\n",
+      "6400.0\n",
+      "6500.0\n",
+      "6600.0\n",
+      "6700.0\n",
+      "6800.0\n",
+      "6900.0\n",
+      "7000.0\n",
+      "7100.0\n",
+      "7200.0\n",
+      "7300.0\n",
+      "7400.0\n",
+      "7500.0\n",
+      "7600.0\n",
+      "7700.0\n",
+      "7800.0\n",
+      "7900.0\n",
+      "8000.0\n",
+      "8100.0\n",
+      "8200.0\n",
+      "8300.0\n",
+      "8400.0\n",
+      "8500.0\n",
+      "8600.0\n",
+      "8700.0\n",
+      "8800.0\n",
+      "8900.0\n",
+      "9000.0\n",
+      "9100.0\n",
+      "9200.0\n",
+      "9300.0\n",
+      "9400.0\n",
+      "9500.0\n",
+      "9600.0\n",
+      "9700.0\n",
+      "9800.0\n",
+      "9900.0\n"
      ]
     }
    ],
@@ -220,31 +219,31 @@
     "    if np.mod(i,1e4) == 0:\n",
     "        print(t_vec[i])\n",
     "\n",
-    "    if t_vec[i] == 0.4*np.max(t_vec):\n",
-    "        V.influx = 0\n",
+    "    if i == 0.4*(nt+1):\n",
+    "        V.set_influx(0.)\n",
     "\n",
-    "    p = rho*g*V.level-0.5*rho*(V.outflux_vel)**2\n",
-    "\n",
-    "    LA_soll = Pegelregler.get_control_variable(V.level)\n",
-    "    T1.change_LA(LA_soll,dt)\n",
+    "    p = V.get_current_pressure()\n",
+    "    Pegelregler.update_control_variable(V.level)\n",
+    "    LA_soll = Pegelregler.get_current_control_variable()\n",
+    "    T1.update_LA(LA_soll)\n",
+    "    T1.set_pressure(p+offset_pressure)\n",
     "    LA_soll_vec[i] = LA_soll\n",
-    "    LA_ist_vec[i]  = T1.LA\n",
-    "    Q_vec[i] = T1.get_Q(p+offset_pressure)\n",
+    "    LA_ist_vec[i]  = T1.get_current_LA()\n",
+    "    Q_vec[i] = T1.get_current_Q()\n",
     "\n",
-    "    V.pressure = p\n",
-    "    V.outflux_vel = 1/V.area_outflux*Q_vec[i]\n",
+    "    \n",
+    "    V.set_outflux(Q_vec[i])\n",
     "\n",
-    "    V.e_RK_4()                                                              \n",
-    "    V.level = V.update_level(V.timestep)                                    \n",
-    "    V.set_volume()                                                                                                                      \n",
-    "    level_vec[i] = V.level \n",
+    "    V.timestep_reservoir_evolution()                                                              \n",
+    "                                                                                                                   \n",
+    "    level_vec[i] = V.get_current_level()\n",
     "    \n",
     "    "
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 39,
+   "execution_count": 6,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -257,7 +256,7 @@
     "axs1[0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
     "axs1[0].set_ylabel(r'$h$ [$\\mathrm{m}$]')\n",
     "axs1[0].plot(t_vec,level_vec)\n",
-    "axs1[0].set_ylim([0.85*initial_level,1.05*initial_level])\n",
+    "axs1[0].set_ylim([0*initial_level,1.5*initial_level])\n",
     "axs1[1].set_title('Flux')\n",
     "axs1[1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
     "axs1[1].set_ylabel(r'$Q$ [$\\mathrm{m} / \\mathrm{s}^3$]')\n",
@@ -265,51 +264,33 @@
     "axs1[1].set_ylim([0,2*initial_influx])\n",
     "axs1[2].set_title('LA')\n",
     "axs1[2].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
-    "axs1[2].set_ylabel(r'$LA$ [\\%]')\n",
+    "axs1[2].set_ylabel(r'$LA$ [%]')\n",
     "axs1[2].plot(t_vec,LA_soll_vec)\n",
     "axs1[2].plot(t_vec,LA_ist_vec)\n",
     "axs1[2].set_ylim([0,1])\n",
     "fig1.tight_layout()\n",
-    "fig1.show()\n",
-    "plt.pause(1)"
+    "fig1.show()\n"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 40,
+   "execution_count": 7,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "[<matplotlib.lines.Line2D at 0x26263b78be0>]"
+       "[<matplotlib.lines.Line2D at 0x1caf15caca0>]"
       ]
      },
-     "execution_count": 40,
+     "execution_count": 7,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
     "fig2 = plt.figure()\n",
-    "plt.plot(t_vec,Pegelregler.error_history[1:])"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 41,
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "[8.         8.         8.         ... 7.21126138 7.21126138 7.21126138]\n"
-     ]
-    }
-   ],
-   "source": [
-    "print(level_vec[:])"
+    "plt.plot(t_vec,Pegelregler.get_error_history())"
    ]
   }
  ],
diff --git a/Regler/Regler_class_file.py b/Regler/Regler_class_file.py
index c205257..dcbf58a 100644
--- a/Regler/Regler_class_file.py
+++ b/Regler/Regler_class_file.py
@@ -2,13 +2,13 @@ import numpy as np
 #based on https://en.wikipedia.org/wiki/PID_controller#Discrete_implementation
 
 def trap_int(vec,timestep):
+    # numerical integration via the trapeziod rule to calculate the performance parameters
     l = np.size(vec)
     int = 0
     for i in range(l-1):
         int = int + (vec[i]+vec[i+1])/2*timestep
     return int
 
-
 def ISE_fun(error_history,timestep):
     # calcuate the integral of square error
     e = np.array(error_history)
@@ -74,61 +74,50 @@ class P_controller_class:
 
 
 class PI_controller_class:
-    def __init__(self,setpoint,deadband,proportionality_constant,Ti, timestep):
+# init
+    def __init__(self,setpoint,deadband,proportionality_constant,Ti,timestep,lower_limit=0.,upper_limit=1.):
         self.SP = setpoint
         self.db = deadband
         self.Kp = proportionality_constant
-        self.Ti = Ti
+        self.Ti = Ti                        # integration time
         self.dt = timestep
-        self.error_history = [0] 
+        # use a list to be able to append more easily - will get converted to np.array when needed
+        self.error_history = [0]            
 
-        self.cv_lower_limit = 0   # default
-        self.cv_upper_limit = +1   # default
+        self.control_variable = -99
 
+        self.cv_lower_limit = lower_limit   # limits for the controll variable
+        self.cv_upper_limit = upper_limit   # limits for the controll variable
 
-    def set_control_variable_limits(self,lower_limit,upper_limit):
-        self.cv_lower_limit = lower_limit
-        self.cv_upper_limit = upper_limit
+# setter
+    def set_setpoint(self,setpoint):
+        self.SP = setpoint
 
-    def calculate_error(self,process_variable):
-        self.error = process_variable-self.SP
-        self.error_history.append(self.error)
+    def set_control_variable(self,control_variable, display_warning=True):
+        if display_warning == True and self.control_variable != -99:
+            print('WARNING! You are setting the control variable of the PI controller manually \
+                and are not using the .update_controll_variable() method')
+        self.control_variable = control_variable
 
-    def get_control_variable(self,process_variable):
-
-        self.calculate_error(process_variable)
-
-        cv = self.control_variable
-        Kp = self.Kp
-        Ti = self.Ti
-        dt = self.dt
-
-        e0 = self.error_history[-1]
-        e1 = self.error_history[-2] 
-        if abs(self.error) > self.db:
-            new_control = cv+Kp*(e0-e1)+dt/Ti*e0
-        else:
-            new_control = cv
-
-        if new_control < self.cv_lower_limit:
-            new_control = self.cv_lower_limit
-
-        if new_control > self.cv_upper_limit:
-            new_control = self.cv_upper_limit
-        self.control_variable = new_control
+# getter
+    def get_current_control_variable(self):
         return self.control_variable
 
+    def get_error_history(self):
+        return self.error_history[1:]
+
     def get_performance_indicators(self,ISE=True,IAE=True,ITSE=True,ITAE=True):
+        # calculate and return the performance indicators of the error history
         ise     = np.nan
         iae     = np.nan
         itse    = np.nan
         itae    = np.nan
 
         # self.error_history[1:] because the first value of the error history is set to [0]
-        #   to avoid special case handling in the calculation of the controll variable
-        if ISE == True:
+        #   to avoid special case handling in the calculation of the control variable
+        if ISE  == True:
             ise = ISE_fun(self.error_history[1:],self.dt)
-        if IAE == True:
+        if IAE  == True:
             iae = IAE_fun(self.error_history[1:],self.dt)
         if ITSE == True:
             itse = ITSE_fun(self.error_history[1:],self.dt)
@@ -137,4 +126,58 @@ class PI_controller_class:
 
         return ise,iae,itse,itae
 
+    def get_info(self):        
+        new_line = '\n'
+        # :<10 pads the self.value to be 10 characters wide
+        print_str = (f"Turbine has the following attributes: {new_line}" 
+            f"----------------------------- {new_line}"
+            f"Type                      =   PI Controller {new_line}"
+            f"Setpoint                  =   {self.SP:<10} {new_line}"
+            f"Deadband                  =   {self.db:<10} {new_line}"
+            f"Proportionality constant  =   {self.Kp:<10} {new_line}"
+            f"Integration time          =   {self.Ti:<10} [s] {new_line}"
+            f"Current control variable  =   {round(self.control_variable,3):<10}  {new_line}"
+            f"Lower limit CV            =   {self.cv_lower_limit:<10}  {new_line}"
+            f"Upper limit CV            =   {self.cv_upper_limit:<10}  {new_line}"
+            f"Simulation timestep       =   {self.dt:<10} [s] {new_line}"
+            f"----------------------------- {new_line}")
+
+        print(print_str)
+
+# methods
+    def calculate_error(self,process_variable):
+        # calculate the error and expand the err history
+        self.error = process_variable-self.SP
+        self.error_history.append(self.error)
+
+    def update_control_variable(self,process_variable):
+        # calculate the current control variable and make sure it does not exceed the limits
+        self.calculate_error(process_variable)
+
+        # initialize some variables
+        cv = self.control_variable
+        Kp = self.Kp
+        Ti = self.Ti
+        dt = self.dt
+
+        e0 = self.error_history[-1]
+        e1 = self.error_history[-2] 
+
+        # test if the error exceeds the deadband range
+            # only if that is the case, change control variable
+        if abs(self.error) > self.db:
+            new_control = cv+Kp*(e0-e1)+dt/Ti*e0
+        else:
+            new_control = cv
+
+        # ensure that the controll variable stays within the predefined limits
+        if new_control < self.cv_lower_limit:
+            new_control = self.cv_lower_limit
+        if new_control > self.cv_upper_limit:
+            new_control = self.cv_upper_limit
+
+        # set the control variable attribute
+        self.set_control_variable(new_control,display_warning=False)
+
+
     
diff --git a/Turbinen/Turbinen_class_file.py b/Turbinen/Turbinen_class_file.py
index 842b345..e1b3d5a 100644
--- a/Turbinen/Turbinen_class_file.py
+++ b/Turbinen/Turbinen_class_file.py
@@ -1,3 +1,4 @@
+from time import time
 import numpy as np
 #importing pressure conversion function
 import sys
@@ -8,35 +9,117 @@ sys.path.append(parent)
 from functions.pressure_conversion import pressure_conversion
 
 class Francis_Turbine:
-    def __init__(self, Q_nenn,p_nenn):
-        self.Q_n    = Q_nenn
-        self.p_n    = p_nenn
-        self.LA_n   = 1. # 100%
-        h           = pressure_conversion(p_nenn,'Pa','MWs')
-        self.A      = Q_nenn/(np.sqrt(2*9.81*h)*0.98)
+ # units 
+    # make sure that units and print units are the same
+    # units are used to label graphs and print units are used to have a bearable format when using pythons print()
+    density_unit        = r'$\mathrm{kg}/\mathrm{m}^3$'
+    flux_unit           = r'$\mathrm{m}^3/\mathrm{s}$'
+    LA_unit             = '%'
+    pressure_unit       = 'Pa'
+    time_unit           = 's'
+    velocity_unit       = r'$\mathrm{m}/\mathrm{s}$'
+    volume_unit         = r'$\mathrm{m}^3$'
 
-    def set_LA(self,LA):
+    density_unit_print      = 'kg/m³'
+    flux_unit_print         = 'm³/s'
+    LA_unit_print           = '%'
+    pressure_unit_print     = 'mWS' 
+    time_unit_print         = 's'
+    velocity_unit_print     = 'm/s'
+    volume_unit_print       = 'm³'
+
+    g = 9.81    # m/s²  gravitational acceleration
+
+ # init
+    def __init__(self, Q_nenn,p_nenn,t_closing=-1.,timestep=-1.):
+        self.Q_n    = Q_nenn                                    # nominal flux 
+        self.p_n    = p_nenn                                    # nominal pressure
+        self.LA_n   = 1. # 100%                                 # nominal Leitapparatöffnung
+        h           = pressure_conversion(p_nenn,'Pa','MWs')    # nominal pressure in terms of hydraulic head
+        self.A      = Q_nenn/(np.sqrt(2*self.g*h)*0.98)         # Ersatzfläche
+
+        self.dt     = timestep                                  # simulation timestep
+        self.t_c = t_closing                                    # closing time
+        self.d_LA_max_dt = 1/t_closing                          # maximal change of LA per second
+
+        # initialize for get_info() - parameters will be converted to display -1 if not overwritten
+        self.p  = pressure_conversion(-1,self.pressure_unit_print,self.pressure_unit)
+        self.Q  = -1.
+        self.LA = -0.01 
+
+
+# setter
+    def set_LA(self,LA,display_warning=True):
+        # set Leitapparatöffnung
         self.LA = LA
+        # warn user, that the .set_LA() method should not be used ot set LA manually
+        if display_warning == True:
+            print('Consider using the .update_LA() method instead of setting LA manually')
+
+    def set_timestep(self,timestep,display_warning=True):
+        # set Leitapparatöffnung
+        self.dt = time
+        # warn user, that the .set_LA() method should not be used ot set LA manually
+        if display_warning == True:
+            print('WARNING: You are changing the timestep of the turbine simulation. This has implications on the simulated closing speed!')
+
     def set_pressure(self,pressure):
+        # set pressure in front of the turbine
         self.p = pressure
 
-    def get_Q(self):
+#getter
+    def get_current_Q(self):
+        # return the flux through the turbine, based on the current pressure in front
+            #  of the turbine and the Leitapparatöffnung
         self.Q = self.Q_n*(self.LA/self.LA_n)*np.sqrt(self.p/self.p_n)
         return self.Q
 
-    def set_closing_time(self,t_closing):
-        self.t_c = t_closing
-        self.d_LA_max_dt = 1/t_closing
+    def get_current_LA(self):
+        return self.LA
 
-    def change_LA(self,LA_soll,timestep):
-        LA_diff = self.LA-LA_soll
-        LA_diff_max = self.d_LA_max_dt*timestep
-        if abs(LA_diff) > LA_diff_max:
-            LA_diff = np.sign(LA_diff)*LA_diff_max
-        self.LA = self.LA-LA_diff
+    def get_info(self, full = False):
+        new_line = '\n'
+        p   = pressure_conversion(self.p,self.pressure_unit,self.pressure_unit_print)
+        p_n = pressure_conversion(self.p_n,self.pressure_unit,self.pressure_unit_print)
+
+        
+        if full == True:
+            # :<10 pads the self.value to be 10 characters wide
+            print_str = (f"Turbine has the following attributes: {new_line}" 
+                f"----------------------------- {new_line}"
+                f"Type                  =       Francis {new_line}"
+                f"Nominal flux          =       {self.Q_n:<10} {self.flux_unit_print} {new_line}"
+                f"Nominal pressure      =       {round(p_n,3):<10} {self.pressure_unit_print}{new_line}"
+                f"Nominal LA            =       {self.LA_n*100:<10} {self.LA_unit_print} {new_line}"
+                f"Closing time          =       {self.t_c:<10} {self.time_unit_print} {new_line}"
+                f"Current flux          =       {self.Q:<10} {self.flux_unit_print} {new_line}" 
+                f"Current pipe pressure =       {round(p,3):<10} {self.pressure_unit_print} {new_line}"
+                f"Current LA            =       {self.LA*100:<10} {self.LA_unit_print} {new_line}"
+                f"Simulation timestep   =       {self.dt:<10} {self.time_unit_print} {new_line}"
+                f"----------------------------- {new_line}")
+        else:
+            # :<10 pads the self.value to be 10 characters wide
+            print_str = (f"The current attributes are: {new_line}" 
+                f"----------------------------- {new_line}"
+                f"Current flux          =       {self.Q:<10} {self.flux_unit_print} {new_line}" 
+                f"Current pipe pressure =       {round(p,3):<10} {self.pressure_unit_print} {new_line}"
+                f"Current LA            =       {self.LA*100:<10} {self.LA_unit_print} {new_line}"
+                f"----------------------------- {new_line}")
+
+        print(print_str)
+
+# methods
+    def update_LA(self,LA_soll):
+        # update the Leitappartöffnung and consider the restrictions of the closing time of the turbine
+        LA_diff = self.LA-LA_soll                   # calculate the difference to the target LA
+        LA_diff_max = self.d_LA_max_dt*self.dt     # calculate the maximum change in LA based on the given timestep  
+        LA_diff = np.sign(LA_diff)*np.min(np.abs([LA_diff,LA_diff_max]))    # calulate the correct change in LA
+        self.set_LA(self.LA-LA_diff,display_warning=False)                # set new LA
 
     def set_steady_state(self,ss_flux,ss_pressure):
+        # calculate and set steady state LA, that allows the flow of ss_flux at ss_pressure through the
+            # turbine at the steady state LA
         ss_LA = self.LA_n*ss_flux/self.Q_n*np.sqrt(self.p_n/ss_pressure)
-        self.set_LA(ss_LA)
         if ss_LA < 0 or ss_LA > 1:
-            print('LA out of range')
+            raise Exception('LA out of range [0;1]')
+        self.set_LA(ss_LA,display_warning=False)
diff --git a/Turbinen/Turbinen_test.ipynb b/Turbinen/Turbinen_test_Kennlinie.ipynb
similarity index 99%
rename from Turbinen/Turbinen_test.ipynb
rename to Turbinen/Turbinen_test_Kennlinie.ipynb
index 7cfe0c3..5fdb785 100644
--- a/Turbinen/Turbinen_test.ipynb
+++ b/Turbinen/Turbinen_test_Kennlinie.ipynb
@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 4,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -23,18 +23,39 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 5,
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Turbine has the following attributes: \n",
+      "----------------------------- \n",
+      "Type                  =       Francis \n",
+      "Nominal flux          =       0.85       m³/s \n",
+      "Nominal pressure      =       108.09     mWS\n",
+      "Nominal LA            =       100.0      % \n",
+      "Closing time          =       -1         s \n",
+      "Current flux          =       -1         m³/s \n",
+      "Current pipe pressure =       -1.0       mWS \n",
+      "Current LA            =       -1.0       % \n",
+      "Simulation timestep   =       -1         s \n",
+      "----------------------------- \n",
+      "\n"
+     ]
+    }
+   ],
    "source": [
     "Q_nenn = 0.85\n",
     "p_nenn = pressure_conversion(10.6,'bar','Pa')\n",
-    "Untertweng1 = Francis_Turbine(Q_nenn,p_nenn)"
+    "Untertweng1 = Francis_Turbine(Q_nenn,p_nenn)\n",
+    "Untertweng1.get_info(full=True)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 6,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -47,7 +68,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 20,
+   "execution_count": 7,
    "metadata": {},
    "outputs": [
     {
@@ -56,14 +77,14 @@
        "Text(0.5, 0, 'Q [m³/s]')"
       ]
      },
-     "execution_count": 20,
+     "execution_count": 7,
      "metadata": {},
      "output_type": "execute_result"
     },
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "7650cd2a1f9a4b87b98e70add29f11a7",
+       "model_id": "75adb3cb47e642e3a5606cb41efedf72",
        "version_major": 2,
        "version_minor": 0
       },
@@ -94,8 +115,8 @@
     "for i in range(n_p):\n",
     "    for j in range(n_LA):\n",
     "        Untertweng1.set_pressure(pp[i,j])\n",
-    "        Untertweng1.set_LA(ll[i,j])\n",
-    "        Q_mat[i,j] = Untertweng1.get_Q()\n",
+    "        Untertweng1.set_LA(ll[i,j],display_warning=False)\n",
+    "        Q_mat[i,j] = Untertweng1.get_current_Q()\n",
     "\n",
     "fig1 = plt.figure()\n",
     "ax1 = plt.axes(projection='3d')\n",
@@ -108,31 +129,23 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 27,
+   "execution_count": 8,
    "metadata": {},
    "outputs": [
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "C:\\Users\\BRANT\\AppData\\Local\\Temp\\9\\ipykernel_7508\\1599598770.py:5: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
-      "  fig = plt.figure()\n"
-     ]
-    },
     {
      "data": {
       "text/plain": [
        "(0.0, 1.275)"
       ]
      },
-     "execution_count": 27,
+     "execution_count": 8,
      "metadata": {},
      "output_type": "execute_result"
     },
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "52152d0b96c74a4ebfb18041a22d8d0e",
+       "model_id": "224f00f9bf85446b845685a08ed27c68",
        "version_major": 2,
        "version_minor": 0
       },
@@ -170,7 +183,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 30,
+   "execution_count": 9,
    "metadata": {},
    "outputs": [
     {
@@ -179,14 +192,14 @@
        "(0.0, 1.275)"
       ]
      },
-     "execution_count": 30,
+     "execution_count": 9,
      "metadata": {},
      "output_type": "execute_result"
     },
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "bf6b0fe84d264693813f6e991600ece1",
+       "model_id": "92741b823c9749c9820ee7b5ba47a6bc",
        "version_major": 2,
        "version_minor": 0
       },
@@ -221,23 +234,6 @@
     "plt.title('P = '+ str(p_test2) + ' [Pa]')\n",
     "plt.ylim([0,1.5*Q_nenn])"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 14,
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "53000.0\n"
-     ]
-    }
-   ],
-   "source": [
-    "print(pp[10,5])"
-   ]
   }
  ],
  "metadata": {
diff --git a/Turbinen/messy.ipynb b/Turbinen/messy.ipynb
new file mode 100644
index 0000000..d9e0b09
--- /dev/null
+++ b/Turbinen/messy.ipynb
@@ -0,0 +1,90 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "from Turbinen_class_file import Francis_Turbine\n",
+    "from mpl_toolkits import mplot3d\n",
+    "import matplotlib.pyplot as plt\n",
+    "%matplotlib widget\n",
+    "\n",
+    "#importing pressure conversion function\n",
+    "import sys\n",
+    "import os\n",
+    "current = os.path.dirname(os.path.realpath('messy.ipynb'))\n",
+    "parent = os.path.dirname(current)\n",
+    "sys.path.append(parent)\n",
+    "from functions.pressure_conversion import pressure_conversion"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "The current attributes are: \n",
+      "----------------------------- \n",
+      "Current flux          =       -1.0       m³/s \n",
+      "Current pipe pressure =       -1.0       mWS \n",
+      "Current LA            =       -1.0       % \n",
+      "----------------------------- \n",
+      "\n",
+      "The current attributes are: \n",
+      "----------------------------- \n",
+      "Current flux          =       -1.0       m³/s \n",
+      "Current pipe pressure =       -1.0       mWS \n",
+      "Current LA            =       -1.0       % \n",
+      "----------------------------- \n",
+      "\n"
+     ]
+    }
+   ],
+   "source": [
+    "Q_nenn = 0.85\n",
+    "p_nenn = pressure_conversion(10.6,'bar','Pa')\n",
+    "Untertweng1 = Francis_Turbine(Q_nenn,p_nenn)\n",
+    "Untertweng2 = Francis_Turbine(Q_nenn,p_nenn)\n",
+    "\n",
+    "\n",
+    "turbines = [Untertweng1,Untertweng2]\n",
+    "for turbine in turbines:\n",
+    "    turbine.get_info()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3.8.13 ('Georg_DT_Slot3')",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.13"
+  },
+  "orig_nbformat": 4,
+  "vscode": {
+   "interpreter": {
+    "hash": "84fb123bdc47ab647d3782661abcbe80fbb79236dd2f8adf4cef30e8755eb2cd"
+   }
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/Untertweng.ipynb b/Untertweng.ipynb
index b2baf53..9eedd6d 100644
--- a/Untertweng.ipynb
+++ b/Untertweng.ipynb
@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -17,15 +17,16 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
     "#define constants\n",
     "\n",
     "#Turbine\n",
-    "Q_nenn = 0.85\n",
-    "p_nenn = pressure_conversion(10.6,'bar','Pa')\n",
+    "Q_nenn          = 0.85  # m³/s\n",
+    "p_nenn          = pressure_conversion(10.6,'bar','Pa')\n",
+    "closing_time    = 30    #s\n",
     "\n",
     "# physics\n",
     "g               = 9.81                          # gravitational acceleration [m/s²]\n",
@@ -37,35 +38,31 @@
     "A_pipe          = D**2/4*np.pi                  # pipeline area\n",
     "h_pipe          = 105                           # hydraulic head without reservoir [m] \n",
     "alpha           = np.arcsin(h_pipe/L)           # Höhenwinkel der Druckrohrleitung \n",
-    "n               = 50                            # number of pipe segments in discretization\n",
+    "n               = 200                            # number of pipe segments in discretization\n",
     "# consider replacing Q0 with a vector be be more flexible in initial conditions\n",
     "# Q0              = Q_nenn                        # initial flow in whole pipe [m³/s]\n",
     "# v0              = Q0/A_pipe                     # initial flow velocity [m/s]\n",
     "f_D             = 0.014                           # Darcy friction factor\n",
     "c               = 500.                          # propagation velocity of the pressure wave [m/s]\n",
     "# consider prescribing a total simulation time and deducting the number of timesteps from that\n",
-    "nt              = 3000                          # number of time steps after initial conditions\n",
+    "nt              = 20000                          # number of time steps after initial conditions\n",
     "\n",
     "# derivatives of the pipeline constants\n",
     "dx              = L/n                           # length of each pipe segment\n",
     "dt              = dx/c                          # timestep according to method of characterisitics\n",
     "nn              = n+1                           # number of nodes\n",
     "initial_level   = 8.                            # water level in upstream reservoir [m]\n",
-    "# p0              = rho*g*initial_level-v0**2*rho/2\n",
-    "pl_vec          = np.arange(0,nn*dx,dx)         # pl = pipe-length. position of the nodes on the pipeline\n",
+    "pl_vec          = np.arange(0,nn,1)*dx          # pl = pipe-length. position of the nodes on the pipeline\n",
     "t_vec           = np.arange(0,nt+1)*dt          # time vector\n",
-    "h_vec           = np.arange(0,n+1)*h_pipe/n     # hydraulic head of pipeline at each node \n",
-    "# v_init          = np.full(nn,Q0/(D**2/4*np.pi)) # initial velocity distribution in pipeline\n",
-    "# p_init          = (rho*g*(initial_level+h_vec)-v_init**2*rho/2)-(f_D*pl_vec/D*rho/2*v_init**2) # ref Wikipedia: Darcy Weisbach\n",
+    "h_vec           = np.arange(0,nn,1)*h_pipe/n    # hydraulic head of pipeline at each node \n",
+    "\n",
     "\n",
     "\n",
     "# reservoir\n",
     "# replace influx by vector\n",
-    "initial_influx              = Q_nenn/1.1        # initial influx  of volume to the reservoir [m³/s]\n",
-    "# initial_outflux             = Q0        # initial outflux of volume from the reservoir to the pipeline [m³/s]\n",
-    "# initial_pipeline_pressure   = p0        # Initial condition for the static pipeline pressure at the reservoir (= hydrostatic pressure - dynamic pressure) \n",
+    "initial_flux                = Q_nenn/1.1        # initial influx  of volume to the reservoir [m³/s]\n",
     "initial_pressure_unit       = 'Pa'      # DO NOT CHANGE! for pressure conversion in print statements and plot labels \n",
-    "conversion_pressure_unit    = 'bar'      # for pressure conversion in print statements and plot labels\n",
+    "conversion_pressure_unit    = 'bar'     # for pressure conversion in print statements and plot labels\n",
     "area_base                   = 74.       # total base are of the cuboid reservoir [m²]   \n",
     "area_outflux                = A_pipe    # outlfux area of the reservoir, given by pipeline area [m²]\n",
     "critical_level_low          = 0.        # for yet-to-be-implemented warnings[m]\n",
@@ -78,43 +75,28 @@
     "\n"
    ]
   },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "#### Ideas for checks after constant definitions: \n",
-    "\n",
-    "- Check that the initial pressure is not negative:\n",
-    "    - may happen, if there is too little hydraulic head to create the initial flow conditions with the given friction\n",
-    "<br>\n",
-    "<br>\n",
-    "- plausbility checks?\n",
-    "    - area > area_outflux ?\n",
-    "    - propable ranges for parameters?\n",
-    "    - angle and height/length fit together?\n",
-    "    "
-   ]
-  },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
     "# create objects\n",
     "\n",
     "V = Ausgleichsbecken_class(area_base,area_outflux,critical_level_low,critical_level_high,simulation_timestep)\n",
-    "V.set_steady_state(initial_influx,initial_level,conversion_pressure_unit)\n",
+    "V.set_steady_state(initial_flux,initial_level,conversion_pressure_unit)\n",
+    "\n",
     "\n",
     "pipe = Druckrohrleitung_class(L,D,n,alpha,f_D)\n",
     "pipe.set_pressure_propagation_velocity(c)\n",
     "pipe.set_number_of_timesteps(nt)\n",
-    "pipe.set_steady_state(initial_influx,V.level,pl_vec,h_vec)\n",
+    "pipe.set_steady_state(initial_flux,initial_level,pl_vec,h_vec)\n",
     "\n",
+    "initial_pressure_turbine = pipe.get_current_pressure_distribution()[-1]\n",
+    "\n",
+    "T1 = Francis_Turbine(Q_nenn,p_nenn,closing_time,timestep = dt)\n",
+    "T1.set_steady_state(initial_flux,initial_pressure_turbine)\n",
     "\n",
-    "T1 = Francis_Turbine(Q_nenn,p_nenn)\n",
-    "T1.set_steady_state(initial_influx,pipe.p0[-1])\n",
-    "T1.set_closing_time(30)\n",
     "\n",
     "# display the attributes of the created reservoir and pipeline object\n",
     "# V.get_info(full=True)\n",
@@ -123,19 +105,19 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
     "# initialization for timeloop\n",
     "\n",
     "# prepare the vectors in which the pressure and velocity distribution in the pipeline from the previous timestep are stored\n",
-    "v_old = pipe.v0.copy()\n",
-    "p_old = pipe.p0.copy()\n",
+    "v_old = pipe.get_current_velocity_distribution()\n",
+    "p_old = pipe.get_current_pressure_distribution()\n",
     "\n",
     "# prepare the vectors in which the temporal evolution of the boundary conditions are stored\n",
-    "    # keep in mind, that the velocity at the turbine and the pressure at the reservoir are set manually and\n",
-    "        # through the time evolution of the reservoir respectively  \n",
+    "    # keep in mind, that the velocity at the turbine and the pressure at the reservoir follow from boundary conditions\n",
+    "        # reservoir level and flow through turbine\n",
     "    # the pressure at the turbine and the velocity at the reservoir are calculated from the method of characteristics\n",
     "v_boundary_res  = np.zeros_like(t_vec)\n",
     "v_boundary_tur  = np.zeros_like(t_vec)\n",
@@ -143,8 +125,7 @@
     "p_boundary_tur  = np.zeros_like(t_vec)\n",
     "\n",
     "# prepare the vectors that store the temporal evolution of the level in the reservoir\n",
-    "level_vec       = np.full(nt+1,V.level)  # level at the end of each pipeline timestep\n",
-    "level_vec_2     = np.zeros([nt_eRK4])   # level throughout each reservoir timestep-used for plotting and overwritten afterwards\n",
+    "level_vec       = np.full(nt+1,initial_level)  # level at the end of each pipeline timestep\n",
     "\n",
     "# set the boundary conditions for the first timestep\n",
     "v_boundary_res[0]   = v_old[0]\n",
@@ -153,14 +134,15 @@
     "p_boundary_tur[0]   = p_old[-1]\n",
     "\n",
     "LA_soll_vec = np.full_like(t_vec,T1.LA)\n",
-    "LA_soll_vec[1500:]= 0\n",
+    "LA_soll_vec[500:]= 0\n",
+    "LA_ist_vec  = np.full_like(t_vec,T1.LA)\n",
     "\n",
     "\n"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -176,16 +158,11 @@
     "axs1[0].set_ylabel(r'$p$ ['+conversion_pressure_unit+']')\n",
     "axs1[1].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
     "axs1[1].set_ylabel(r'$v$ [$\\mathrm{m} / \\mathrm{s}$]')\n",
-    "lo_00, = axs1[0].plot(pl_vec,pressure_conversion(pipe.p_old,initial_pressure_unit, conversion_pressure_unit),marker='.')\n",
-    "lo_01, = axs1[1].plot(pl_vec,pipe.v_old,marker='.')\n",
+    "lo_00, = axs1[0].plot(pl_vec,pressure_conversion(p_old,initial_pressure_unit, conversion_pressure_unit),marker='.')\n",
+    "lo_01, = axs1[1].plot(pl_vec,v_old,marker='.')\n",
     "axs1[0].autoscale()\n",
     "axs1[1].autoscale()\n",
-    "# displaying the reservoir level within each pipeline timestep\n",
-    "# axs1[2].set_title('Level reservoir')\n",
-    "# axs1[2].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
-    "# axs1[2].set_ylabel(r'$h$ [m]')\n",
-    "# lo_02, = axs1[2].plot(level_vec_2)\n",
-    "# axs1[2].autoscale()\n",
+    "\n",
     "fig1.tight_layout()\n",
     "fig1.show()\n",
     "plt.pause(1)\n"
@@ -193,60 +170,62 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
     "# loop through time steps of the pipeline\n",
     "for it_pipe in range(1,pipe.nt+1):\n",
     "\n",
+    "    if t_vec[it_pipe]>20:\n",
+    "        if V.get_current_influx() > 0:\n",
+    "            V.set_influx(np.max([V.get_current_influx()-initial_flux*5*1e-3,0.]))\n",
+    "\n",
     "# for each pipeline timestep, execute nt_eRK4 timesteps of the reservoir code\n",
     "    # set initial conditions for the reservoir time evolution calculted with e-RK4\n",
-    "    V.pressure      = p_old[0]\n",
-    "    V.outflux_vel   = v_old[0]\n",
+    "    V.set_pressure      = p_old[0]\n",
+    "    V.set_outflux       = v_old[0]*area_outflux\n",
     "    # calculate the time evolution of the reservoir level within each pipeline timestep to avoid runaway numerical error\n",
     "    for it_res in range(nt_eRK4):\n",
-    "        V.e_RK_4()                                                              # call e-RK4 to update outflux\n",
-    "        V.level = V.update_level(V.timestep)                                    # \n",
-    "        V.update_volume()                                                          # update volume in reservoir\n",
-    "        level_vec_2[it_res] = V.level                                           # save for plotting\n",
-    "        if (V.level < critical_level_low) or (V.level > critical_level_high):   # make sure to never exceed critical levels\n",
-    "            i_max = it_pipe                                                     # for plotting only calculated values\n",
-    "            break                                                               \n",
-    "    level_vec[it_pipe] = V.level                                                  \n",
+    "        V.timestep_reservoir_evolution()                                                             \n",
+    "    level_vec[it_pipe] = V.get_current_level()                                                 \n",
+    "\n",
+    "    # change the Leitapparatöffnung  based on the target value\n",
+    "    T1.update_LA(LA_soll_vec[it_pipe])\n",
+    "    T1.set_pressure(p_old[-1])\n",
+    "\n",
+    "    LA_ist_vec[it_pipe] = T1.LA\n",
     "\n",
     "    # set boundary conditions for the next timestep of the characteristic method\n",
-    "    p_boundary_res[it_pipe] = rho*g*V.level-V.outflux_vel**2*rho/2\n",
-    "\n",
-    "    T1.change_LA(LA_soll_vec[it_pipe],dt)\n",
-    "    v_boundary_tur[it_pipe] = 1/A_pipe*T1.get_Q(p_old[-1])\n",
+    "    p_boundary_res[it_pipe] = V.get_current_pressure()\n",
+    "    v_boundary_tur[it_pipe] = 1/A_pipe*T1.get_current_Q()\n",
     "\n",
     "    # the the boundary conditions in the pipe.object and thereby calculate boundary pressure at turbine\n",
     "    pipe.set_boundary_conditions_next_timestep(p_boundary_res[it_pipe],v_boundary_tur[it_pipe])\n",
-    "    p_boundary_tur[it_pipe] = pipe.p_boundary_tur\n",
-    "    v_boundary_res[it_pipe] = pipe.v_boundary_res\n",
+    "    p_boundary_tur[it_pipe] = pipe.get_current_pressure_distribution()[-1]\n",
+    "    v_boundary_res[it_pipe] = pipe.get_current_velocity_distribution()[0]\n",
     "\n",
     "    # perform the next timestep via the characteristic method\n",
     "    pipe.timestep_characteristic_method()\n",
     "\n",
+    "    # prepare for next loop\n",
+    "    p_old = pipe.get_current_pressure_distribution()\n",
+    "    v_old = pipe.get_current_velocity_distribution()\n",
+    "\n",
     "    # plot some stuff\n",
     "        # remove line-objects to autoscale axes (there is definetly a better way, but this works ¯\\_(ツ)_/¯ )\n",
     "    lo_00.remove()\n",
     "    lo_01.remove()\n",
     "    # lo_02.remove()\n",
     "        # plot new pressure and velocity distribution in the pipeline\n",
-    "    lo_00, = axs1[0].plot(pl_vec,pressure_conversion(pipe.p_old,initial_pressure_unit, conversion_pressure_unit),marker='.',c='blue')\n",
-    "    lo_01, = axs1[1].plot(pl_vec,pipe.v_old,marker='.',c='blue')\n",
+    "    lo_00, = axs1[0].plot(pl_vec,pressure_conversion(p_old,initial_pressure_unit, conversion_pressure_unit),marker='.',c='blue')\n",
+    "    lo_01, = axs1[1].plot(pl_vec,v_old,marker='.',c='blue')\n",
     "    # lo_02, = axs1[2].plot(level_vec_2,c='blue')\n",
     "    fig1.suptitle(str(round(t_vec[it_pipe],2))+ ' s / '+str(round(t_vec[-1],2)) + ' s' )\n",
     "    fig1.canvas.draw()\n",
     "    fig1.tight_layout()\n",
     "    fig1.show()\n",
-    "    plt.pause(0.1)    \n",
-    "\n",
-    "    # prepare for next loop\n",
-    "    p_old = pipe.p_old\n",
-    "    v_old = pipe.v_old  \n",
+    "    plt.pause(0.001)    \n",
     "\n",
     "        \n",
     "        "
@@ -254,35 +233,44 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 9,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
     "# plot time evolution of boundary pressure and velocity as well as the reservoir level\n",
     "\n",
     "fig2,axs2 = plt.subplots(3,2)\n",
-    "axs2[0,0].plot(t_vec,pressure_conversion(p_boundary_res,initial_pressure_unit, conversion_pressure_unit))\n",
-    "axs2[0,1].plot(t_vec,v_boundary_res)\n",
-    "axs2[0,1].set_ylim(-2*Q_nenn,+2*Q_nenn)\n",
-    "axs2[1,0].plot(t_vec,pressure_conversion(p_boundary_tur,initial_pressure_unit, conversion_pressure_unit))\n",
-    "axs2[1,1].plot(t_vec,v_boundary_tur)\n",
-    "axs2[2,0].plot(t_vec,level_vec)\n",
     "axs2[0,0].set_title('Pressure reservoir')\n",
-    "axs2[0,1].set_title('Velocity reservoir')\n",
-    "axs2[1,0].set_title('Pressure turbine')\n",
-    "axs2[1,1].set_title('Velocity turbine')\n",
-    "axs2[2,0].set_title('Level reservoir')\n",
+    "axs2[0,0].plot(t_vec,pressure_conversion(p_boundary_res,initial_pressure_unit, conversion_pressure_unit))\n",
     "axs2[0,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
     "axs2[0,0].set_ylabel(r'$p$ ['+conversion_pressure_unit+']')\n",
+    "\n",
+    "axs2[0,1].set_title('Velocity reservoir')\n",
+    "axs2[0,1].plot(t_vec,v_boundary_res)\n",
+    "axs2[0,1].set_ylim(-2*Q_nenn,+2*Q_nenn)\n",
     "axs2[0,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
     "axs2[0,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n",
+    "\n",
+    "axs2[1,0].set_title('Pressure turbine')\n",
+    "axs2[1,0].plot(t_vec,pressure_conversion(p_boundary_tur,initial_pressure_unit, conversion_pressure_unit))\n",
     "axs2[1,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
     "axs2[1,0].set_ylabel(r'$p$ ['+conversion_pressure_unit+']')\n",
+    "\n",
+    "axs2[1,1].set_title('Velocity turbine')\n",
+    "axs2[1,1].plot(t_vec,v_boundary_tur)\n",
     "axs2[1,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
     "axs2[1,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n",
+    "\n",
+    "axs2[2,0].set_title('Level reservoir')\n",
+    "axs2[2,0].plot(t_vec,level_vec)\n",
     "axs2[2,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
     "axs2[2,0].set_ylabel(r'$h$ [m]')\n",
-    "axs2[2,1].axis('off')\n",
+    "\n",
+    "axs2[2,1].set_title('LA')\n",
+    "axs2[2,1].plot(t_vec,100*LA_soll_vec)\n",
+    "axs2[2,1].plot(t_vec,100*LA_ist_vec)\n",
+    "axs2[2,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
+    "axs2[2,1].set_ylabel(r'$LA$ [%]')\n",
     "fig2.tight_layout()\n",
     "plt.show()"
    ]
diff --git a/Untertweng_mit_Pegelregler.ipynb b/Untertweng_mit_Pegelregler.ipynb
index 0d5d5e7..8c4669e 100644
--- a/Untertweng_mit_Pegelregler.ipynb
+++ b/Untertweng_mit_Pegelregler.ipynb
@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 22,
+   "execution_count": 18,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -18,20 +18,29 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 23,
+   "execution_count": 19,
    "metadata": {},
    "outputs": [],
    "source": [
     "#define constants\n",
     "\n",
     "#Turbine\n",
-    "Q_nenn = 0.85\n",
-    "p_nenn,_ = pressure_conversion(10.6,'bar','Pa')\n",
+    "Q_nenn          = 0.85  # m³/s\n",
+    "p_nenn          = pressure_conversion(10.6,'bar','Pa')\n",
+    "closing_time    = 70    #s\n",
     "\n",
     "# physics\n",
     "g               = 9.81                          # gravitational acceleration [m/s²]\n",
     "rho             = 1000.                         # density of water [kg/m³]\n",
     "\n",
+    "\n",
+    "# define controller constants\n",
+    "target_level = 8.  # m\n",
+    "Kp = 0.1\n",
+    "Ti = 10.\n",
+    "deadband_range = 0.05 # m\n",
+    "\n",
+    "\n",
     "# pipeline\n",
     "L               = 535.+478.                     # length of pipeline [m]\n",
     "D               = 0.9                           # pipe diameter [m]\n",
@@ -42,34 +51,29 @@
     "f_D             = 0.014                           # Darcy friction factor\n",
     "c               = 500.                          # propagation velocity of the pressure wave [m/s]\n",
     "# consider prescribing a total simulation time and deducting the number of timesteps from that\n",
-    "nt              = 1500                          # number of time steps after initial conditions\n",
+    "nt              = 1000                          # number of time steps after initial conditions\n",
     "\n",
     "# derivatives of the pipeline constants\n",
     "dx              = L/n                           # length of each pipe segment\n",
     "dt              = dx/c                          # timestep according to method of characterisitics\n",
     "nn              = n+1                           # number of nodes\n",
-    "initial_level   = 8.                            # water level in upstream reservoir [m]\n",
-    "pl_vec          = np.arange(0,nn*dx,dx)         # pl = pipe-length. position of the nodes on the pipeline\n",
+    "initial_level   = target_level                  # water level in upstream reservoir [m]\n",
+    "pl_vec          = np.arange(0,nn,1)*dx          # pl = pipe-length. position of the nodes on the pipeline\n",
     "t_vec           = np.arange(0,nt+1)*dt          # time vector\n",
-    "h_vec           = np.arange(0,n+1)*h_pipe/n     # hydraulic head of pipeline at each node \n",
+    "h_vec           = np.arange(0,nn,1)*h_pipe/n    # hydraulic head of pipeline at each node \n",
+    "\n",
+    "\n",
     "\n",
     "# reservoir\n",
     "# replace influx by vector\n",
-    "initial_influx              = Q_nenn/1.1        # initial influx  of volume to the reservoir [m³/s]\n",
+    "initial_flux                = Q_nenn/1.1        # initial influx  of volume to the reservoir [m³/s]\n",
     "initial_pressure_unit       = 'Pa'      # DO NOT CHANGE! for pressure conversion in print statements and plot labels \n",
-    "conversion_pressure_unit    = 'bar'      # for pressure conversion in print statements and plot labels\n",
+    "conversion_pressure_unit    = 'bar'     # for pressure conversion in print statements and plot labels\n",
     "area_base                   = 74.       # total base are of the cuboid reservoir [m²]   \n",
     "area_outflux                = A_pipe    # outlfux area of the reservoir, given by pipeline area [m²]\n",
     "critical_level_low          = 0.        # for yet-to-be-implemented warnings[m]\n",
     "critical_level_high         = np.inf    # for yet-to-be-implemented warnings[m]\n",
     "\n",
-    "\n",
-    "# define controller constants\n",
-    "target_level = initial_level  # m\n",
-    "Kp = 2\n",
-    "Ti = 10\n",
-    "deadband_range = 0.05 # m\n",
-    "\n",
     "# make sure e-RK4 method of reservoir has a small enough timestep to avoid runaway numerical error\n",
     "nt_eRK4                     = 1000      # number of simulation steps of reservoir in between timesteps of pipeline              \n",
     "simulation_timestep         = dt/nt_eRK4\n",
@@ -77,75 +81,55 @@
     "\n"
    ]
   },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "#### Ideas for checks after constant definitions: \n",
-    "\n",
-    "- Check that the initial pressure is not negative:\n",
-    "    - may happen, if there is too little hydraulic head to create the initial flow conditions with the given friction\n",
-    "<br>\n",
-    "<br>\n",
-    "- plausbility checks?\n",
-    "    - area > area_outflux ?\n",
-    "    - propable ranges for parameters?\n",
-    "    - angle and height/length fit together?\n",
-    "    "
-   ]
-  },
   {
    "cell_type": "code",
-   "execution_count": 24,
+   "execution_count": 20,
    "metadata": {},
    "outputs": [],
    "source": [
     "# create objects\n",
     "\n",
     "V = Ausgleichsbecken_class(area_base,area_outflux,critical_level_low,critical_level_high,simulation_timestep)\n",
-    "V.set_steady_state(initial_influx,initial_level,initial_pressure_unit,conversion_pressure_unit)\n",
+    "V.set_steady_state(initial_flux,initial_level,conversion_pressure_unit)\n",
+    "\n",
     "\n",
     "pipe = Druckrohrleitung_class(L,D,n,alpha,f_D)\n",
     "pipe.set_pressure_propagation_velocity(c)\n",
     "pipe.set_number_of_timesteps(nt)\n",
-    "pipe.set_steady_state(initial_influx,V.level,pl_vec,h_vec,initial_pressure_unit,conversion_pressure_unit)\n",
+    "pipe.set_steady_state(initial_flux,initial_level,pl_vec,h_vec)\n",
     "\n",
+    "initial_pressure_turbine = pipe.get_current_pressure_distribution()[-1]\n",
     "\n",
-    "T1 = Francis_Turbine(Q_nenn,p_nenn)\n",
-    "T1.set_steady_state(initial_influx,pipe.p0[-1])\n",
-    "T1.set_closing_time(5)\n",
+    "T1 = Francis_Turbine(Q_nenn,p_nenn,closing_time,timestep = dt)\n",
+    "T1.set_steady_state(initial_flux,initial_pressure_turbine)\n",
     "\n",
     "Pegelregler = PI_controller_class(target_level,deadband_range,Kp,Ti,dt)\n",
-    "\n",
-    "# display the attributes of the created reservoir and pipeline object\n",
-    "# V.get_info(full=True)\n",
-    "# pipe.get_info()"
+    "Pegelregler.control_variable = T1.get_current_LA()\n"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 25,
+   "execution_count": 22,
    "metadata": {},
    "outputs": [],
    "source": [
     "# initialization for timeloop\n",
     "\n",
     "# prepare the vectors in which the pressure and velocity distribution in the pipeline from the previous timestep are stored\n",
-    "v_old = pipe.v0.copy()\n",
-    "p_old = pipe.p0.copy()\n",
+    "v_old = pipe.get_current_velocity_distribution()\n",
+    "p_old = pipe.get_current_pressure_distribution()\n",
     "\n",
     "# prepare the vectors in which the temporal evolution of the boundary conditions are stored\n",
-    "    # keep in mind, that the velocity at the turbine and the pressure at the reservoir are set manually and\n",
-    "        # through the time evolution of the reservoir respectively  \n",
+    "    # keep in mind, that the velocity at the turbine and the pressure at the reservoir follow from boundary conditions\n",
+    "        # reservoir level and flow through turbine\n",
     "    # the pressure at the turbine and the velocity at the reservoir are calculated from the method of characteristics\n",
-    "v_boundary_res  = np.empty_like(t_vec)\n",
-    "v_boundary_tur  = np.empty_like(t_vec)\n",
-    "p_boundary_res  = np.empty_like(t_vec)\n",
-    "p_boundary_tur  = np.empty_like(t_vec)\n",
+    "v_boundary_res  = np.zeros_like(t_vec)\n",
+    "v_boundary_tur  = np.zeros_like(t_vec)\n",
+    "p_boundary_res  = np.zeros_like(t_vec)\n",
+    "p_boundary_tur  = np.zeros_like(t_vec)\n",
     "\n",
     "# prepare the vectors that store the temporal evolution of the level in the reservoir\n",
-    "level_vec       = np.full(nt+1,V.level)  # level at the end of each pipeline timestep\n",
-    "level_vec_2     = np.empty([nt_eRK4])   # level throughout each reservoir timestep-used for plotting and overwritten afterwards\n",
+    "level_vec       = np.full(nt+1,initial_level)  # level at the end of each pipeline timestep\n",
     "\n",
     "# set the boundary conditions for the first timestep\n",
     "v_boundary_res[0]   = v_old[0]\n",
@@ -154,15 +138,14 @@
     "p_boundary_tur[0]   = p_old[-1]\n",
     "\n",
     "LA_soll_vec = np.full_like(t_vec,T1.LA)\n",
-    "Pegelregler.control_variable = T1.LA\n",
-    "\n",
+    "LA_ist_vec  = np.full_like(t_vec,T1.LA)\n",
     "\n",
     "\n"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 26,
+   "execution_count": 23,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -178,16 +161,11 @@
     "axs1[0].set_ylabel(r'$p$ ['+conversion_pressure_unit+']')\n",
     "axs1[1].set_xlabel(r'$x$ [$\\mathrm{m}$]')\n",
     "axs1[1].set_ylabel(r'$v$ [$\\mathrm{m} / \\mathrm{s}$]')\n",
-    "lo_00, = axs1[0].plot(pl_vec,pressure_conversion(pipe.p_old,initial_pressure_unit, conversion_pressure_unit)[0],marker='.')\n",
-    "lo_01, = axs1[1].plot(pl_vec,pipe.v_old,marker='.')\n",
+    "lo_00, = axs1[0].plot(pl_vec,pressure_conversion(p_old,initial_pressure_unit, conversion_pressure_unit),marker='.')\n",
+    "lo_01, = axs1[1].plot(pl_vec,v_old,marker='.')\n",
     "axs1[0].autoscale()\n",
-    "axs1[1].set_ylim([0,2])\n",
-    "# displaying the reservoir level within each pipeline timestep\n",
-    "# axs1[2].set_title('Level reservoir')\n",
-    "# axs1[2].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
-    "# axs1[2].set_ylabel(r'$h$ [m]')\n",
-    "# lo_02, = axs1[2].plot(level_vec_2)\n",
-    "# axs1[2].autoscale()\n",
+    "axs1[1].autoscale()\n",
+    "\n",
     "fig1.tight_layout()\n",
     "fig1.show()\n",
     "plt.pause(1)\n"
@@ -195,65 +173,83 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 27,
+   "execution_count": 24,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "1.2146505196687856\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(initial_flux/area_outflux)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 25,
    "metadata": {},
    "outputs": [],
    "source": [
     "# loop through time steps of the pipeline\n",
     "for it_pipe in range(1,pipe.nt+1):\n",
     "\n",
-    "    if it_pipe == 150:\n",
-    "        V.influx = 0\n",
     "\n",
+    "    if it_pipe > 0.015*(nt+1):\n",
+    "        if V.get_current_influx() > 0:\n",
+    "            V.set_influx(np.max([V.get_current_influx()-initial_flux*5*1e-3,0.]))\n",
     "# for each pipeline timestep, execute nt_eRK4 timesteps of the reservoir code\n",
     "    # set initial conditions for the reservoir time evolution calculted with e-RK4\n",
-    "    V.pressure      = p_old[0]\n",
-    "    V.outflux_vel   = v_old[0]\n",
+    "    V.set_pressure      = p_old[0]\n",
+    "    V.set_outflux       = v_old[0]\n",
     "    # calculate the time evolution of the reservoir level within each pipeline timestep to avoid runaway numerical error\n",
     "    for it_res in range(nt_eRK4):\n",
-    "        V.e_RK_4()                                                              # call e-RK4 to update outflux\n",
-    "        V.level = V.update_level(V.timestep)                                    # \n",
-    "        V.set_volume()                                                          # update volume in reservoir\n",
-    "        level_vec_2[it_res] = V.level                                           # save for plotting\n",
-    "        if (V.level < critical_level_low) or (V.level > critical_level_high):   # make sure to never exceed critical levels\n",
-    "            i_max = it_pipe                                                     # for plotting only calculated values\n",
-    "            break                                                               \n",
-    "    level_vec[it_pipe] = V.level                                                  \n",
+    "        V.timestep_reservoir_evolution()                                                             \n",
+    "    level_vec[it_pipe] = V.get_current_level()                                                 \n",
+    "\n",
+    "    # get the control variable\n",
+    "    Pegelregler.update_control_variable(level_vec[it_pipe])\n",
+    "    LA_soll_vec[it_pipe] = Pegelregler.get_current_control_variable()\n",
+    "    \n",
+    "    # change the Leitapparatöffnung  based on the target value\n",
+    "    T1.update_LA(LA_soll_vec[it_pipe])\n",
+    "    T1.set_pressure(p_old[-1])\n",
+    "\n",
+    "    LA_ist_vec[it_pipe] = T1.LA\n",
     "\n",
     "    # set boundary conditions for the next timestep of the characteristic method\n",
-    "    p_boundary_res[it_pipe] = rho*g*V.level-V.outflux_vel**2*rho/2\n",
-    "    v_boundary_res[it_pipe] = v_old[1]+1/(rho*c)*(p_boundary_res[it_pipe]-p_old[1])-f_D*dt/(2*D)*abs(v_old[1])*v_old[1] \\\n",
-    "        +dt*g*np.sin(alpha)\n",
-    "\n",
-    "    LA_soll_vec[it_pipe] = Pegelregler.get_control_variable(V.level)\n",
-    "    T1.change_LA(LA_soll_vec[it_pipe],dt)\n",
-    "    v_boundary_tur[it_pipe] = 1/A_pipe*T1.get_Q(p_old[-1])\n",
+    "    p_boundary_res[it_pipe] = V.get_current_pressure()\n",
+    "    v_boundary_tur[it_pipe] = 1/A_pipe*T1.get_current_Q()\n",
     "\n",
     "    # the the boundary conditions in the pipe.object and thereby calculate boundary pressure at turbine\n",
-    "    pipe.set_boundary_conditions_next_timestep(v_boundary_res[it_pipe],p_boundary_res[it_pipe],v_boundary_tur[it_pipe])\n",
-    "    p_boundary_tur[it_pipe] = pipe.p_boundary_tur\n",
+    "    pipe.set_boundary_conditions_next_timestep(p_boundary_res[it_pipe],v_boundary_tur[it_pipe])\n",
+    "    p_boundary_tur[it_pipe] = pipe.get_current_pressure_distribution()[-1]\n",
+    "    v_boundary_res[it_pipe] = pipe.get_current_velocity_distribution()[0]\n",
     "\n",
     "    # perform the next timestep via the characteristic method\n",
     "    pipe.timestep_characteristic_method()\n",
     "\n",
+    "    # prepare for next loop\n",
+    "    p_old = pipe.get_current_pressure_distribution()\n",
+    "    v_old = pipe.get_current_velocity_distribution()\n",
+    "\n",
     "    # plot some stuff\n",
     "        # remove line-objects to autoscale axes (there is definetly a better way, but this works ¯\\_(ツ)_/¯ )\n",
     "    lo_00.remove()\n",
     "    lo_01.remove()\n",
     "    # lo_02.remove()\n",
     "        # plot new pressure and velocity distribution in the pipeline\n",
-    "    lo_00, = axs1[0].plot(pl_vec,pressure_conversion(pipe.p_old,initial_pressure_unit, conversion_pressure_unit)[0],marker='.',c='blue')\n",
-    "    lo_01, = axs1[1].plot(pl_vec,pipe.v_old,marker='.',c='blue')\n",
+    "    lo_00, = axs1[0].plot(pl_vec,pressure_conversion(p_old,initial_pressure_unit, conversion_pressure_unit),marker='.',c='blue')\n",
+    "    lo_01, = axs1[1].plot(pl_vec,v_old,marker='.',c='blue')\n",
     "    # lo_02, = axs1[2].plot(level_vec_2,c='blue')\n",
     "    fig1.suptitle(str(round(t_vec[it_pipe],2))+ ' s / '+str(round(t_vec[-1],2)) + ' s' )\n",
     "    fig1.canvas.draw()\n",
     "    fig1.tight_layout()\n",
     "    fig1.show()\n",
-    "    plt.pause(0.1)    \n",
-    "\n",
-    "    # prepare for next loop\n",
-    "    p_old = pipe.p_old\n",
-    "    v_old = pipe.v_old  \n",
+    "    plt.pause(0.001)    \n",
     "\n",
     "        \n",
     "        "
@@ -261,37 +257,64 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 28,
+   "execution_count": 26,
    "metadata": {},
    "outputs": [],
    "source": [
     "# plot time evolution of boundary pressure and velocity as well as the reservoir level\n",
     "\n",
     "fig2,axs2 = plt.subplots(3,2)\n",
-    "axs2[0,0].plot(t_vec,pressure_conversion(p_boundary_res,initial_pressure_unit, conversion_pressure_unit)[0])\n",
-    "axs2[0,1].plot(t_vec,v_boundary_res)\n",
-    "axs2[1,0].plot(t_vec,pressure_conversion(p_boundary_tur,initial_pressure_unit, conversion_pressure_unit)[0])\n",
-    "axs2[1,1].plot(t_vec,v_boundary_tur)\n",
-    "axs2[2,0].plot(t_vec,level_vec)\n",
     "axs2[0,0].set_title('Pressure reservoir')\n",
-    "axs2[0,1].set_title('Velocity reservoir')\n",
-    "axs2[1,0].set_title('Pressure turbine')\n",
-    "axs2[1,1].set_title('Velocity turbine')\n",
-    "axs2[2,0].set_title('Level reservoir')\n",
+    "axs2[0,0].plot(t_vec,pressure_conversion(p_boundary_res,initial_pressure_unit, conversion_pressure_unit))\n",
     "axs2[0,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
     "axs2[0,0].set_ylabel(r'$p$ ['+conversion_pressure_unit+']')\n",
+    "\n",
+    "axs2[0,1].set_title('Velocity reservoir')\n",
+    "axs2[0,1].plot(t_vec,v_boundary_res)\n",
+    "axs2[0,1].set_ylim(-2*Q_nenn,+2*Q_nenn)\n",
     "axs2[0,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
     "axs2[0,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n",
+    "\n",
+    "axs2[1,0].set_title('Pressure turbine')\n",
+    "axs2[1,0].plot(t_vec,pressure_conversion(p_boundary_tur,initial_pressure_unit, conversion_pressure_unit))\n",
     "axs2[1,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
     "axs2[1,0].set_ylabel(r'$p$ ['+conversion_pressure_unit+']')\n",
+    "\n",
+    "axs2[1,1].set_title('Velocity turbine')\n",
+    "axs2[1,1].plot(t_vec,v_boundary_tur)\n",
     "axs2[1,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
     "axs2[1,1].set_ylabel(r'$v$ [$\\mathrm{m}/\\mathrm{s}$]')\n",
+    "\n",
+    "axs2[2,0].set_title('Level reservoir')\n",
+    "axs2[2,0].plot(t_vec,level_vec)\n",
     "axs2[2,0].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
     "axs2[2,0].set_ylabel(r'$h$ [m]')\n",
-    "axs2[2,1].axis('off')\n",
+    "\n",
+    "axs2[2,1].set_title('LA')\n",
+    "axs2[2,1].plot(t_vec,100*LA_soll_vec)\n",
+    "axs2[2,1].plot(t_vec,100*LA_ist_vec)\n",
+    "axs2[2,1].set_xlabel(r'$t$ [$\\mathrm{s}$]')\n",
+    "axs2[2,1].set_ylabel(r'$LA$ [%]')\n",
     "fig2.tight_layout()\n",
     "plt.show()"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 27,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "0.0\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(np.sum(v_boundary_res[2500:])*area_outflux)"
+   ]
   }
  ],
  "metadata": {