diff --git a/Ausgleichsbecken/Ausgleichsbecken_class_file.py b/Ausgleichsbecken/Ausgleichsbecken_class_file.py index 2eab201..69652dc 100644 --- a/Ausgleichsbecken/Ausgleichsbecken_class_file.py +++ b/Ausgleichsbecken/Ausgleichsbecken_class_file.py @@ -14,7 +14,7 @@ def FODE_function(x_out,h,A,A_a,p,rho,g): # based on the outflux formula by Andreas Malcherek # https://www.youtube.com/watch?v=8HO2LwqOhqQ # adapted for a pressurized pipeline into which the reservoir effuses - # and flow direction + # and flow direction # x_out ... effusion velocity # h ... level in the reservoir # A_a ... Area_outflux @@ -27,14 +27,14 @@ def FODE_function(x_out,h,A,A_a,p,rho,g): class Ausgleichsbecken_class: # units - # make sure that units and print units are the same + # make sure that units and display units are the same # units are used to label graphs and disp units are used to have a bearable format when using pythons print() area_unit = r'$\mathrm{m}^2$' area_outflux_unit = r'$\mathrm{m}^2$' density_unit = r'$\mathrm{kg}/\mathrm{m}^3$' flux_unit = r'$\mathrm{m}^3/\mathrm{s}$' level_unit = 'm' - pressure_unit = 'Pa' + pressure_unit = 'Pa' # DONT CHANGE needed for pressure conversion time_unit = 's' velocity_unit = r'$\mathrm{m}/\mathrm{s}$' volume_unit = r'$\mathrm{m}^3$' @@ -44,6 +44,7 @@ class Ausgleichsbecken_class: density_unit_disp = 'kg/m³' flux_unit_disp = 'm³/s' level_unit_disp = 'm' + # pressure_unit_disp will be set within the __init__() method time_unit_disp = 's' velocity_unit_disp = 'm/s' volume_unit_disp = 'm³' @@ -53,26 +54,37 @@ class Ausgleichsbecken_class: # init def __init__(self,area,area_outflux,timestep,pressure_unit_disp,level_min=0,level_max=np.inf,rho = 1000.): + """ + Creates a reservoir with given attributes in this order: \n + Base Area [m²] \n + Outflux Area [m²] \n + Simulation timestep [s] \n + Pressure unit for displaying [string] \n + Minimal level [m] \n + Maximal level [m] \n + Density of the liquid [kg/m³] \n + """ + #set initial attributes self.area = area # base area of the cuboid reservoir self.area_out = area_outflux # area of the outlet towards the pipeline self.density = rho # density of the liquid in the system - self.level_min = level_min # lowest allowed water level - self.level_max = level_max # highest allowed water level - self.pressure_unit_disp = pressure_unit_disp # pressure unit for displaying + self.level_min = level_min # lowest allowed water level - warning yet to be implemented + self.level_max = level_max # highest allowed water level - warning yet to be implemented + self.pressure_unit_disp = pressure_unit_disp # pressure unit for displaying self.timestep = timestep # timestep in the time evolution method - # initialize for get_info - self.influx = "--" - self.outflux = "--" - self.level = "--" - self.pressure = "--" - self.volume = "--" + # initialize for get_info() (if get_info() gets called before set_steady_state() is executed) + self.influx = -np.inf + self.outflux = -np.inf + self.level = -np.inf + self.pressure = -np.inf + self.volume = -np.inf -# setter +# setter - set attributes def set_initial_level(self,initial_level): # sets the initial level in the reservoir and should only be called during initialization - if self.level == '--': + if self.level == -np.inf: self.level = initial_level self.update_volume(set_flag=True) else: @@ -80,7 +92,7 @@ class Ausgleichsbecken_class: def set_initial_pressure(self,initial_pressure): # sets the initial static pressure present at the outlet of the reservoir and should only be called during initialization - if self.pressure == '--': + if self.pressure == -np.inf: self.pressure = initial_pressure else: raise Exception('Initial pressure was already set once. Use the .update_pressure(self) method to update pressure based current level.') @@ -96,7 +108,7 @@ class Ausgleichsbecken_class: if display_warning == True: print('You are setting the outflux from the reservoir manually. \n \ This is not an intended use of this method. \n \ - Refer to the timestep_reservoir_evolution() method instead.') + Refer to the timestep_reservoir_evolution() or set_steady_state() method instead.') self.outflux = outflux def set_level(self,level,display_warning=True): @@ -104,7 +116,7 @@ class Ausgleichsbecken_class: if display_warning == True: print('You are setting the level of the reservoir manually. \n \ This is not an intended use of this method. \n \ - Refer to the update_level() method instead.') + Refer to the update_level() or set_steady_state() method instead.') self.level = level def set_pressure(self,pressure,display_warning=True): @@ -112,48 +124,53 @@ class Ausgleichsbecken_class: if display_warning == True: print('You are setting the pressure below the reservoir manually. \n \ This is not an intended use of this method. \n \ - Refer to the update_pressure() method instead.') + Refer to the update_pressure() or set_steady_state() method instead.') self.pressure = pressure def set_volume(self,volume,display_warning=True): + # sets volume in reservoir if display_warning == True: print('You are setting the volume in the reservoir manually. \n \ This is not an intended use of this method. \n \ - Refer to the .update_volume() or set_initial_level() method instead.') + Refer to the .update_volume() or set_initial_level() or set_steady_state() method instead.') self.volume = volume def set_steady_state(self,ss_influx,ss_level): - # set the steady state (ss) condition in which the net flux is zero + # set the reservoir to 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_influx_vel = abs(ss_influx/self.area) - ss_outflux_vel = abs(ss_outflux/self.area_out) + ss_outflux = ss_influx + ss_influx_vel = abs(ss_influx/self.area) + ss_outflux_vel = abs(ss_outflux/self.area_out) + # see confluence doc for explaination on how to arrive at the ss pressure formula ss_pressure = self.density*self.g*ss_level+self.density*ss_outflux_vel*(ss_influx_vel-ss_outflux_vel) + # use setter methods to set the attributes to their steady state values self.set_influx(ss_influx) self.set_initial_level(ss_level) self.set_initial_pressure(ss_pressure) self.set_outflux(ss_outflux,display_warning=False) -# getter +# getter - return attributes def get_info(self, full = False): + # prints out the info on the current state of the reservoir new_line = '\n' - p = pressure_conversion(self.pressure,self.pressure_unit,self.pressure_unit_disp) - outflux_vel = self.outflux/self.area_out - + if self.pressure != np.inf: + p = pressure_conversion(self.pressure,self.pressure_unit,self.pressure_unit_disp) + if self.outflux != np.inf: + outflux_vel = self.outflux/self.area_out if full == True: # :<10 pads the self.value to be 10 characters wide print_str = (f"The cuboid reservoir has the following attributes: {new_line}" f"----------------------------- {new_line}" f"Base area = {self.area:<10} {self.area_unit_disp} {new_line}" - f"Outflux area = {round(self.area_out,3):<10} {self.area_out_unit_disp} {new_line}" + f"Outflux area = {round(self.area_out,3):<10} {self.area_outflux_unit_disp} {new_line}" f"Current level = {self.level:<10} {self.level_unit_disp}{new_line}" f"Critical level low = {self.level_min:<10} {self.level_unit_disp} {new_line}" f"Critical level high = {self.level_max:<10} {self.level_unit_disp} {new_line}" f"Volume in reservoir = {self.volume:<10} {self.volume_unit_disp} {new_line}" - f"Current influx = {self.influx:<10} {self.flux_unit_disp} {new_line}" - f"Current outflux = {self.outflux:<10} {self.flux_unit_disp} {new_line}" + f"Current influx = {round(self.influx,3):<10} {self.flux_unit_disp} {new_line}" + f"Current outflux = {round(self.outflux,3):<10} {self.flux_unit_disp} {new_line}" f"Current outflux vel = {round(outflux_vel,3):<10} {self.velocity_unit_disp} {new_line}" f"Current pipe pressure = {round(p,3):<10} {self.pressure_unit_disp} {new_line}" f"Simulation timestep = {self.timestep:<10} {self.time_unit_disp} {new_line}" @@ -165,8 +182,8 @@ class Ausgleichsbecken_class: f"----------------------------- {new_line}" f"Current level = {self.level:<10} {self.level_unit_disp}{new_line}" f"Current volume = {self.volume:<10} {self.volume_unit_disp} {new_line}" - f"Current influx = {self.influx:<10} {self.flux_unit_disp} {new_line}" - f"Current outflux = {self.outflux:<10} {self.flux_unit_disp} {new_line}" + f"Current influx = {round(self.influx,3):<10} {self.flux_unit_disp} {new_line}" + f"Current outflux = {round(self.outflux,3):<10} {self.flux_unit_disp} {new_line}" f"Current outflux vel = {round(outflux_vel,3):<10} {self.velocity_unit_disp} {new_line}" f"Current pipe pressure = {round(p,3):<10} {self.pressure_unit_disp} {new_line}" f"----------------------------- {new_line}") @@ -188,22 +205,26 @@ class Ausgleichsbecken_class: def get_current_volume(self): return self.volume -# update methods +# update methods - update attributes based on some parameter def update_level(self,timestep,set_flag=False): # 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 net_flux = self.influx-self.outflux delta_level = net_flux*timestep/self.area level_new = (self.level+delta_level) + # set flag is necessary because update_level() is used to get a halfstep value in the time evoultion if set_flag == True: self.set_level(level_new,display_warning=False) elif set_flag == False: return level_new def update_pressure(self,set_flag=False): - influx_vel = abs(self.influx/self.area) + # update pressure based on level and flux velocities + # see confluence doc for explaination + influx_vel = abs(self.influx/self.area) outflux_vel = abs(self.outflux/self.area_out) p_new = self.density*self.g*self.level+self.density*outflux_vel*(influx_vel-outflux_vel) + # set flag for consistency with update_level() if set_flag ==True: self.set_pressure(p_new,display_warning=False) elif set_flag == False: @@ -211,6 +232,7 @@ class Ausgleichsbecken_class: def update_volume(self,set_flag=False): volume_new = self.level*self.area + # set flag for consistency with update_level() if set_flag == True: self.set_volume(volume_new,display_warning=False) elif set_flag == False: @@ -218,7 +240,9 @@ class Ausgleichsbecken_class: #methods def timestep_reservoir_evolution(self): - # update outflux and outflux velocity based on current pipeline pressure and waterlevel in reservoir + # update outflux, level, pressure and volume based on current pipeline pressure and waterlevel in reservoir + + # get some variables dt = self.timestep rho = self.density g = self.g @@ -229,7 +253,8 @@ class Ausgleichsbecken_class: h_hs = self.update_level(dt/2) p = self.pressure p_hs = self.pressure + rho*g*(h_hs-h) - # explicit 4 step Runge Kutta + + # perform explicit 4 step Runge Kutta Y1 = yn 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) @@ -237,6 +262,7 @@ class Ausgleichsbecken_class: 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)) + # set/update the attributes to their new values self.set_outflux(ynp1*A_a,display_warning=False) self.update_level(dt,set_flag=True) self.update_volume(set_flag=True) diff --git a/Ausgleichsbecken/Ausgleichsbecken_test_steady_state.ipynb b/Ausgleichsbecken/Ausgleichsbecken_test_steady_state.ipynb index 918aa20..449ab39 100644 --- a/Ausgleichsbecken/Ausgleichsbecken_test_steady_state.ipynb +++ b/Ausgleichsbecken/Ausgleichsbecken_test_steady_state.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 29, + "execution_count": 1, "metadata": {}, "outputs": [], "source": [ @@ -21,7 +21,7 @@ }, { "cell_type": "code", - "execution_count": 30, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ @@ -33,12 +33,10 @@ "pUnit_calc = 'Pa' # [text] DO NOT CHANGE! for pressure conversion in print statements and plot labels \n", "pUnit_conv = 'mWS' # [text] for pressure conversion in print statements and plot labels\n", "\n", - "\n", " # for Turbine\n", - "Tur_Q_nenn = 0.85 # [m³/s] nominal flux of turbine \n", - "Tur_p_nenn = pressure_conversion(10.6,'bar',pUnit_calc) # [Pa] nominal pressure of turbine \n", - "Tur_closingTime = 90. # [s] closing time of turbine\n", - "\n", + "Tur_Q_nenn = 0.85 # [m³/s] nominal flux of turbine \n", + "Tur_p_nenn = pressure_conversion(10.6,'bar',pUnit_calc) # [Pa] nominal pressure of turbine \n", + "Tur_closingTime = 90. # [s] closing time of turbine\n", "\n", " # for PI controller\n", "Con_targetLevel = 8. # [m]\n", @@ -46,26 +44,24 @@ "Con_T_i = 10. # [s] timespan in which a steady state error is corrected by the intergal term\n", "Con_deadbandRange = 0.05 # [m] Deadband range around targetLevel for which the controller does NOT intervene\n", "\n", - "\n", " # for pipeline\n", - "Pip_length = (535.+478.) # [m] length of pipeline\n", - "Pip_dia = 0.9 # [m] diameter of pipeline\n", - "Pip_area = Pip_dia**2/4*np.pi # [m²] crossectional area of pipeline\n", - "Pip_head = 105. # [m] hydraulic head of pipeline without reservoir\n", - "Pip_angle = np.arcsin(Pip_head/Pip_length) # [rad] elevation angle of pipeline \n", - "Pip_n_seg = 50 # [-] number of pipe segments in discretization\n", - "Pip_f_D = 0.014 # [-] Darcy friction factor\n", - "Pip_pw_vel = 500. # [m/s] propagation velocity of the pressure wave (pw) in the given pipeline\n", + "Pip_length = (535.+478.) # [m] length of pipeline\n", + "Pip_dia = 0.9 # [m] diameter of pipeline\n", + "Pip_area = Pip_dia**2/4*np.pi # [m²] crossectional area of pipeline\n", + "Pip_head = 105. # [m] hydraulic head of pipeline without reservoir\n", + "Pip_angle = np.arcsin(Pip_head/Pip_length) # [rad] elevation angle of pipeline \n", + "Pip_n_seg = 50 # [-] number of pipe segments in discretization\n", + "Pip_f_D = 0.014 # [-] Darcy friction factor\n", + "Pip_pw_vel = 500. # [m/s] propagation velocity of the pressure wave (pw) in the given pipeline\n", " # derivatives of the pipeline constants\n", - "Pip_dx = Pip_length/Pip_n_seg # [m] length of each pipe segment\n", - "Pip_dt = Pip_dx/Pip_pw_vel # [s] timestep according to method of characteristics\n", - "Pip_nn = Pip_n_seg+1 # [1] number of nodes\n", - "Pip_x_vec = np.arange(0,Pip_nn,1)*Pip_dx # [m] vector holding the distance of each node from the upstream reservoir along the pipeline\n", - "Pip_h_vec = np.arange(0,Pip_nn,1)*Pip_head/Pip_n_seg # [m] vector holding the vertival distance of each node from the upstream reservoir\n", - "\n", + "Pip_dx = Pip_length/Pip_n_seg # [m] length of each pipe segment\n", + "Pip_dt = Pip_dx/Pip_pw_vel # [s] timestep according to method of characteristics\n", + "Pip_nn = Pip_n_seg+1 # [1] number of nodes\n", + "Pip_x_vec = np.arange(0,Pip_nn,1)*Pip_dx # [m] vector holding the distance of each node from the upstream reservoir along the pipeline\n", + "Pip_h_vec = np.arange(0,Pip_nn,1)*Pip_head/Pip_n_seg # [m] vector holding the vertival distance of each node from the upstream reservoir\n", "\n", " # for reservoir\n", - "Res_area_base = 5. # [m²] total base are of the cuboid reservoir \n", + "Res_area_base = 74. # [m²] total base are of the cuboid reservoir \n", "Res_area_out = Pip_area # [m²] outflux area of the reservoir, given by pipeline area\n", "Res_level_crit_lo = 0. # [m] for yet-to-be-implemented warnings\n", "Res_level_crit_hi = np.inf # [m] for yet-to-be-implemented warnings\n", @@ -74,25 +70,65 @@ "Res_dt = Pip_dt/Res_nt # [s] harmonised timestep of reservoir time evolution\n", "\n", " # for general simulation\n", - "flux_init = Tur_Q_nenn/1.1 # [m³/s] initial flux through whole system for steady state initialization \n", - "level_init = Con_targetLevel # [m] initial water level in upstream reservoir for steady state initialization\n", - "simTime_target = 600. # [s] target for total simulation time (will vary slightly to fit with Pip_dt)\n", - "nt = int(simTime_target//Pip_dt) # [1] Number of timesteps of the whole system\n", - "t_vec = np.arange(0,nt+1,1)*Pip_dt # [s] time vector. At each step of t_vec the system parameters are stored\n" + "flux_init = Tur_Q_nenn/1.1 # [m³/s] initial flux through whole system for steady state initialization \n", + "level_init = Con_targetLevel # [m] initial water level in upstream reservoir for steady state initialization\n", + "simTime_target = 600. # [s] target for total simulation time (will vary slightly to fit with Pip_dt)\n", + "nt = int(simTime_target//Pip_dt) # [1] Number of timesteps of the whole system\n", + "t_vec = np.arange(0,nt+1,1)*Pip_dt # [s] time vector. At each step of t_vec the system parameters are stored\n" ] }, { "cell_type": "code", - "execution_count": 31, + "execution_count": 3, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The cuboid reservoir has the following attributes: \n", + "----------------------------- \n", + "Base area = 74.0 m² \n", + "Outflux area = 0.636 m² \n", + "Current level = -inf m\n", + "Critical level low = 0.0 m \n", + "Critical level high = inf m \n", + "Volume in reservoir = -inf m³ \n", + "Current influx = -inf m³/s \n", + "Current outflux = -inf m³/s \n", + "Current outflux vel = -inf m/s \n", + "Current pipe pressure = -inf mWS \n", + "Simulation timestep = 0.001013 s \n", + "Density of liquid = 1000.0 kg/m³ \n", + "----------------------------- \n", + "\n", + "The cuboid reservoir has the following attributes: \n", + "----------------------------- \n", + "Base area = 74.0 m² \n", + "Outflux area = 0.636 m² \n", + "Current level = 8.0 m\n", + "Critical level low = 0.0 m \n", + "Critical level high = inf m \n", + "Volume in reservoir = 592.0 m³ \n", + "Current influx = 0.773 m³/s \n", + "Current outflux = 0.773 m³/s \n", + "Current outflux vel = 1.215 m/s \n", + "Current pipe pressure = 7.854 mWS \n", + "Simulation timestep = 0.001013 s \n", + "Density of liquid = 1000.0 kg/m³ \n", + "----------------------------- \n", + "\n" + ] + } + ], "source": [ "# create objects\n", "\n", "# Upstream reservoir\n", - "reservoir = Ausgleichsbecken_class(Res_area_base,Res_area_out,Res_dt,Res_level_crit_lo,Res_level_crit_hi,rho)\n", + "reservoir = Ausgleichsbecken_class(Res_area_base,Res_area_out,Res_dt,pUnit_conv,Res_level_crit_lo,Res_level_crit_hi,rho)\n", + "# print(reservoir.__init__.__doc__)\n", + "reservoir.get_info(full=True)\n", "reservoir.set_steady_state(flux_init,level_init)\n", - "\n", "reservoir.get_info(full=True)\n", "\n", "# initialize vectors\n", @@ -108,9 +144,26 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 4, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The current attributes are: \n", + "----------------------------- \n", + "Current level = 8.0 m\n", + "Current volume = 592.0 m³ \n", + "Current influx = 0.773 m³/s \n", + "Current outflux = 0.773 m³/s \n", + "Current outflux vel = 1.215 m/s \n", + "Current pipe pressure = 7.854 mWS \n", + "----------------------------- \n", + "\n" + ] + } + ], "source": [ "# time loop\n", "for i in range(1,nt+1):\n", @@ -125,16 +178,42 @@ " level_vec[i] = reservoir.get_current_level()\n", " pressure_vec[i] = reservoir.get_current_pressure()\n", "\n", - " reservoir.get_info()" + "reservoir.get_info()" ] }, { "cell_type": "code", - "execution_count": 32, + "execution_count": 5, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "3987c69b8123480db34b021e1a5393a0", + "version_major": 2, + "version_minor": 0 + }, + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoAAAAPoCAYAAABOHU+gAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAACns0lEQVR4nOzdeVhV5f7//9dmBhVwQAZFHDI1cbZMTc0yPU5pnkwtBxwqP1lpVidNcyoPjaZ1UsscU9PM4Vg5oanpsXICx3JIFNKNpCmYAwjcvz/6un/tAAXdsIX9fFzXuq72Wve613vdYL6817AtxhgjAAAAuAw3ZxcAAACAwkUABAAAcDEEQAAAABdDAAQAAHAxBEAAAAAXQwAEAABwMQRAAAAAF0MABAAAcDEEQAAAABdDAAQAAHAxBEAAAAAXQwAEAABwMQRAAAAAF0MABAAAcDEEQAAAABdDAAQAAHAxBEAAAAAXQwAEAABwMQRAAAAAF0MABFzABx98IIvFosjISGeXYrNp0yZZLBZt2rSp0PYdN26cLBaLzpw5k+9j5lbDl19+ect9FYXjAiheCICAC5g1a5Yk6cCBA/rxxx+dXM2ta9iwob7//ns1bNjQ2aUAQJFEAASKuZ07d2rPnj3q2LGjJGnmzJlOrujW+fv7695775W/v7+zSwGAIokACBRz1wLfm2++qWbNmmnRokW6dOmSbXtul1OPHz8ui8WiOXPm2NYdO3ZMPXv2VFhYmLy9vRUcHKwHH3xQcXFxtjZpaWl68cUXFRISIj8/P7Vs2VK7du1S5cqVFRUVdcN6d+7cqYcfflhlypSRj4+PGjRooC+++MKuTW41//jjj+rcubPKli0rHx8fVatWTcOGDct2jNOnT6tXr14KCAhQcHCwBgwYoJSUFLs2S5YsUZMmTRQQECA/Pz9VrVpVAwYMyNbXlStXNHz4cIWEhMjX11etWrVSbGzsTZ2XJJ08eVJPPfWUwsPD5eXlpbCwMD366KM6ffp0rmOWmpqqdu3aKTg4WNu3b5ckpaen64033lDNmjXl7e2toKAg9e/fX7/99pvdvpUrV1anTp20Zs0aNWzYUL6+vqpZs6Zt1hhA8eTh7AIAFJzLly/r888/1913363IyEgNGDBAgwYN0pIlS9SvX79899ehQwdlZmbq7bffVqVKlXTmzBlt27ZN58+ft7Xp37+/Fi9erH/961964IEHdPDgQT3yyCNKTU29Yf8bN27UP/7xDzVp0kTTp09XQECAFi1apB49eujSpUvXDZBr165V586dVatWLU2aNEmVKlXS8ePHtW7dumxt//nPf6pHjx4aOHCg9u3bp5EjR0r6/y+Vf//99+rRo4d69OihcePGycfHRydOnNC3336bra9XX31VDRs21KeffqqUlBSNGzdO999/v2JjY1W1atV8ndfJkyd199136+rVq3r11VdVt25dnT17VmvXrtW5c+cUHByc7fi//vqrOnTooPT0dH3//feqWrWqsrKy1KVLF23ZskX/+te/1KxZM504cUJjx47V/fffr507d8rX19fWx549e/Tiiy9qxIgRCg4O1qeffqqBAwfqjjvuUMuWLW/4cwNQBBkAxda8efOMJDN9+nRjjDEXLlwwJUuWNC1atLC12bhxo5FkNm7caLdvfHy8kWRmz55tjDHmzJkzRpKZPHlyrsc7cOCAkWReeeUVu/Wff/65kWT69et33ePWrFnTNGjQwFy9etVu/06dOpnQ0FCTmZmZ677VqlUz1apVM5cvX861vrFjxxpJ5u2337Zb/8wzzxgfHx+TlZVljDHm3XffNZLM+fPnc+3rWg0NGza07WeMMcePHzeenp5m0KBB+T6vAQMGGE9PT3Pw4MEbHnfJkiUmNjbWhIWFmRYtWpizZ8/a2lwb76VLl9rtu2PHDiPJTJ061bYuIiLC+Pj4mBMnTtjWXb582ZQpU8Y8/fTTudYBoGjjEjBQjM2cOVO+vr7q2bOnJKlkyZLq3r27tmzZoiNHjuSrrzJlyqhatWp65513NGnSJMXGxiorK8uuzebNmyVJjz32mN36Rx99VB4e17/gcPToUf3888964oknJEkZGRm2pUOHDrJarTp06FCO+x4+fFi//PKLBg4cKB8fnxuey8MPP2z3uW7durpy5YqSk5MlSXfffbftPL744gudPHky174ef/xxWSwW2+eIiAg1a9ZMGzduzPd5rV69Wq1bt1atWrVueA5r165VixYt1LJlS8XExKhMmTK2bV9//bUCAwPVuXNnu+PVr19fISEh2S6d169fX5UqVbJ99vHx0Z133qkTJ07csA4ARRMBECimjh49qu+++04dO3aUMUbnz5/X+fPn9eijj0pSvu/xslgs2rBhg9q1a6e3335bDRs2VFBQkJ5//nlduHBBknT27FlJynap0sPDQ2XLlr1u/9fucXvppZfk6elptzzzzDOSlOvrW67d11axYsU8ncvfa/H29pb05yVzSWrZsqVWrFihjIwM9e3bVxUrVlRkZKQ+//zzbH2FhITkuO7aWOTnvH777bc8n8OKFSt0+fJl/d///Z+t/mtOnz6t8+fPy8vLK9sxk5KSso1jTj8bb29v23gAKH64BxAopmbNmiVjjL788ssc3xk3d+5cvfHGG7YZs7S0NLvtOYWtiIgI20Mlhw8f1hdffKFx48YpPT1d06dPtwWJ06dPq0KFCrb9MjIybIEoN+XKlZMkjRw5Ut26dcuxTY0aNXJcHxQUJOnP++EcpUuXLurSpYvS0tL0ww8/KDo6Wo8//rgqV66spk2b2tolJSVl2zcpKck2Fvk5r6CgoDyfw/vvv69Fixapffv2Wr58udq2bWvbVq5cOZUtW1Zr1qzJcd9SpUrl6RgAii8CIFAMZWZmau7cuapWrZo+/fTTbNu//vprvffee1q9erUaN24sSdq7d6/atWtna7Ny5crrHuPOO+/U6NGjtXTpUu3evVuSbA8MLF682O4dfV9++aUyMjKu21+NGjVUvXp17dmzR//+97/zdqJ/qaVatWqaNWuWhg8fnm1G7FZ4e3urVatWCgwM1Nq1axUbG2sXAD///HMNHz7cdhn4xIkT2rZtm/r27Zvv82rfvr0+++wzHTp0KNewe42Pj4+WL1+u3r176+GHH9bixYvVpUsXSVKnTp20aNEiZWZmqkmTJrdy+gCKKQIgUAytXr1ap06d0ltvvaX7778/2/bIyEj95z//0cyZM9WpUye1adNG0dHRKl26tCIiIrRhwwYtW7bMbp+9e/fq2WefVffu3VW9enV5eXnp22+/1d69ezVixAhJUu3atdWrVy+99957cnd31wMPPKADBw7ovffeU0BAgNzcrn/Xyccff6z27durXbt2ioqKUoUKFfT777/rp59+0u7du7VkyZJc9/3oo4/UuXNn3XvvvXrhhRdUqVIlJSQkaO3atVqwYEG+xm/MmDH69ddf9eCDD6pixYo6f/68pkyZIk9PT7Vq1cqubXJysh555BE9+eSTSklJ0dixY+Xj42N7sjg/5zVhwgStXr1aLVu21Kuvvqo6dero/PnzWrNmjYYPH66aNWvaHdvT01Off/65Bg0apEcffVTz5s1Tr1691LNnTy1YsEAdOnTQ0KFDdc8998jT01O//vqrNm7cqC5duuiRRx7J15gAKGac/RQKAMfr2rWr8fLyMsnJybm26dmzp/Hw8DBJSUnGarWaRx991JQpU8YEBASY3r17m507d9o9BXz69GkTFRVlatasaUqUKGFKlixp6tata95//32TkZFh6/fKlStm+PDhpnz58sbHx8fce++95vvvvzcBAQHmhRdesLXL7enjPXv2mMcee8yUL1/eeHp6mpCQEPPAAw/YnmS+3r7ff/+9ad++vQkICDDe3t6mWrVqdse89hTwb7/9Zrff7NmzjSQTHx9vjDHm66+/Nu3btzcVKlQwXl5epnz58qZDhw5my5Yt2Wr47LPPzPPPP2+CgoKMt7e3adGihdm5c2e28c7LeRljTGJiohkwYIAJCQkxnp6eJiwszDz22GPm9OnTdsddsmSJbZ+srCzz/PPPGzc3NzNjxgxjjDFXr1417777rqlXr57x8fExJUuWNDVr1jRPP/20OXLkiG3fiIgI07Fjx2z1tmrVyrRq1SrbegDFg8UYY5yYPwG4gG3btql58+ZasGCBHn/8cWeXAwAujwAIwKFiYmL0/fffq1GjRvL19dWePXv05ptvKiAgQHv37s3Ta1oAAAWLewABOJS/v7/WrVunyZMn68KFCypXrpzat2+v6Ohowh8A3CaYAQQAAHAxvAgaAADAxRAAAQAAXAwBEAAAwMUQAAEAAFwMARAAAMDFEAABAABcDAEQAADAxRAAAQAAXAzfBOIAWVlZOnXqlEqVKiWLxeLscgAAQDFnjNGFCxcUFhYmN7f8z+cRAB3g1KlTCg8Pd3YZAADAxSQmJqpixYr53o8A6AClSpWS9OcPwd/f38nVAACA4i41NVXh4eG2DJJfBEAHuHbZ19/fnwAIAAAKzc3eesZDIAAAAC6GAAgAAOBiCIAAAAAuhnsAAQC4zWVmZurq1avOLgOFyNPTU+7u7gXWf7EKgBkZGRo3bpwWLFigpKQkhYaGKioqSqNHj87TO3L+97//qVWrVoqMjFRcXFzBFwwAwHUYY5SUlKTz5887uxQ4QWBgoEJCQgrkHcPFKgC+9dZbmj59uubOnavatWtr586d6t+/vwICAjR06NDr7puSkqK+ffvqwQcf1OnTpwupYgAAcnct/JUvX15+fn582YCLMMbo0qVLSk5OliSFhoY6/BjFKgB+//336tKlizp27ChJqly5sj7//HPt3Lnzhvs+/fTTevzxx+Xu7q4VK1YUcKUAAFxfZmamLfyVLVvW2eWgkPn6+kqSkpOTVb58eYdfDi5WD4Hcd9992rBhgw4fPixJ2rNnj7Zu3aoOHTpcd7/Zs2frl19+0dixYwujTAAAbujaPX9+fn5OrgTOcu1nXxD3fxarGcBXXnlFKSkpqlmzptzd3ZWZmamJEyeqV69eue5z5MgRjRgxQlu2bJGHR96GIy0tTWlpabbPqampt1w7AAA54bKv6yrIn32xmgFcvHix5s+fr4ULF2r37t2aO3eu3n33Xc2dOzfH9pmZmXr88cc1fvx43XnnnXk+TnR0tAICAmwL3wMMAEDRdvz4cVkslkJ/CNRisTjl1rNiFQBffvlljRgxQj179lSdOnXUp08fvfDCC4qOjs6x/YULF7Rz5049++yz8vDwkIeHhyZMmKA9e/bIw8ND3377bY77jRw5UikpKbYlMTGxIE8LAIAiY/r06SpVqpQyMjJs6/744w95enqqRYsWdm23bNkii8Viu3XrejZt2iSLxcIT0Q5SrC4BX7p0KdvrXtzd3ZWVlZVje39/f+3bt89u3dSpU/Xtt9/qyy+/VJUqVXLcz9vbW97e3o4pGgCAYqR169b6448/tHPnTt17772S/gx6ISEh2rFjhy5dumS7t23Tpk0KCwvL11W4W2WMUWZmZp5v+yquitUMYOfOnTVx4kR98803On78uJYvX65JkybpkUcesbUZOXKk+vbtK0lyc3NTZGSk3VK+fHn5+PgoMjJSJUqUcNapAABQJNWoUUNhYWHatGmTbd2mTZvUpUsXVatWTdu2bbNb37p1a0nS/Pnz1bhxY5UqVUohISF6/PHHba9BOX78uK1d6dKlZbFYFBUVJenPQPf222+ratWq8vX1Vb169fTll1/aHcNisWjt2rVq3LixvL29tWXLljydy8GDB9WhQweVLFlSwcHB6tOnj86cOSNJ+vjjj1WhQoVsk0wPP/yw+vXrZ/v81VdfqVGjRvLx8VHVqlU1fvx4u9lRZylWAfDDDz/Uo48+qmeeeUa1atXSSy+9pKefflqvv/66rY3ValVCQoITqwQAoHi7//77tXHjRtvnjRs36v7771erVq1s69PT0/X999/bgl16erpef/117dmzRytWrFB8fLwt5IWHh2vp0qWSpEOHDslqtWrKlCmSpNGjR2v27NmaNm2aDhw4oBdeeEG9e/fW5s2b7Wr617/+pejoaP3000+qW7fuDc/BarWqVatWql+/vnbu3Kk1a9bo9OnTeuyxxyRJ3bt315kzZ+zO89y5c1q7dq2eeOIJSdLatWvVu3dvPf/88zp48KA+/vhjzZkzRxMnTryZYXUsg1uWkpJiJJmUlBRnlwIAKCYuX75sDh48aC5fvmxbl5WVZS6mXXXKkpWVlefaP/nkE1OiRAlz9epVk5qaajw8PMzp06fNokWLTLNmzYwxxmzevNlIMr/88kuOfWzfvt1IMhcuXDDGGLNx40YjyZw7d87W5o8//jA+Pj5m27ZtdvsOHDjQ9OrVy26/FStWXLfm+Ph4I8nExsYaY4x57bXXTNu2be3aJCYmGknm0KFDxhhjHn74YTNgwADb9o8//tiEhISYjIwMY4wxLVq0MP/+97/t+vjss89MaGio7bMks3z58hxryul34JpbzR6ufQEcAIAi5PLVTN01Zq1Tjn1wQjv5eeUtNrRu3VoXL17Ujh07dO7cOd15550qX768WrVqpT59+ujixYvatGmTKlWqpKpVq0qSYmNjNW7cOMXFxen333+3XVpNSEjQXXfdlXNNBw/qypUreuihh+zWp6enq0GDBnbrGjdunK/z3bVrlzZu3KiSJUtm2/bLL7/ozjvv1BNPPKGnnnpKU6dOlbe3txYsWKCePXvaXtq8a9cu7dixw27GLzMzU1euXLG7F9IZCIAAAMCh7rjjDlWsWFEbN27UuXPn1KpVK0lSSEiIqlSpov/973/auHGjHnjgAUnSxYsX1bZtW7Vt21bz589XUFCQEhIS1K5dO6Wnp+d6nGsh8ZtvvlGFChXstv39Yc383teflZWlzp0766233sq27dpXs3Xu3FlZWVn65ptvdPfdd2vLli2aNGmSXR/jx49Xt27dsvXh4+OTr3ocjQAIAEAR4evproMT2jnt2PnRunVrbdq0SefOndPLL79sW9+qVSutXbtWP/zwg/r37y9J+vnnn3XmzBm9+eabtnfr/v1rXL28vCT9OYN2zV133SVvb28lJCTYQqajNGzYUEuXLlXlypVzfWLY19dX3bp104IFC3T06FHdeeedatSokV0fhw4d0h133OHQ2hyBAAgAQBFhsVjyfBnW2Vq3bq0hQ4bo6tWrduGsVatW+r//+z9duXLF9gBIpUqV5OXlpQ8//FCDBw/W/v377R7glKSIiAhZLBZ9/fXX6tChg3x9fVWqVCm99NJLeuGFF5SVlaX77rtPqamp2rZtm0qWLGn3NG5+DRkyRDNmzFCvXr308ssvq1y5cjp69KgWLVqkGTNm2C7zPvHEE+rcubMOHDig3r172/UxZswYderUSeHh4erevbvc3Ny0d+9e7du3T2+88cZN1+YIxeopYAAAcHto3bq1Ll++rDvuuEPBwcG29a1atdKFCxdUrVo122xfUFCQ5syZoyVLluiuu+7Sm2++qXfffdeuvwoVKmj8+PEaMWKEgoOD9eyzz0qSXn/9dY0ZM0bR0dGqVauW2rVrp6+++irXd/nmVVhYmP73v/8pMzNT7dq1U2RkpIYOHaqAgAC7dw4/8MADKlOmjA4dOqTHH3/cro927drp66+/VkxMjO6++27de++9mjRpkiIiIm6pNkew/L8nUHALUlNTFRAQoJSUFPn7+zu7HABAMXDlyhXFx8erSpUqTr9fDM5xvd+BW80ezAACAAC4GAIgAACAiyEAAgAAuBgCIAAAgIshAAIAALgYAiAAALcxXtbhugryZ08ABADgNuTp6SlJunTpkpMrgbNc+9lf+11wpKLxOnEAAFyMu7u7AgMDlZycLEny8/OTxWJxclUoDMYYXbp0ScnJyQoMDLR964gjEQABALhNhYSESJItBMK1BAYG2n4HHI0ACADAbcpisSg0NFTly5fX1atXnV0OCpGnp2eBzPxdQwAEAOA25+7uXqBhAK6Hh0AAAABcDAEQAADAxRSrAJiRkaHRo0erSpUq8vX1VdWqVTVhwgRlZWXlus/WrVvVvHlzlS1bVr6+vqpZs6bef//9QqwaAACgcBWrewDfeustTZ8+XXPnzlXt2rW1c+dO9e/fXwEBARo6dGiO+5QoUULPPvus6tatqxIlSmjr1q16+umnVaJECT311FOFfAYAAAAFz2KK0SvGO3XqpODgYM2cOdO27p///Kf8/Pz02Wef5bmfbt26qUSJEnneJzU1VQEBAUpJSZG/v3++6wYAAMiPW80exeoS8H333acNGzbo8OHDkqQ9e/Zo69at6tChQ577iI2N1bZt29SqVatc26SlpSk1NdVuAQAAKCqK1SXgV155RSkpKapZs6bc3d2VmZmpiRMnqlevXjfct2LFivrtt9+UkZGhcePGadCgQbm2jY6O1vjx4x1ZOgAAQKEpVjOAixcv1vz587Vw4ULt3r1bc+fO1bvvvqu5c+fecN8tW7Zo586dmj59uiZPnqzPP/8817YjR45USkqKbUlMTHTkaQAAABSoYnUPYHh4uEaMGKEhQ4bY1r3xxhuaP3++fv755zz388Ybb+izzz7ToUOH8tSeewABAEBh4h7Av7h06ZLc3OxPyd3d/bqvgcmJMUZpaWmOLA0AAOC2UazuAezcubMmTpyoSpUqqXbt2oqNjdWkSZM0YMAAW5uRI0fq5MmTmjdvniTpo48+UqVKlVSzZk1Jf74X8N1339Vzzz3nlHMAAAAoaMUqAH744Yd67bXX9Mwzzyg5OVlhYWF6+umnNWbMGFsbq9WqhIQE2+esrCyNHDlS8fHx8vDwULVq1fTmm2/q6aefdsYpAAAAFLhidQ+gs3APIAAAKEzcAwgAAIB8IQACAAC4GAIgAACAiyEAAgAAuBgCIAAAgIshAAIAALgYAiAAAICLIQACAAC4GAIgAACAiyEAAgAAuBgCIAAAgIshAAIAALgYAiAAAICLIQACAAC4GAIgAACAiyEAAgAAuBgCIAAAgIshAAIAALiYYhUAMzIyNHr0aFWpUkW+vr6qWrWqJkyYoKysrFz3WbZsmR566CEFBQXJ399fTZs21dq1awuxagAAgMJVrALgW2+9penTp+s///mPfvrpJ7399tt655139OGHH+a6z3fffaeHHnpIq1at0q5du9S6dWt17txZsbGxhVg5AABA4bEYY4yzi3CUTp06KTg4WDNnzrSt++c//yk/Pz999tlnee6ndu3a6tGjh8aMGZOn9qmpqQoICFBKSor8/f3zXTcAAEB+3Gr2KFYzgPfdd582bNigw4cPS5L27NmjrVu3qkOHDnnuIysrSxcuXFCZMmVybZOWlqbU1FS7BQAAoKjwcHYBjvTKK68oJSVFNWvWlLu7uzIzMzVx4kT16tUrz3289957unjxoh577LFc20RHR2v8+PGOKBkAAKDQFasZwMWLF2v+/PlauHChdu/erblz5+rdd9/V3Llz87T/559/rnHjxmnx4sUqX758ru1GjhyplJQU25KYmOioUwAAAChwxWoG8OWXX9aIESPUs2dPSVKdOnV04sQJRUdHq1+/ftfdd/HixRo4cKCWLFmiNm3aXLett7e3vL29HVY3AABAYSpWM4CXLl2Sm5v9Kbm7u1/3NTDSnzN/UVFRWrhwoTp27FiQJQIAADhdsZoB7Ny5syZOnKhKlSqpdu3aio2N1aRJkzRgwABbm5EjR+rkyZOaN2+epD/DX9++fTVlyhTde++9SkpKkiT5+voqICDAKecBAABQkIrVa2AuXLig1157TcuXL1dycrLCwsLUq1cvjRkzRl5eXpKkqKgoHT9+XJs2bZIk3X///dq8eXO2vvr166c5c+bk6bi8BgYAABSmW80exSoAOgsBEAAAFCbeAwgAAIB8IQACAAC4GAIgAACAiyEAAgAAuBgCIAAAgIshAAIAALgYAiAAAICLIQACAAC4GAIgAACAiyEAAgAAuBgCIAAAgIshAAIAALgYAiAAAICLIQACAAC4GAIgAACAiyEAAgAAuBgCIAAAgIshAAIAALiYYhUAMzIyNHr0aFWpUkW+vr6qWrWqJkyYoKysrFz3sVqtevzxx1WjRg25ublp2LBhhVcwAACAE3g4uwBHeuuttzR9+nTNnTtXtWvX1s6dO9W/f38FBARo6NChOe6TlpamoKAgjRo1Su+//34hVwwAAFD4ilUA/P7779WlSxd17NhRklS5cmV9/vnn2rlzZ677VK5cWVOmTJEkzZo1q1DqBAAAcKZidQn4vvvu04YNG3T48GFJ0p49e7R161Z16NDByZUBAADcPorVDOArr7yilJQU1axZU+7u7srMzNTEiRPVq1cvhx4nLS1NaWlpts+pqakO7R8AAKAgFasZwMWLF2v+/PlauHChdu/erblz5+rdd9/V3LlzHXqc6OhoBQQE2Jbw8HCH9g8AAFCQLMYY4+wiHCU8PFwjRozQkCFDbOveeOMNzZ8/Xz///PMN97///vtVv359TZ48+brtcpoBDA8PV0pKivz9/W+6fgAAgLxITU1VQEDATWePYnUJ+NKlS3Jzs5/UdHd3v+5rYG6Gt7e3vL29HdonAABAYSlWAbBz586aOHGiKlWqpNq1ays2NlaTJk3SgAEDbG1GjhypkydPat68ebZ1cXFxkqQ//vhDv/32m+Li4uTl5aW77rqrsE8BAACgwBWrS8AXLlzQa6+9puXLlys5OVlhYWHq1auXxowZIy8vL0lSVFSUjh8/rk2bNtn2s1gs2fqKiIjQ8ePH83TcW52GBQAAyI9bzR7FKgA6CwEQAAAUplvNHsXqKWAAAADcGAEQAADAxRAAAQAAXAwBEAAAwMUQAAEAAFxMob8HcOXKlfne56GHHpKvr28BVAMAAOB6Cj0Adu3aNV/tLRaLjhw5oqpVqxZMQQAAAC7GKZeAk5KSlJWVlafFz8/PGSUCAAAUW4UeAPv165evy7m9e/fm5coAAAAOxDeBOADfBAIAAAoT3wQCAACAfCn0h0D+7sqVK9q7d6+Sk5OVlZVlt+3hhx92UlUAAADFl1MD4Jo1a9S3b1+dOXMm2zaLxaLMzEwnVAUAAFC8OfUS8LPPPqvu3bvLarVme/qX8AcAAFAwnBoAk5OTNXz4cAUHBzuzDAAAAJfi1AD46KOPatOmTc4sAQAAwOU49TUwly5dUvfu3RUUFKQ6derI09PTbvvzzz/vpMryh9fAAACAwnSr2cOpD4EsXLhQa9eula+vrzZt2iSLxWLbZrFYikwABAAAKEqcegl49OjRmjBhglJSUnT8+HHFx8fblmPHjuW7v4yMDI0ePVpVqlSRr6+vqlatqgkTJmR7vczfbd68WY0aNZKPj4+qVq2q6dOn3+wpAQAA3PacOgOYnp6uHj16yM3NMTn0rbfe0vTp0zV37lzVrl1bO3fuVP/+/RUQEKChQ4fmuE98fLw6dOigJ598UvPnz9f//vc/PfPMMwoKCtI///lPh9QFAABwO3HqPYAvvPCCgoKC9Oqrrzqkv06dOik4OFgzZ860rfvnP/8pPz8/ffbZZznu88orr2jlypX66aefbOsGDx6sPXv26Pvvv8/TcbkHEAAAFKYifQ9gZmam3n77ba1du1Z169bN9hDIpEmT8tXffffdp+nTp+vw4cO68847tWfPHm3dulWTJ0/OdZ/vv/9ebdu2tVvXrl07zZw5U1evXs1WkzMkpVzRH2kZzi4DAADkQ7C/t0r5OD9H5MSpAXDfvn1q0KCBJGn//v122/76QEhevfLKK0pJSVHNmjXl7u6uzMxMTZw4Ub169cp1n6SkpGzvIQwODlZGRobOnDmj0NDQbPukpaUpLS3N9jk1NTXftebHhK8PaNW+pAI9BgAAcKwPejXQw/XCnF1GjpwaADdu3OjQ/hYvXqz58+dr4cKFql27tuLi4jRs2DCFhYWpX79+ue7397B57ap4biE0Ojpa48ePd1zhN+Dn5aEA39vzXxAAACBnXu75n8wqLIV+D+DevXsVGRmZ5wc/Dhw4oBo1asjD48ZZNTw8XCNGjNCQIUNs69544w3Nnz9fP//8c477tGzZUg0aNNCUKVNs65YvX67HHntMly5dyvEScE4zgOHh4dwDCAAACsWt3gNY6K+BadCggc6ePZvn9k2bNlVCQkKe2l66dClbsHR3d7/ua2CaNm2qmJgYu3Xr1q1T48aNc73/z9vbW/7+/nYLAABAUVHol4CNMXrttdfk5+eXp/bp6el57rtz586aOHGiKlWqpNq1ays2NlaTJk3SgAEDbG1GjhypkydPat68eZL+fOL3P//5j4YPH64nn3xS33//vWbOnKnPP/88fycGAABQRBR6AGzZsqUOHTqU5/ZNmzaVr69vntp++OGHeu211/TMM88oOTlZYWFhevrppzVmzBhbG6vVajejWKVKFa1atUovvPCCPvroI4WFhemDDz7gHYAAAKDYcup7AIsL3gMIAAAKU5G7BxAAAADORQAEAABwMQRAAAAAF0MABAAAcDEEQAAAABfj1K+C27Fjh0aMGKHffvtNd9xxh+rXr29bKlWq5MzSAAAAii2nzgD26dNH7u7uGjx4sKpWrarNmzerf//+qly5ssqWLevM0gAAAIotp84AJiYm6ptvvlG1atXs1p84cUJxcXHOKQoAAKCYc2oAbN68uRITE7MFwIiICEVERDipKgAAgOKt0ANgly5dVK9ePdWrV0+DBw/WhAkTVKdOHS75AgAAFJJCD4DVq1fXtm3bNG3aNJ09e1aSVKNGDXXp0kVNmzZVgwYNVKdOHXl5eRV2aQAAAC7Bqd8F/OuvvyouLs5uiY+Pl7u7u2rWrKm9e/c6q7R84buAAQBAYbrV7OHUewArVqyoihUrqlOnTrZ1f/zxh2JjY4tM+AMAAChqnDoDWFwwAwgAAArTrWYPvgkEAADAxRAAAQAAXAwBEAAAwMUQAAEAAFwMARAAAMDFFKsAWLlyZVkslmzLkCFDct3no48+Uq1ateTr66saNWpo3rx5hVgxAABA4XPqewAdbceOHcrMzLR93r9/vx566CF17949x/bTpk3TyJEjNWPGDN19993avn27nnzySZUuXVqdO3curLIBAAAKVbF+D+CwYcP09ddf68iRI7JYLNm2N2vWTM2bN9c777xjt8/OnTu1devWPB+H9wACAIDCVKS/CaQgpaena/78+Ro+fHiO4U+S0tLS5OPjY7fO19dX27dv19WrV+Xp6ZnrfmlpabbPqampjiscAACggBWrewD/asWKFTp//ryioqJybdOuXTt9+umn2rVrl4wx2rlzp2bNmqWrV6/qzJkzue4XHR2tgIAA2xIeHl4AZwAAAFAwiu0l4Hbt2snLy0tfffVVrm0uX76sIUOG6LPPPpMxRsHBwerdu7fefvttnT59WuXLl89xv5xmAMPDw7kEDAAACgVfBZeDEydOaP369Ro0aNB12/n6+mrWrFm6dOmSjh8/roSEBFWuXFmlSpVSuXLlct3P29tb/v7+dgsAAEBRUSzvAZw9e7bKly+vjh075qm9p6enKlasKElatGiROnXqJDe3YpmNAQAAil8AzMrK0uzZs9WvXz95eNif3siRI3Xy5Enbu/4OHz6s7du3q0mTJjp37pwmTZqk/fv3a+7cuc4oHQAAoFAUuwC4fv16JSQkaMCAAdm2Wa1WJSQk2D5nZmbqvffe06FDh+Tp6anWrVtr27Ztqly5ciFWDAAAULiK7UMghYn3AAIAgMLEQyAAAADIFwIgAACAiyEAAgAAuBgCIAAAgIshAAIAALgYAiAAAICLIQACAAC4GAIgAACAiyEAAgAAuBgCIAAAgIshAAIAALgYAiAAAICLIQACAAC4GAIgAACAiyEAAgAAuBgCIAAAgIshAAIAALgYAiAAAICLKVYBsHLlyrJYLNmWIUOG5LrPggULVK9ePfn5+Sk0NFT9+/fX2bNnC7FqAACAwlWsAuCOHTtktVptS0xMjCSpe/fuObbfunWr+vbtq4EDB+rAgQNasmSJduzYoUGDBhVm2QAAAIXKw9kFOFJQUJDd5zfffFPVqlVTq1atcmz/ww8/qHLlynr++eclSVWqVNHTTz+tt99+u8BrBQAAcJZiNQP4V+np6Zo/f74GDBggi8WSY5tmzZrp119/1apVq2SM0enTp/Xll1+qY8eO1+07LS1NqampdgsAAEBRUWwD4IoVK3T+/HlFRUXl2qZZs2ZasGCBevToIS8vL4WEhCgwMFAffvjhdfuOjo5WQECAbQkPD3dw9QAAAAXHYowxzi6iILRr105eXl766quvcm1z8OBBtWnTRi+88ILatWsnq9Wql19+WXfffbdmzpyZ635paWlKS0uzfU5NTVV4eLhSUlLk7+/v0PMAAAD4u9TUVAUEBNx09iiWAfDEiROqWrWqli1bpi5duuTark+fPrpy5YqWLFliW7d161a1aNFCp06dUmhoaJ6Od6s/BAAAgPy41exRLC8Bz549W+XLl7/hvXyXLl2Sm5v9ELi7u0uSimEuBgAAkFQMA2BWVpZmz56tfv36ycPD/iHnkSNHqm/fvrbPnTt31rJlyzRt2jQdO3ZM//vf//T888/rnnvuUVhYWGGXDgAAUCiK1WtgJGn9+vVKSEjQgAEDsm2zWq1KSEiwfY6KitKFCxf0n//8Ry+++KICAwP1wAMP6K233irMkgEAAApVsbwHsLBxDyAAAChM3AMIAACAfCEAAgAAuBgCIAAAgIshAAIAALgYAiAAAICLIQACAAC4GAIgAACAiyEAAgAAuBgCIAAAgIshAAIAALgYAiAAAICLIQACAAC4GAIgAACAiyEAAgAAuBgCIAAAgIshAAIAALgYAiAAAICLIQACAAC4mGIVACtXriyLxZJtGTJkSI7to6Kicmxfu3btQq4cAACg8BSrALhjxw5ZrVbbEhMTI0nq3r17ju2nTJli1z4xMVFlypTJtT0AAEBx4OHsAhwpKCjI7vObb76patWqqVWrVjm2DwgIUEBAgO3zihUrdO7cOfXv379A6wQAAHCmYjUD+Ffp6emaP3++BgwYIIvFkqd9Zs6cqTZt2igiIqKAqwMAAHCeYjUD+FcrVqzQ+fPnFRUVlaf2VqtVq1ev1sKFC2/YNi0tTWlpabbPqampN1smAABAoSu2M4AzZ85U+/btFRYWlqf2c+bMUWBgoLp27XrDttHR0bbLxwEBAQoPD7/FagEAAApPsQyAJ06c0Pr16zVo0KA8tTfGaNasWerTp4+8vLxu2H7kyJFKSUmxLYmJibdaMgAAQKEplpeAZ8+erfLly6tjx455ar9582YdPXpUAwcOzFN7b29veXt730qJAAAATlPsZgCzsrI0e/Zs9evXTx4e9vl25MiR6tu3b7Z9Zs6cqSZNmigyMrKwygQAAHCaYhcA169fr4SEBA0YMCDbNqvVqoSEBLt1KSkpWrp0aZ5n/wAAAIo6izHGOLuIoi41NVUBAQFKSUmRv7+/s8sBAADF3K1mj2I3AwgAAIDrIwACAAC4GAIgAACAiyEAAgAAuBgCIAAAgIshAAIAALgYAiAAAICLIQACAAC4GAIgAACAiyEAAgAAuBgCIAAAgIshAAIAALgYAiAAAICLIQACAAC4GAIgAACAiyEAAgAAuBgCIAAAgIshAAIAALgYAiAAAICLKVYBsHLlyrJYLNmWIUOG5LpPWlqaRo0apYiICHl7e6tatWqaNWtWIVYNAABQuDycXYAj7dixQ5mZmbbP+/fv10MPPaTu3bvnus9jjz2m06dPa+bMmbrjjjuUnJysjIyMwigXAADAKYpVAAwKCrL7/Oabb6patWpq1apVju3XrFmjzZs369ixYypTpoykP2cRAQAAirNidQn4r9LT0zV//nwNGDBAFoslxzYrV65U48aN9fbbb6tChQq688479dJLL+ny5cvX7TstLU2pqal2CwAAQFFRrGYA/2rFihU6f/68oqKicm1z7Ngxbd26VT4+Plq+fLnOnDmjZ555Rr///vt17wOMjo7W+PHjC6BqAACAgmcxxhhnF1EQ2rVrJy8vL3311Ve5tmnbtq22bNmipKQkBQQESJKWLVumRx99VBcvXpSvr2+O+6WlpSktLc32OTU1VeHh4UpJSZG/v79jTwQAAOBvUlNTFRAQcNPZo1jOAJ44cULr16/XsmXLrtsuNDRUFSpUsIU/SapVq5aMMfr1119VvXr1HPfz9vaWt7e3Q2sGAAAoLMXyHsDZs2erfPny6tix43XbNW/eXKdOndIff/xhW3f48GG5ubmpYsWKBV0mAACAUxS7AJiVlaXZs2erX79+8vCwn+AcOXKk+vbta/v8+OOPq2zZsurfv78OHjyo7777Ti+//LIGDBiQ6+VfAACAoq7YBcD169crISFBAwYMyLbNarUqISHB9rlkyZKKiYnR+fPn1bhxYz3xxBPq3LmzPvjgg8IsGQAAoFAV24dACtOt3ogJAACQH7eaPYrdDCAAAACujwAIAADgYgiAAAAALoYACAAA4GKK5YugC9u152j4TmAAAFAYrmWOm32WlwDoABcuXJAkhYeHO7kSAADgSi5cuGD3jWZ5xWtgHCArK0unTp1SqVKlZLFYHN7/te8aTkxM5DUzt4ixdAzG0TEYR8dhLB2DcXSMwhhHY4wuXLigsLAwubnl/44+ZgAdoLC+Os7f358/kA7CWDoG4+gYjKPjMJaOwTg6RkGP483M/F3DQyAAAAAuhgAIAADgYgiARYC3t7fGjh0rb29vZ5dS5DGWjsE4Ogbj6DiMpWMwjo5RFMaRh0AAAABcDDOAAAAALoYACAAA4GIIgAAAAC6GAFgETJ06VVWqVJGPj48aNWqkLVu2OLuk28p3332nzp07KywsTBaLRStWrLDbbozRuHHjFBYWJl9fX91///06cOCAXZu0tDQ999xzKleunEqUKKGHH35Yv/76ayGehfNFR0fr7rvvVqlSpVS+fHl17dpVhw4dsmvDWN7YtGnTVLduXdv7v5o2barVq1fbtjOGNyc6OloWi0XDhg2zrWMs82bcuHGyWCx2S0hIiG0745h3J0+eVO/evVW2bFn5+fmpfv362rVrl217kRpLg9vaokWLjKenp5kxY4Y5ePCgGTp0qClRooQ5ceKEs0u7baxatcqMGjXKLF261Egyy5cvt9v+5ptvmlKlSpmlS5eaffv2mR49epjQ0FCTmppqazN48GBToUIFExMTY3bv3m1at25t6tWrZzIyMgr5bJynXbt2Zvbs2Wb//v0mLi7OdOzY0VSqVMn88ccftjaM5Y2tXLnSfPPNN+bQoUPm0KFD5tVXXzWenp5m//79xhjG8GZs377dVK5c2dStW9cMHTrUtp6xzJuxY8ea2rVrG6vValuSk5Nt2xnHvPn9999NRESEiYqKMj/++KOJj48369evN0ePHrW1KUpjSQC8zd1zzz1m8ODBdutq1qxpRowY4aSKbm9/D4BZWVkmJCTEvPnmm7Z1V65cMQEBAWb69OnGGGPOnz9vPD09zaJFi2xtTp48adzc3MyaNWsKrfbbTXJyspFkNm/ebIxhLG9F6dKlzaeffsoY3oQLFy6Y6tWrm5iYGNOqVStbAGQs827s2LGmXr16OW5jHPPulVdeMffdd1+u24vaWHIJ+DaWnp6uXbt2qW3btnbr27Ztq23btjmpqqIlPj5eSUlJdmPo7e2tVq1a2cZw165dunr1ql2bsLAwRUZGuvQ4p6SkSJLKlCkjibG8GZmZmVq0aJEuXryopk2bMoY3YciQIerYsaPatGljt56xzJ8jR44oLCxMVapUUc+ePXXs2DFJjGN+rFy5Uo0bN1b37t1Vvnx5NWjQQDNmzLBtL2pjSQC8jZ05c0aZmZkKDg62Wx8cHKykpCQnVVW0XBun641hUlKSvLy8VLp06VzbuBpjjIYPH6777rtPkZGRkhjL/Ni3b59Kliwpb29vDR48WMuXL9ddd93FGObTokWLtGvXLkVHR2fbxljmXZMmTTRv3jytXbtWM2bMUFJSkpo1a6azZ88yjvlw7NgxTZs2TdWrV9fatWs1ePBgPf/885o3b56kovc76VGoR8NNsVgsdp+NMdnW4fpuZgxdeZyfffZZ7d27V1u3bs22jbG8sRo1aiguLk7nz5/X0qVL1a9fP23evNm2nTG8scTERA0dOlTr1q2Tj49Pru0Yyxtr37697b/r1Kmjpk2bqlq1apo7d67uvfdeSYxjXmRlZalx48b697//LUlq0KCBDhw4oGnTpqlv3762dkVlLJkBvI2VK1dO7u7u2f5VkJycnO1fGMjZtSfdrjeGISEhSk9P17lz53Jt40qee+45rVy5Uhs3blTFihVt6xnLvPPy8tIdd9yhxo0bKzo6WvXq1dOUKVMYw3zYtWuXkpOT1ahRI3l4eMjDw0ObN2/WBx98IA8PD9tYMJb5V6JECdWpU0dHjhzhdzIfQkNDddddd9mtq1WrlhISEiQVvf9HEgBvY15eXmrUqJFiYmLs1sfExKhZs2ZOqqpoqVKlikJCQuzGMD09XZs3b7aNYaNGjeTp6WnXxmq1av/+/S41zsYYPfvss1q2bJm+/fZbValSxW47Y3nzjDFKS0tjDPPhwQcf1L59+xQXF2dbGjdurCeeeEJxcXGqWrUqY3mT0tLS9NNPPyk0NJTfyXxo3rx5tldjHT58WBEREZKK4P8jC/WRE+TbtdfAzJw50xw8eNAMGzbMlChRwhw/ftzZpd02Lly4YGJjY01sbKyRZCZNmmRiY2Ntr8p58803TUBAgFm2bJnZt2+f6dWrV46P5VesWNGsX7/e7N692zzwwAMu94qD//u//zMBAQFm06ZNdq+LuHTpkq0NY3ljI0eONN99952Jj483e/fuNa+++qpxc3Mz69atM8Ywhrfir08BG8NY5tWLL75oNm3aZI4dO2Z++OEH06lTJ1OqVCnb3yOMY95s377deHh4mIkTJ5ojR46YBQsWGD8/PzN//nxbm6I0lgTAIuCjjz4yERERxsvLyzRs2ND2Wg78aePGjUZStqVfv37GmD8fzR87dqwJCQkx3t7epmXLlmbfvn12fVy+fNk8++yzpkyZMsbX19d06tTJJCQkOOFsnCenMZRkZs+ebWvDWN7YgAEDbH9eg4KCzIMPPmgLf8Ywhrfi7wGQscyba++i8/T0NGFhYaZbt27mwIEDtu2MY9599dVXJjIy0nh7e5uaNWuaTz75xG57URpLizHGFO6cIwAAAJyJewABAABcDAEQAADAxRAAAQAAXAwBEAAAwMUQAAEAAFwMARAAAMDFEAABAABcDAEQAADAxRAAAQAAXAwBEAAK2f333y+LxSKLxaK4uLg87RMVFWXbZ8WKFQVaH4DijwAIAA42bNgwde3a9bptnnzySVmtVkVGRuapzylTpshqtTqgOgAgAAKAw+3YsUP33HPPddv4+fkpJCREHh4eeeozICBAISEhjigPAAiAAOAoV69elZeXl7Zt26ZRo0bJYrGoSZMmed7/yy+/VJ06deTr66uyZcuqTZs2unjxYgFWDMBV5e2fngCAG3J3d9fWrVvVpEkTxcXFKTg4WD4+Pnna12q1qlevXnr77bf1yCOP6MKFC9qyZYuMMQVcNQBXRAAEAAdxc3PTqVOnVLZsWdWrVy9f+1qtVmVkZKhbt26KiIiQJNWpU6cgygQALgEDgCPFxsbmO/xJUr169fTggw+qTp066t69u2bMmKFz584VQIUAQAAEAIeKi4u7qQDo7u6umJgYrV69WnfddZc+/PBD1ahRQ/Hx8QVQJQBXRwAEAAfat2+f6tate1P7WiwWNW/eXOPHj1dsbKy8vLy0fPlyB1cIANwDCAAOlZWVpb179+rUqVMqUaKEAgIC8rTfjz/+qA0bNqht27YqX768fvzxR/3222+qVatWAVcMwBUxAwgADvTGG29o8eLFqlChgiZMmJDn/fz9/fXdd9+pQ4cOuvPOOzV69Gi99957at++fQFWC8BVMQMIAA7Uu3dv9e7dO9/71apVS2vWrCmAigAgO2YAAcAJpk6dqpIlS2rfvn15aj948GCVLFmygKsC4CoshreMAkChOnnypC5fvixJqlSpkry8vG64T3JyslJTUyVJoaGhKlGiRIHWCKB4IwACAAC4GC4BAwAAuBgCIAAAgIshAAIAALgYAiAAAICLIQACAAC4GAIgAACAiyEAAgAAuBgCIAAAgIshAAIAALgYAiAAAICLIQACAAC4GAIgAACAiyEAAgAAuBgCIAAAgIvxcHYBxUFWVpZOnTqlUqVKyWKxOLscAABQzBljdOHCBYWFhcnNLf/zeQRABzh16pTCw8OdXQYAAHAxiYmJqlixYr73IwA6QKlSpST9+UPw9/d3cjUAAKC4S01NVXh4uC2D5BcB0AGuXfb19/cnAAIAgEJzs7ee8RAIAACAiyEAAgAAuBgCIAAAgIvhHkAAAGAnMzNTV69edXYZLs3T01Pu7u4F1j8BEAAASPrz3XJJSUk6f/68s0uBpMDAQIWEhBTIO4YJgAAAQJJs4a98+fLy8/Pjyw2cxBijS5cuKTk5WZIUGhrq8GMQAAEAgDIzM23hr2zZss4ux+X5+vpKkpKTk1W+fHmHXw7mIRAAAGC758/Pz8/JleCaaz+LgrgfkwAIAABsuOx7+yjIn0WRC4BTp05VlSpV5OPjo0aNGmnLli3Xbb9gwQLVq1dPfn5+Cg0NVf/+/XX27Fnb9gMHDuif//ynKleuLIvFosmTJxfwGQAAADhXkQqAixcv1rBhwzRq1CjFxsaqRYsWat++vRISEnJsv3XrVvXt21cDBw7UgQMHtGTJEu3YsUODBg2ytbl06ZKqVq2qN998UyEhIYV1KgAAoAgZN26cgoODZbFYtGLFCkVFRalr167OLuumFakAOGnSJA0cOFCDBg1SrVq1NHnyZIWHh2vatGk5tv/hhx9UuXJlPf/886pSpYruu+8+Pf3009q5c6etzd1336133nlHPXv2lLe3d2GdCgAAcKDExEQNHDhQYWFh8vLyUkREhIYOHWp31e9Gjh8/LovFori4OLv1P/30k8aPH6+PP/5YVqtV7du3d3D1ha/IBMD09HTt2rVLbdu2tVvftm1bbdu2Lcd9mjVrpl9//VWrVq2SMUanT5/Wl19+qY4dO95SLWlpaUpNTbVbAACAcxw7dkyNGzfW4cOH9fnnn+vo0aOaPn26NmzYoKZNm+r333+/pf5/+eUXSVKXLl0UEhJSLCaMikwAPHPmjDIzMxUcHGy3Pjg4WElJSTnu06xZMy1YsEA9evSQl5eXQkJCFBgYqA8//PCWaomOjlZAQIBtCQ8Pv6X+AADAzRsyZIi8vLy0bt06tWrVSpUqVVL79u21fv16nTx5UqNGjZIk2+XbvwoMDNScOXMkSVWqVJEkNWjQQBaLRffff7/GjRunzp07S5Lc3NxyfTCjcuXK2Z4jqF+/vsaNGydJ2rRpk7y8vOyeXXjvvfdUrlw5Wa3WWxyB/CsyAfCavw+8MSbXH8bBgwf1/PPPa8yYMdq1a5fWrFmj+Ph4DR48+JZqGDlypFJSUmxLYmLiLfUHAMDtyBijS+kZhb4YY/Jc4++//661a9fqmWeesb0775qQkBA98cQTWrx4cZ763L59uyRp/fr1slqtWrZsmV566SXNnj1bkmS1Wm86rN1///0aNmyY+vTpo5SUFO3Zs0ejRo3SjBkzCuRFzzdSZF4EXa5cObm7u2eb7UtOTs42K3hNdHS0mjdvrpdfflmSVLduXZUoUUItWrTQG2+8cdMD7u3tXSymfwEAuJ7LVzN115i1hX7cgxPayc8rbxHlyJEjMsaoVq1aOW6vVauWzp07p99+++2GfQUFBUmSypYta/dgaGBgoCTd8sOib7zxhtavX6+nnnpKBw4cUJ8+ffTII4/cUp83q8jMAHp5ealRo0aKiYmxWx8TE6NmzZrluM+lS5fk5mZ/itfepJ2ff10AAICi6drf97fD+w29vLw0f/58LV26VJcvX3bqq+eKzAygJA0fPlx9+vRR48aN1bRpU33yySdKSEiwXdIdOXKkTp48qXnz5kmSOnfurCeffFLTpk1Tu3btZLVaNWzYMN1zzz0KCwuT9OfDJQcPHrT998mTJxUXF6eSJUvqjjvucM6JAgBwG/D1dNfBCe2ccty8uuOOO2SxWHTw4MEcX8vy888/q3Tp0ipXrpwsFku2CSBHfcuGm5tbnvq+9uDq77//rt9//10lSpRwyPHzq0gFwB49eujs2bOaMGGCrFarIiMjtWrVKkVEREj689r8X98JGBUVpQsXLug///mPXnzxRQUGBuqBBx7QW2+9ZWtz6tQpNWjQwPb53Xff1bvvvqtWrVpp06ZNhXZuAADcbiwWS54vxTpL2bJl9dBDD2nq1Kl64YUX7O4DTEpK0oIFC9S3b19ZLBYFBQXZ3cN35MgRXbp0yfbZy8tL0p/fi5xff+87NTVV8fHxdm1++eUXvfDCC5oxY4a++OIL9e3bVxs2bMh2tbIw3N4/1Rw888wzeuaZZ3Lcdu0pnr967rnn9Nxzz+XaX+XKlbkcDABAEfaf//xHzZo1U7t27fTGG2+oSpUqOnDggF5++WVVqFBBEydOlCQ98MAD+s9//qN7771XWVlZeuWVV+Tp6Wnrp3z58vL19dWaNWtUsWJF+fj4KCAgIE81PPDAA5ozZ446d+6s0qVL67XXXrPddib9GSr79Omjtm3bqn///mrfvr3q1Kmj9957z/asQmEqMvcAAgAA5KR69erauXOnqlWrph49eqhatWp66qmn1Lp1a33//fcqU6aMpD9fuxIeHq6WLVvq8ccf10svvSQ/Pz9bPx4eHvrggw/08ccfKywsTF26dMlzDSNHjlTLli3VqVMndejQQV27dlW1atVs2ydOnKjjx4/rk08+kfTnAyWffvqpRo8ene3F04XBYpj+umWpqakKCAhQSkqK/P39nV0OAAD5duXKFcXHx6tKlSry8fFxdjnQ9X8mt5o9mAEEAABwMQRAAAAAF0MABAAAcDEEQAAAABdDAAQAAHAxBEAAAGCTlZXl7BLw/xTkz6LIvQgaAAA4npeXl9zc3HTq1CkFBQXJy8vrtvj+XFdkjFF6erp+++03ubm52b6hxJEIgAAAQG5ubqpSpYqsVqtOnTrl7HIgyc/PT5UqVSqQr4ojAAIAAEl/zgJWqlRJGRkZN/V9uHAcd3d3eXh4FNgsLAEQAADYWCwWeXp62n1HLoofHgIBAABwMQRAAAAAF0MABAAAcDEEQAAAABdDAAQAAHAxBEAAAAAXQwAEAABwMUUuAE6dOlVVqlSRj4+PGjVqpC1btly3/YIFC1SvXj35+fkpNDRU/fv319mzZ+3aLF26VHfddZe8vb111113afny5QV5CgAAAE5VpALg4sWLNWzYMI0aNUqxsbFq0aKF2rdvr4SEhBzbb926VX379tXAgQN14MABLVmyRDt27NCgQYNsbb7//nv16NFDffr00Z49e9SnTx899thj+vHHHwvrtAAAAAqVxRhjnF1EXjVp0kQNGzbUtGnTbOtq1aqlrl27Kjo6Olv7d999V9OmTdMvv/xiW/fhhx/q7bffVmJioiSpR48eSk1N1erVq21t/vGPf6h06dL6/PPP81RXamqqAgIClJKSIn9//5s9PQAAgDy51exRZGYA09PTtWvXLrVt29Zufdu2bbVt27Yc92nWrJl+/fVXrVq1SsYYnT59Wl9++aU6duxoa/P9999n67Ndu3a59gkAAFDUFZkAeObMGWVmZio4ONhufXBwsJKSknLcp1mzZlqwYIF69OghLy8vhYSEKDAwUB9++KGtTVJSUr76lKS0tDSlpqbaLQAAAEVFkQmA11gsFrvPxphs6645ePCgnn/+eY0ZM0a7du3SmjVrFB8fr8GDB990n5IUHR2tgIAA2xIeHn6TZwMAAFD4ikwALFeunNzd3bPNzCUnJ2ebwbsmOjpazZs318svv6y6deuqXbt2mjp1qmbNmiWr1SpJCgkJyVefkjRy5EilpKTYlmv3EwIAABQFRSYAenl5qVGjRoqJibFbHxMTo2bNmuW4z6VLl+TmZn+K7u7ukv6c5ZOkpk2bZutz3bp1ufYpSd7e3vL397dbAAAAigoPZxeQH8OHD1efPn3UuHFjNW3aVJ988okSEhJsl3RHjhypkydPat68eZKkzp0768knn9S0adPUrl07Wa1WDRs2TPfcc4/CwsIkSUOHDlXLli311ltvqUuXLvrvf/+r9evXa+vWrU47TwAAgIJUpAJgjx49dPbsWU2YMEFWq1WRkZFatWqVIiIiJElWq9XunYBRUVG6cOGC/vOf/+jFF19UYGCgHnjgAb311lu2Ns2aNdOiRYs0evRovfbaa6pWrZoWL16sJk2aFPr5AQAAFIYi9R7A2xXvAQQAAIXJZd4DCAAAAMcgAAIAALgYAiAAAICLIQACAAC4GAIgAACAiyEAAgAAuBgCIAAAgIshAAIAALgYAiAAAICLIQACAAC4GAIgAACAiyEAAgAAuBgCIAAAgIshAAIAALgYAiAAAICLIQACAAC4GAIgAACAiyEAAgAAuBgCIAAAgIspcgFw6tSpqlKlinx8fNSoUSNt2bIl17ZRUVGyWCzZltq1a9vaXL16VRMmTFC1atXk4+OjevXqac2aNYVxKgAAAE5RpALg4sWLNWzYMI0aNUqxsbFq0aKF2rdvr4SEhBzbT5kyRVar1bYkJiaqTJky6t69u63N6NGj9fHHH+vDDz/UwYMHNXjwYD3yyCOKjY0trNMCAAAoVBZjjHF2EXnVpEkTNWzYUNOmTbOtq1Wrlrp27aro6Ogb7r9ixQp169ZN8fHxioiIkCSFhYVp1KhRGjJkiK1d165dVbJkSc2fPz9PdaWmpiogIEApKSny9/fP51kBAADkz61mjyIzA5ienq5du3apbdu2duvbtm2rbdu25amPmTNnqk2bNrbwJ0lpaWny8fGxa+fr66utW7feetEAAAC3IQ9nF5BXZ86cUWZmpoKDg+3WBwcHKykp6Yb7W61WrV69WgsXLrRb365dO02aNEktW7ZUtWrVtGHDBv33v/9VZmZmrn2lpaUpLS3N9jk1NTWfZwMAAOA8RWYG8BqLxWL32RiTbV1O5syZo8DAQHXt2tVu/ZQpU1S9enXVrFlTXl5eevbZZ9W/f3+5u7vn2ld0dLQCAgJsS3h4+E2dCwAAgDMUmQBYrlw5ubu7Z5vtS05OzjYr+HfGGM2aNUt9+vSRl5eX3bagoCCtWLFCFy9e1IkTJ/Tzzz+rZMmSqlKlSq79jRw5UikpKbYlMTHx5k8MAACgkBWZAOjl5aVGjRopJibGbn1MTIyaNWt23X03b96so0ePauDAgbm28fHxUYUKFZSRkaGlS5eqS5cuubb19vaWv7+/3QIAAFBUFJl7ACVp+PDh6tOnjxo3bqymTZvqk08+UUJCggYPHizpz5m5kydPat68eXb7zZw5U02aNFFkZGS2Pn/88UedPHlS9evX18mTJzVu3DhlZWXpX//6V6GcEwAAQGErUgGwR48eOnv2rCZMmCCr1arIyEitWrXK9lSv1WrN9k7AlJQULV26VFOmTMmxzytXrmj06NE6duyYSpYsqQ4dOuizzz5TYGBgQZ8OAACAUxSp9wDerngPIAAAKEwu8x5AAAAAOAYBEAAAwMUQAAEAAFwMARAAAMDFEAABAABcDAEQAADAxRAAAQAAXAwBEAAAwMUQAAEAAFwMARAAAMDFOOS7gFeuXJnvfR566CH5+vo64vAAAADIB4cEwK5du+arvcVi0ZEjR1S1alVHHB4AAAD54LBLwElJScrKysrT4ufn56jDAgAAIJ8cEgD79euXr8u5vXv3lr+/vyMODQAAgHyyGGOMs4so6lJTUxUQEKCUlBSCLQAAKHC3mj0c/hTw5cuXdenSJdvnEydOaPLkyVq3bp2jDwUAAICb4PAA2KVLF82bN0+SdP78eTVp0kTvvfeeunTpomnTpjn6cAAAAMgnhwfA3bt3q0WLFpKkL7/8UsHBwTpx4oTmzZunDz74wNGHAwAAQD45PABeunRJpUqVkiStW7dO3bp1k5ubm+69916dOHHC0YcDAABAPjk8AN5xxx1asWKFEhMTtXbtWrVt21aSlJyczAMSAAAAtwGHB8AxY8bopZdeUuXKldWkSRM1bdpU0p+zgQ0aNLjl/qdOnaoqVarIx8dHjRo10pYtW3JtGxUVJYvFkm2pXbu2XbvJkyerRo0a8vX1VXh4uF544QVduXLllmsFAAC4HRXIa2CSkpJktVpVr149ubn9mTG3b98uf39/1axZ86b7Xbx4sfr06aOpU6eqefPm+vjjj/Xpp5/q4MGDqlSpUrb2KSkpunz5su1zRkaG6tWrp+eee07jxo2TJC1YsEADBw7UrFmz1KxZMx0+fFhRUVHq0aOH3n///TzVxWtgAABAYbrV7OGwAPjqq6+qa9euuueeexzRXY6aNGmihg0b2j1NXKtWLXXt2lXR0dE33H/FihXq1q2b4uPjFRERIUl69tln9dNPP2nDhg22di+++KK2b99+3dnFvyIAAgCAwnTbvAfQarWqU6dOCg0N1VNPPaVvvvlGaWlpjupe6enp2rVrl+2ewmvatm2rbdu25amPmTNnqk2bNrbwJ0n33Xefdu3ape3bt0uSjh07plWrVqljx4659pOWlqbU1FS7BQAAoKhwWACcPXu2Tp8+rS+++EKBgYF68cUXVa5cOXXr1k1z5szRmTNnbqn/M2fOKDMzU8HBwXbrg4ODlZSUdMP9rVarVq9erUGDBtmt79mzp15//XXdd9998vT0VLVq1dS6dWuNGDEi176io6MVEBBgW8LDw2/upAAAAJzAoQ+BWCwWtWjRQm+//bZ+/vlnbd++Xffee69mzJihChUqqGXLlnr33Xd18uTJWzrGXxljsq3LyZw5cxQYGKiuXbvard+0aZMmTpyoqVOnavfu3Vq2bJm+/vprvf7667n2NXLkSKWkpNiWxMTEmzoXAAAAZ/AoyM5r1aqlWrVq6V//+pd+++03rVy5UitXrpQkvfTSS/nqq1y5cnJ3d88225ecnJxtVvDvjDGaNWuW+vTpIy8vL7ttr732mvr06WObGaxTp44uXryop556SqNGjbI9xPJX3t7e8vb2zlf9AAAAtwuHvwYmN0FBQRo4cKD++9//5jv8SZKXl5caNWqkmJgYu/UxMTFq1qzZdffdvHmzjh49qoEDB2bbdunSpWwhz93dXcYYFcAD0gAAAE53yzOA586dkzFGZcqU0W+//abvvvtONWrUUGRkpCPqszN8+HD16dNHjRs3VtOmTfXJJ58oISFBgwcPlvTnpdmTJ0/avov4mpkzZ6pJkyY51tS5c2dNmjRJDRo0UJMmTXT06FG99tprevjhh+Xu7u7wcwAAAHC2WwqAn376qaKjo5WVlaV//etfWrBggerWrauxY8fq+eef11NPPeWoOiVJPXr00NmzZzVhwgRZrVZFRkZq1apVtqd6rVarEhIS7PZJSUnR0qVLNWXKlBz7HD16tCwWi0aPHq2TJ08qKChInTt31sSJEx1aOwAAwO3ilt4DWK9ePf3444+6dOmSKlWqpPj4eAUFBSk1NVUtW7ZUXFycA0u9ffEeQAAAUJhuNXvc0gygu7u7fHx85OPjozvuuENBQUGSJH9//zw9mQsAAIDCd0sPgXh4eNi+M3fz5s229RcuXLi1qgAAAFBgbikAfvvtt7bXoQQEBNjWX758WTNnzry1ygAAAFAgbukScMmSJXNc7+/vr4yMDH399dfKysqy2/bwww/fyiEBAABwixz+Iug1a9aob9++OX71m8ViUWZmpqMPCQAAgHxw+Iugn332WXXv3l1Wq1VZWVl2C+EPAADA+RweAJOTkzV8+PAbfj0bAAAAnMPhAfDRRx/Vpk2bHN0tAAAAHOSWXgSdk0uXLql79+4KCgpSnTp15Onpabf9+eefd+Thbgu8CBoAABQmp74IOicLFy7U2rVr5evrq02bNtm9ENpisRTLAAgAAFCUODwAjh49WhMmTNCIESPk5ubwK8wAAAC4RQ5PaOnp6erRowfhDwAA4Dbl8JTWr18/LV682NHdAgAAwEEcfgk4MzNTb7/9ttauXau6detmewhk0qRJjj4kAAAA8sHhAXDfvn1q0KCBJGn//v122/76QAjy7o2vD2rr0ezfrAIAAG5fr7SvqdY1yju7jBw5PABu3LjR0V26vFMpl/Vz0gVnlwEAAPLhwpUMZ5eQK4cHQDje8w9W1xNNIpxdBgAAyIfqwSWdXUKuHBIA9+7dq8jIyDw/+XvgwAHVqFFDHh7kz7yoGcLLpQEAgOM45CngBg0a6OzZs3lu37RpUyUkJDji0AAAAMgnh0zBGWP02muvyc/PL0/t09PTb/pYU6dO1TvvvCOr1aratWtr8uTJatGiRY5to6KiNHfu3Gzr77rrLh04cECSdP/992vz5s3Z2nTo0EHffPPNTdcJAABwu3JIAGzZsqUOHTqU5/ZNmzaVr69vvo+zePFiDRs2TFOnTlXz5s318ccfq3379jp48KAqVaqUrf2UKVP05ptv2j5nZGSoXr166t69u23dsmXL7ALp2bNns7UBAAAoTizGGOPsIvKqSZMmatiwoaZNm2ZbV6tWLXXt2lXR0dE33H/FihXq1q2b4uPjFRGR80MVkydP1pgxY2S1WlWiRIk81XWrX8gMAACQH7eaPYrM97Wlp6dr165datu2rd36tm3batu2bXnqY+bMmWrTpk2u4e9am549e143/KWlpSk1NdVuAQAAKCqKTAA8c+aMMjMzFRwcbLc+ODhYSUlJN9zfarVq9erVGjRoUK5ttm/frv3791+3jSRFR0crICDAtoSHh+ftJAAAAG4DRSYAXvP3bxMxxuTpG0bmzJmjwMBAde3aNdc2M2fOVGRkpO65557r9jVy5EilpKTYlsTExDzVDgAAcDsolBfxpaWlydvb+5b6KFeunNzd3bPN9iUnJ2ebFfw7Y4xmzZqlPn36yMvLK8c2ly5d0qJFizRhwoQb1uLt7X3L5wMAAOAshTID2KxZs2zrDh8+nK8+vLy81KhRI8XExNitj4mJybH/v9q8ebOOHj2qgQMH5trmiy++UFpamnr37p2vugAAAIqaAp0B/Prrr/Xzzz/r4sWLOnXqlMLCwmzbunfvrj179uSrv+HDh6tPnz5q3LixmjZtqk8++UQJCQkaPHiwpD8vzZ48eVLz5s2z22/mzJlq0qSJIiMjc+175syZ6tq1q8qWLZuvmgAAAIqaAg2AtWvXVkJCgpKTk9WrVy8lJiaqYsWKCgsLk7u7e77769Gjh86ePasJEybIarUqMjJSq1atsj3Va7Vas33DSEpKipYuXaopU6bk2u/hw4e1detWrVu3Lt81AQAAFDWF8h7A7777Ti1btpQknTx5UvHx8YqMjFRgYGBBH7pQ8B5AAABQmG41exTKQyDXwp8kVahQQRUqVCiMwwIAACAHhRIAo6KiFBkZqdq1aysyMpL35gEAADhRoTwF/NRTT6lkyZL66quv9M9//lOBgYFq2rRpYRwaAAAAf1MoM4DNmjWze1XLd999pw0bNhTGoQEAAPA3hTIDmJKSYve5ZcuW+uWXXwrj0AAAAPibQnsI5OLFi7rzzjsVGRkpHx8f7d27tzAODQAAgL8plAC4Z88eZWZm6tChQ9q/f79+//13rVy5sjAODQAAgL8plACYkZGhRYsW6bffftNdd92l7t27y2KxFMahAQAA8DeFcg9gr169tHXrVlksFn355Zdq0KBBvr8LGAAAAI5RKDOAhw4dsrvnb/fu3Xrqqae0adOmwjg8AAAA/qJQZgBLlixp99Rvw4YN9fvvvxfGoQEAAPA3hTID+PHHH6tr165q3769atWqpZ9++kmVKlUqjEMDAADgbxw+A3j48GEtWbJEy5cv17FjxyRJderU0c6dO9WoUSOdOHFC1apV0xdffOHoQwMAACAPHDYDmJGRof79+2vhwoUyxkiSLBaLmjdvrg8++ED169dXjx49HHU4AAAA3CSHzQBOnDhRq1at0owZM/TLL79o//79mjNnjtLT09WiRQutX7/eUYcCAADALbCYa9N1t+iOO+7QmDFj1Ldv32zb3nvvPY0bN05HjhyRr6+vdu/erdatWzvisLeF1NRUBQQEKCUlRf7+/s4uBwAAFHO3mj0cNgOYmJioFi1a5LjtxRdfVM+ePTVw4EA1atRIP/zwg6MOCwAAgHxyWAAsU6aMzp07l+v2QYMGafXq1WrTpo1efPFFRx0WAAAA+eSwAHj//fdr/vz5uW4PDg6Wh4eHpk+fLi8vL0cdFgAAAPnksAD4yiuv6KOPPso1BO7cuVMVK1a85eNMnTpVVapUkY+Pjxo1aqQtW7bk2jYqKkoWiyXbUrt2bbt258+f15AhQxQaGiofHx/VqlVLq1atuuVaAQAAbkcOC4D169fXtGnTFBUVpS5dumjdunU6ffq0UlJStHLlSr3wwgu3/BqYxYsXa9iwYRo1apRiY2PVokULtW/fXgkJCTm2nzJliqxWq21JTExUmTJl1L17d1ub9PR0PfTQQzp+/Li+/PJLHTp0SDNmzFCFChVuqVYAAIDblcOeAr5m8+bNGj58uGJjY2WxWCRJxhj94x//0LJly+Tj43PTfTdp0kQNGzbUtGnTbOtq1aqlrl27Kjo6+ob7r1ixQt26dVN8fLwiIiIkSdOnT9c777yjn3/+WZ6enjdVF08BAwCAwnSr2cPhAfCa/fv3Ky4uTunp6apbt64aN258S/2lp6fLz89PS5Ys0SOPPGJbP3ToUMXFxWnz5s037KNz585KS0vTunXrbOs6dOigMmXKyM/PT//9738VFBSkxx9/XK+88orc3d3zVBsBEAAAFKZbzR4F9l3AkZGRioyMdFh/Z86cUWZmpoKDg+3WBwcHKykp6Yb7W61WrV69WgsXLrRbf+zYMX377bd64okntGrVKh05ckRDhgxRRkaGxowZk2NfaWlpSktLs31OTU29iTMCAABwDod/F3BBu3ZZ+RpjTLZ1OZkzZ44CAwPVtWtXu/VZWVkqX768PvnkEzVq1Eg9e/bUqFGj7C4z/110dLQCAgJsS3h4+E2dCwAAgDMUmQBYrlw5ubu7Z5vtS05OzjYr+HfGGM2aNUt9+vTJ9gqa0NBQ3XnnnXaXe2vVqqWkpCSlp6fn2N/IkSOVkpJiWxITE2/yrAAAAApfkQmAXl5eatSokWJiYuzWx8TEqFmzZtfdd/PmzTp69KgGDhyYbVvz5s119OhRZWVl2dYdPnxYoaGhub6v0NvbW/7+/nYLAABAUVFkAqAkDR8+XJ9++qlmzZqln376SS+88IISEhI0ePBgSX/OzOX0XcQzZ85UkyZNcrwn8f/+7/909uxZDR06VIcPH9Y333yjf//73xoyZEiBnw8AAIAzFNhDIAWhR48eOnv2rCZMmCCr1arIyEitWrXK9koXq9Wa7Z2AKSkpWrp0qaZMmZJjn+Hh4Vq3bp1eeOEF1a1bVxUqVNDQoUP1yiuvFPj5AAAAOEOBvQbGlfAaGAAAUJhuNXsUqUvAAAAAuHUEQAAAABdDAAQAAHAxBEAAAAAXQwAEAABwMQRAAAAAF0MABAAAcDEEQAAAABdDAAQAAHAxBEAAAAAXQwAEAABwMQRAAAAAF0MABAAAcDEEQAAAABdDAAQAAHAxBEAAAAAXQwAEAABwMQRAAAAAF0MABAAAcDEEQAAAABdT5ALg1KlTVaVKFfn4+KhRo0basmVLrm2joqJksViyLbVr17a1mTNnTo5trly5UhinAwAAUOiKVABcvHixhg0bplGjRik2NlYtWrRQ+/btlZCQkGP7KVOmyGq12pbExESVKVNG3bt3t2vn7+9v185qtcrHx6cwTgkAAKDQFakAOGnSJA0cOFCDBg1SrVq1NHnyZIWHh2vatGk5tg8ICFBISIht2blzp86dO6f+/fvbtbNYLHbtQkJCCuN0AAAAnKLIBMD09HTt2rVLbdu2tVvftm1bbdu2LU99zJw5U23atFFERITd+j/++EMRERGqWLGiOnXqpNjY2Ov2k5aWptTUVLsFAACgqCgyAfDMmTPKzMxUcHCw3frg4GAlJSXdcH+r1arVq1dr0KBBdutr1qypOXPmaOXKlfr888/l4+Oj5s2b68iRI7n2FR0drYCAANsSHh5+cycFAADgBEUmAF5jsVjsPhtjsq3LyZw5cxQYGKiuXbvarb/33nvVu3dv1atXTy1atNAXX3yhO++8Ux9++GGufY0cOVIpKSm2JTEx8abOBQAAwBk8nF1AXpUrV07u7u7ZZvuSk5OzzQr+nTFGs2bNUp8+feTl5XXdtm5ubrr77ruvOwPo7e0tb2/vvBcPAABwGykyM4BeXl5q1KiRYmJi7NbHxMSoWbNm19138+bNOnr0qAYOHHjD4xhjFBcXp9DQ0FuqFwAA4HZVZGYAJWn48OHq06ePGjdurKZNm+qTTz5RQkKCBg8eLOnPS7MnT57UvHnz7PabOXOmmjRposjIyGx9jh8/Xvfee6+qV6+u1NRUffDBB4qLi9NHH31UKOcEAABQ2IpUAOzRo4fOnj2rCRMmyGq1KjIyUqtWrbI91Wu1WrO9EzAlJUVLly7VlClTcuzz/Pnzeuqpp5SUlKSAgAA1aNBA3333ne65554CPx8AAABnsBhjjLOLKOpSU1MVEBCglJQU+fv7O7scAABQzN1q9igy9wACAADAMQiAAAAALoYACAAA4GIIgAAAAC6GAAgAAOBiCIAAAAAuhgAIAADgYgiAAAAALoYACAAA4GIIgAAAAC6GAAgAAOBiCIAAAAAuhgAIAADgYgiAAAAALoYACAAA4GIIgAAAAC6GAAgAAOBiCIAAAAAuhgAIAADgYopcAJw6daqqVKkiHx8fNWrUSFu2bMm1bVRUlCwWS7aldu3aObZftGiRLBaLunbtWkDVAwAAOF+RCoCLFy/WsGHDNGrUKMXGxqpFixZq3769EhIScmw/ZcoUWa1W25KYmKgyZcqoe/fu2dqeOHFCL730klq0aFHQpwEAAOBURSoATpo0SQMHDtSgQYNUq1YtTZ48WeHh4Zo2bVqO7QMCAhQSEmJbdu7cqXPnzql///527TIzM/XEE09o/Pjxqlq1amGcCgAAgNMUmQCYnp6uXbt2qW3btnbr27Ztq23btuWpj5kzZ6pNmzaKiIiwWz9hwgQFBQVp4MCBeeonLS1NqampdgsAAEBR4eHsAvLqzJkzyszMVHBwsN364OBgJSUl3XB/q9Wq1atXa+HChXbr//e//2nmzJmKi4vLcy3R0dEaP358ntsDAADcTorMDOA1FovF7rMxJtu6nMyZM0eBgYF2D3hcuHBBvXv31owZM1SuXLk81zBy5EilpKTYlsTExDzvCwAA4GxFZgawXLlycnd3zzbbl5ycnG1W8O+MMZo1a5b69OkjLy8v2/pffvlFx48fV+fOnW3rsrKyJEkeHh46dOiQqlWrlq0/b29veXt738rpAAAAOE2RmQH08vJSo0aNFBMTY7c+JiZGzZo1u+6+mzdv1tGjR7Pd41ezZk3t27dPcXFxtuXhhx9W69atFRcXp/DwcIefBwAAgLMVmRlASRo+fLj69Omjxo0bq2nTpvrkk0+UkJCgwYMHS/rz0uzJkyc1b948u/1mzpypJk2aKDIy0m69j49PtnWBgYGSlG09AABAcVGkAmCPHj109uxZTZgwQVarVZGRkVq1apXtqV6r1ZrtnYApKSlaunSppkyZ4oySAQAAbjsWY4xxdhFFXWpqqgICApSSkiJ/f39nlwMAAIq5W80eReYeQAAAADgGARAAAMDFEAABAABcDAEQAADAxRAAAQAAXAwBEAAAwMUQAAEAAFwMARAAAMDFEAABAABcDAEQAADAxRAAAQAAXAwBEAAAwMUQAAEAAFwMARAAAMDFEAABAABcDAEQAADAxRAAAQAAXAwBEAAAwMUQAAEAAFwMARAAAMDFEAABAABcDAEQAADAxXg4u4DiwBgjSUpNTXVyJQAAwBVcyxzXMkh+EQAd4MKFC5Kk8PBwJ1cCAABcyYULFxQQEJDv/SzmZqMjbLKysnTq1CmVKlVKFovF4f2npqYqPDxciYmJ8vf3d3j/roSxdAzG0TEYR8dhLB2DcXSMwhhHY4wuXLigsLAwubnl/44+ZgAdwM3NTRUrVizw4/j7+/MH0kEYS8dgHB2DcXQcxtIxGEfHKOhxvJmZv2t4CAQAAMDFEAABAABcDAGwCPD29tbYsWPl7e3t7FKKPMbSMRhHx2AcHYexdAzG0TGKwjjyEAgAAICLYQYQAADAxRAAAQAAXAwBEAAAwMUQAAEAAFwMAbAImDp1qqpUqSIfHx81atRIW7ZscXZJt5XvvvtOnTt3VlhYmCwWi1asWGG33RijcePGKSwsTL6+vrr//vt14MABuzZpaWl67rnnVK5cOZUoUUIPP/ywfv3110I8C+eLjo7W3XffrVKlSql8+fLq2rWrDh06ZNeGsbyxadOmqW7durYXwDZt2lSrV6+2bWcMb050dLQsFouGDRtmW8dY5s24ceNksVjslpCQENt2xjHvTp48qd69e6ts2bLy8/NT/fr1tWvXLtv2IjWWBre1RYsWGU9PTzNjxgxz8OBBM3ToUFOiRAlz4sQJZ5d221i1apUZNWqUWbp0qZFkli9fbrf9zTffNKVKlTJLly41+/btMz169DChoaEmNTXV1mbw4MGmQoUKJiYmxuzevdu0bt3a1KtXz2RkZBTy2ThPu3btzOzZs83+/ftNXFyc6dixo6lUqZL5448/bG0YyxtbuXKl+eabb8yhQ4fMoUOHzKuvvmo8PT3N/v37jTGM4c3Yvn27qVy5sqlbt64ZOnSobT1jmTdjx441tWvXNlar1bYkJyfbtjOOefP777+biIgIExUVZX788UcTHx9v1q9fb44ePWprU5TGkgB4m7vnnnvM4MGD7dbVrFnTjBgxwkkV3d7+HgCzsrJMSEiIefPNN23rrly5YgICAsz06dONMcacP3/eeHp6mkWLFtnanDx50ri5uZk1a9YUWu23m+TkZCPJbN682RjDWN6K0qVLm08//ZQxvAkXLlww1atXNzExMaZVq1a2AMhY5t3YsWNNvXr1ctzGOObdK6+8Yu67775ctxe1seQS8G0sPT1du3btUtu2be3Wt23bVtu2bXNSVUVLfHy8kpKS7MbQ29tbrVq1so3hrl27dPXqVbs2YWFhioyMdOlxTklJkSSVKVNGEmN5MzIzM7Vo0SJdvHhRTZs2ZQxvwpAhQ9SxY0e1adPGbj1jmT9HjhxRWFiYqlSpop49e+rYsWOSGMf8WLlypRo3bqzu3burfPnyatCggWbMmGHbXtTGkgB4Gztz5owyMzMVHBxstz44OFhJSUlOqqpouTZO1xvDpKQkeXl5qXTp0rm2cTXGGA0fPlz33XefIiMjJTGW+bFv3z6VLFlS3t7eGjx4sJYvX6677rqLMcynRYsWadeuXYqOjs62jbHMuyZNmmjevHlau3atZsyYoaSkJDVr1kxnz55lHPPh2LFjmjZtmqpXr661a9dq8ODBev755zVv3jxJRe930qNQj4abYrFY7D4bY7Ktw/XdzBi68jg/++yz2rt3r7Zu3ZptG2N5YzVq1FBcXJzOnz+vpUuXql+/ftq8ebNtO2N4Y4mJiRo6dKjWrVsnHx+fXNsxljfWvn1723/XqVNHTZs2VbVq1TR37lzde++9khjHvMjKylLjxo3173//W5LUoEEDHThwQNOmTVPfvn1t7YrKWDIDeBsrV66c3N3ds/2rIDk5Odu/MJCza0+6XW8MQ0JClJ6ernPnzuXaxpU899xzWrlypTZu3KiKFSva1jOWeefl5aU77rhDjRs3VnR0tOrVq6cpU6Ywhvmwa9cuJScnq1GjRvLw8JCHh4c2b96sDz74QB4eHraxYCzzr0SJEqpTp46OHDnC72Q+hIaG6q677rJbV6tWLSUkJEgqev+PJADexry8vNSoUSPFxMTYrY+JiVGzZs2cVFXRUqVKFYWEhNiNYXp6ujZv3mwbw0aNGsnT09OujdVq1f79+11qnI0xevbZZ7Vs2TJ9++23qlKlit12xvLmGWOUlpbGGObDgw8+qH379ikuLs62NG7cWE888YTi4uJUtWpVxvImpaWl6aefflJoaCi/k/nQvHnzbK/GOnz4sCIiIiQVwf9HFuojJ8i3a6+BmTlzpjl48KAZNmyYKVGihDl+/LizS7ttXLhwwcTGxprY2FgjyUyaNMnExsbaXpXz5ptvmoCAALNs2TKzb98+06tXrxwfy69YsaJZv3692b17t3nggQdc7hUH//d//2cCAgLMpk2b7F4XcenSJVsbxvLGRo4cab777jsTHx9v9u7da1599VXj5uZm1q1bZ4xhDG/FX58CNoaxzKsXX3zRbNq0yRw7dsz88MMPplOnTqZUqVK2v0cYx7zZvn278fDwMBMnTjRHjhwxCxYsMH5+fmb+/Pm2NkVpLAmARcBHH31kIiIijJeXl2nYsKHttRz408aNG42kbEu/fv2MMX8+mj927FgTEhJivL29TcuWLc2+ffvs+rh8+bJ59tlnTZkyZYyvr6/p1KmTSUhIcMLZOE9OYyjJzJ4929aGsbyxAQMG2P68BgUFmQcffNAW/oxhDG/F3wMgY5k3195F5+npacLCwky3bt3MgQMHbNsZx7z76quvTGRkpPH29jY1a9Y0n3zyid32ojSWFmOMKdw5RwAAADgT9wACAAC4GAIgAACAiyEAAgAAuBgCIAAAgIshAAIAALgYAiAAAICLIQACAAC4GAIgAACAiyEAAgAAuBgCIAAUsvvvv18Wi0UWi0VxcXF52icqKsq2z4oVKwq0PgDFHwEQABxs2LBh6tq163XbPPnkk7JarYqMjMxTn1OmTJHVanVAdQBAAAQAh9uxY4fuueee67bx8/NTSEiIPDw88tRnQECAQkJCHFEeABAAAcBRrl69Ki8vL23btk2jRo2SxWJRkyZN8rz/l19+qTp16sjX11dly5ZVmzZtdPHixQKsGICryts/PQEAN+Tu7q6tW7eqSZMmiouLU3BwsHx8fPK0r9VqVa9evfT222/rkUce0YULF7RlyxYZYwq4agCuiAAIAA7i5uamU6dOqWzZsqpXr16+9rVarcrIyFC3bt0UEREhSapTp05BlAkAXAIGAEeKjY3Nd/iTpHr16unBBx9UnTp11L17d82YMUPnzp0rgAoBgAAIAA4VFxd3UwHQ3d1dMTExWr16te666y59+OGHqlGjhuLj4wugSgCujgAIAA60b98+1a1b96b2tVgsat68ucaPH6/Y2Fh5eXlp+fLlDq4QALgHEAAcKisrS3v37tWpU6dUokQJBQQE5Gm/H3/8URs2bFDbtm1Vvnx5/fjjj/rtt99Uq1atAq4YgCtiBhAAHOiNN97Q4sWLVaFCBU2YMCHP+/n7++u7775Thw4ddOedd2r06NF677331L59+wKsFoCrYgYQAByod+/e6t27d773q1WrltasWVMAFQFAdswAAoATTJ06VSVLltS+ffvy1H7w4MEqWbJkAVcFwFVYDG8ZBYBCdfLkSV2+fFmSVKlSJXl5ed1wn+TkZKWmpkqSQkNDVaJEiQKtEUDxRgAEAABwMVwCBgAAcDEEQAAAABdDAAQAAHAxBEAAAAAXQwAEAABwMQRAAAAAF0MABAAAcDEEQAAAABdDAAQAAHAxBEAAAAAXQwAEAABwMQRAAAAAF0MABAAAcDEEQAAAABfj4ewCioOsrCydOnVKpUqVksVicXY5AACgmDPG6MKFCwoLC5ObW/7n8wiADnDq1CmFh4c7uwwAAOBiEhMTVbFixXzvRwB0gFKlSkn684fg7+/v5GoAAEBxl5qaqvDwcFsGyS8CoANcu+zr7+9PAAQAAIXmZm894yEQAAAAF0MABAAAcDEEQAAAABfDPYAAUMRkZmbq6tWrzi4DQAFyd3eXh4dHgb1ejgAIAEXIH3/8oV9//VXGGGeXAqCA+fn5KTQ0VF5eXg7vmwAIAEVEZmamfv31V/n5+SkoKIgXzwPFlDFG6enp+u233xQfH6/q1avf1Muer4cACABFxNWrV2WMUVBQkHx9fZ1dDoAC5OvrK09PT504cULp6eny8fFxaP88BAIARQwzf4BrcPSsn13fBdYzAAAAbksEQACAU40bN07169d3aJ+bNm2SxWLR+fPnJUlz5sxRYGCgQ48B5IfFYtGKFSucXYYNARAAUKCioqJksVhksVjk6empqlWr6qWXXtLFixclSS+99JI2bNhQoDX06NFDhw8fLtBjIG9utyBUWKxWq9q3b+/sMmx4CAQAUOD+8Y9/aPbs2bp69aq2bNmiQYMG6eLFi5o2bZpKliypkiVLFujxfX19b9sHZ65evSpPT09nl2GHmhx/zJCQkEI7Vl4wAwgAKHDe3t4KCQlReHi4Hn/8cT3xxBO2WaC/XwKOiopS165dNX78eJUvX17+/v56+umnlZ6ebmtjjNHbb7+tqlWrytfXV/Xq1dOXX36Z6/H/fgn42jE/++wzVa5cWQEBAerZs6cuXLhw08eQpMqVK+v111/X448/rpIlSyosLEwffvihXRuLxaLp06erS5cuKlGihN544w1J0ldffaVGjRrJx8dHVatW1fjx45WRkWFXc6VKleTt7a2wsDA9//zztm1Tp05V9erV5ePjo+DgYD366KN2NU2ePNmuhvr162vcuHG3XNPf7dixQw899JDKlSungIAAtWrVSrt377arRZIeeeQRWSwW2+e/O378uCwWi7744gvdf//98vHx0fz58yVJs2fPVq1ateTj46OaNWtq6tSptv3S09P17LPPKjQ0VD4+PqpcubKio6Nt21NSUvTUU0/Zfq8eeOAB7dmzx26M69evr1mzZqlq1ary9vbWxx9/rAoVKigrK8uuxocfflj9+vWzfZ42bZqqVasmLy8v1ahRQ5999pld+7/OfF7v/AqNwS1LSUkxkkxKSoqzSwFQjF2+fNkcPHjQXL582RhjTFZWlrmYdtUpS1ZWVp7r7tevn+nSpYvduueee86ULVvWGGPM2LFjTb169ezalyxZ0vTo0cPs37/ffP311yYoKMi8+uqrtjavvvqqqVmzplmzZo355ZdfzOzZs423t7fZtGmTMcaYjRs3Gknm3LlzxhhjZs+ebQICAmz7jx071pQsWdJ069bN7Nu3z3z33XcmJCQkX8fISUREhClVqpSJjo42hw4dMh988IFxd3c369ats7WRZMqXL29mzpxpfvnlF3P8+HGzZs0a4+/vb+bMmWN++eUXs27dOlO5cmUzbtw4Y4wxS5YsMf7+/mbVqlXmxIkT5scffzSffPKJMcaYHTt2GHd3d7Nw4UJz/Phxs3v3bjNlyhS7mt5//327OuvVq2fGjh17SzXlZMOGDeazzz4zBw8eNAcPHjQDBw40wcHBJjU11RhjTHJyspFkZs+ebaxWq0lOTs6xn/j4eCPJVK5c2SxdutQcO3bMnDx50nzyyScmNDTUtm7p0qWmTJkyZs6cOcYYY9555x0THh5uvvvuO3P8+HGzZcsWs3DhQmPMn39emjdvbjp37mx27NhhDh8+bF588UVTtmxZc/bsWdvvRYkSJUy7du3M7t27zZ49e8yZM2eMl5eXWb9+va2+33//3Xh5eZm1a9caY4xZtmyZ8fT0NB999JE5dOiQee+994y7u7v59ttv7cZ4+fLl1z2/v/v7n/m/utXswSVgACiiLl/N1F1j1jrl2AcntJOf1839FbJ9+3YtXLhQDz74YK5tvLy8NGvWLPn5+al27dqaMGGCXn75Zb3++uu6fPmyJk2apG+//VZNmzaVJFWtWlVbt27Vxx9/rFatWuWpjqysLM2ZM0elSpWSJPXp00cbNmzQxIkTdfHixZs+RvPmzTVixAhJ0p133qn//e9/ev/99/XQQw/Z2jz++OMaMGCA7XOfPn00YsQI24xS1apV9frrr+tf//qXxo4dq4SEBIWEhKhNmzby9PRUpUqVdM8990iSEhISVKJECXXq1EmlSpVSRESEGjRokKcx+Kv81pSTBx54wO7zxx9/rNKlS2vz5s3q1KmTgoKCJEmBgYE3vCQqScOGDVO3bt1sn19//XW99957tnVVqlTRwYMH9fHHH6tfv35KSEhQ9erVdd9998lisSgiIsK278aNG7Vv3z4lJyfL29tbkvTuu+9qxYoV+vLLL/XUU09J+nMW8bPPPrPVKv15C8Nff2eXLFmiMmXK2D6/++67ioqK0jPPPCNJGj58uH744Qe9++67at26dZ7PrzARAAEABe7rr79WyZIllZGRoatXr6pLly7ZLo3+Vb169eTn52f73LRpU/3xxx9KTExUcnKyrly5YheopD//4s5P8KlcubIt/ElSaGiokpOTJUkHDx686WNcC4x//fz3S7CNGze2+7xr1y7t2LFDEydOtK3LzMzUlStXdOnSJXXv3l2TJ09W1apV9Y9//EMdOnRQ586d5eHhoYceekgRERG2bf/4xz/0yCOP2I1fXuS3ppz6T05O1pgxY/Ttt9/q9OnTyszM1KVLl5SQkJCvWnKq6bffflNiYqIGDhyoJ5980rY+IyNDAQEBkv68feChhx5SjRo19I9//EOdOnVS27Ztbefzxx9/qGzZsnbHuHz5sn755Rfb54iICLvwJ0lPPPGEnnrqKU2dOlXe3t5asGCBevbsKXd3d0nSTz/9ZAuQ1zRv3lxTpkzJ8/kVNgIgABRRvp7uOjihndOOnR+tW7fWtGnT5OnpqbCwsJu+2d1isdjuxfrmm29UoUIFu+3XZnby4u81/LVvRx3jr33/VYkSJew+Z2Vlafz48TnOBvn4+Cg8PFyHDh1STEyM1q9fr2eeeUbvvPOONm/erFKlSmn37t3atGmT1q1bpzFjxmjcuHHasWOHAgMD5ebmlu27o69evZrtOPmtKSdRUVH67bffNHnyZEVERMjb21tNmza1u38zP/5a07WfyYwZM9SkSRO7dteCWMOGDRUfH6/Vq1dr/fr1euyxx9SmTRt9+eWXysrKUmhoqDZt2pTtOH+9P/Tv4yBJnTt3VlZWlr755hvdfffd2rJliyZNmmTX5u8/Y2PMDV/antOxCgsBEACKKIvFctOXYQtbiRIldMcdd+S5/Z49e3T58mXbk7s//PCDSpYsqYoVK6p06dLy9vZWQkJCni/35tddd91108f44Ycfsn2uWbPmdfdp2LChDh06dN0x8vX11cMPP6yHH35YQ4YMUc2aNbVv3z41bNhQHh4eatOmjdq0aaOxY8cqMDBQ3377rbp166agoCBZrVZbP6mpqYqPj7/heeSlpr/bsmWLpk6dqg4dOkiSEhMTdebMGbs2np6eyszMzHOf1wQHB6tChQo6duyYnnjiiVzb+fv7q0ePHurRo4ceffRR/eMf/9Dvv/+uhg0bKikpSR4eHrk+fJIbX19fdevWTQsWLNDRo0d15513qlGjRrbttWrV0tatW9W3b1/bum3btqlWrVr5Ps/CUjT+z5FHGRkZGjdunBYsWKCkpCSFhoYqKipKo0ePzvXrVJYtW6Zp06YpLi5OaWlpql27tsaNG6d27Zzzr2oAwJ+XWgcOHKjRo0frxIkTGjt2rJ599lm5ubmpVKlSeumll/TCCy8oKytL9913n1JTU7Vt2zaVLFnS7snMm3Urx/jf//6nt99+W127dlVMTIyWLFmib7755rrHGzNmjDp16qTw8HB1795dbm5u2rt3r/bt26c33nhDc+bMUWZmppo0aSI/Pz999tln8vX1VUREhL7++msdO3ZMLVu2VOnSpbVq1SplZWWpRo0akv68L2/OnDnq3LmzSpcurddee802Y3YrNeXkjjvu0GeffabGjRsrNTVVL7/8crbX71SuXFkbNmxQ8+bN5e3trdKlS9+wlmvGjRun559/Xv7+/mrfvr3S0tK0c+dOnTt3TsOHD9f777+v0NBQ1a9fX25ublqyZIlCQkIUGBioNm3aqGnTpurataveeust1ahRQ6dOndKqVavUtWvXG16OfeKJJ9S5c2cdOHBAvXv3ttv28ssv67HHHlPDhg314IMP6quvvtKyZcu0fv36PJ9bobupR0duU2+88YYpW7as+frrr018fLxZsmSJKVmypJk8eXKu+wwdOtS89dZbZvv27ebw4cNm5MiRxtPT0+zevTvPx+UpYACF4XpPBN7OcnoK+K9yegq4S5cuZsyYMaZs2bKmZMmSZtCgQebKlSu2NllZWWbKlCmmRo0axtPT0wQFBZl27dqZzZs3G2Py9hTwX49pjDHvv/++iYiIyPMxchIREWHGjx9vHnvsMePn52eCg4Oz/R2kvzwN+ldr1qwxzZo1M76+vsbf39/cc889tid9ly9fbpo0aWL8/f1NiRIlzL333mt7KnXLli2mVatWpnTp0sbX19fUrVvXLF682NZvSkqKeeyxx4y/v78JDw83c+bMyfEp4PzWlJPdu3ebxo0bG29vb1O9enWzZMmSbE8hr1y50txxxx3Gw8PDbrz/6tpTsrGxsdm2LViwwNSvX994eXmZ0qVLm5YtW5ply5YZY4z55JNPTP369U2JEiWMv7+/efDBB+3+Pk9NTTXPPfecCQsLM56eniY8PNw88cQTJiEhwRiT8+/FNRkZGSY0NNRIMr/88ku27VOnTjVVq1Y1np6e5s477zTz5s2z264cngLO6fz+qiCfArb8v6KKhU6dOik4OFgzZ860rfvnP/9p+9dSXtWuXVs9evTQmDFj8tQ+NTVVAQEBSklJkb+/f77rBoC8uHLliuLj41WlSpVc78EqDqKionT+/Pki+W0RlStX1rBhwzRs2DBnl4Ji4Hp/5m81exSrF0Hfd9992rBhg+3rfvbs2aOtW7fa7kXIi6ysLF24cEFlypTJtU1aWppSU1PtFgAAgKKiWN0D+MorryglJUU1a9aUu7u7MjMzNXHiRPXq1SvPfbz33nu6ePGiHnvssVzbREdHa/z48Y4oGQAAoNAVqwC4ePFizZ8/XwsXLlTt2rUVFxenYcOGKSwsLE83BX/++ecaN26c/vvf/6p8+fK5ths5cqSGDx9u+5yamqrw8HCHnAMAuLo5c+Y4u4Sbdvz4cWeXAORJsQqAL7/8skaMGKGePXtKkurUqaMTJ04oOjr6hgFw8eLFGjhwoJYsWaI2bdpct623t/dNvQcKAADgdlCs7gG8dOlStte9uLu7Z/sC57/7/PPPFRUVpYULF6pjx44FWSIAAIDTFasZwM6dO2vixImqVKmSateurdjYWE2aNMnuuw1HjhypkydPat68eZL+DH99+/bVlClTdO+99yopKUnSny99vPbVMgBwOylGL28AcB0F+We9WM0Afvjhh3r00Uf1zDPPqFatWnrppZf09NNP6/XXX7e1sVqtdt9J+PHHHysjI0NDhgxRaGiobRk6dKgzTgEAcnXt5b03+7VaAIqWS5cuScr+tYWOUKzeA+gsvAcQQGEwxighIUFXr15VWFhYrt9wBKBoM8bo0qVLSk5OVmBgoEJDQ7O1udXsUawuAQNAcWaxWBQaGqr4+HidOHHC2eUAKGCBgYEKCQkpkL4JgABQhHh5eal69epcBgaKOU9Pzzx9Z/PNIgACQBHj5uZWrL8KDkDB4wYSAAAAF0MABAAAcDEEQAAAABdDAAQAAHAxBEAAAAAXQwAEAABwMQRAAAAAF0MABAAAcDEEQAAAABdDAAQAAHAxBEAAAAAXQwAEAABwMQRAAAAAF0MABAAAcDEEQAAAABdDAAQAAHAxxSoAZmRkaPTo0apSpYp8fX1VtWpVTZgwQVlZWbnuY7Va9fjjj6tGjRpyc3PTsGHDCq9gAAAAJ/BwdgGO9NZbb+n/a+/eo6qq8/+Pv45cjqCBeUFAUZC8BeiUdtGakCzNa42jMzqjglST3zQvTS0lK7UZB7PvWDY59s0L6TB99ZsiS0fTpBGVtLwioKU2XjCEyIpLJQeV/fujn2cNIyqXwz7Cfj7W2mu19/589n6fz7J6+Tl7f87bb7+tlStXKiIiQvv379eECRPk7++vqVOnVtnH4XCoTZs2mjVrll5//XWTKwYAADBfowqAe/bs0aOPPqohQ4ZIkkJDQ/W///u/2r9//zX7hIaGatGiRZKkFStWmFInAACAOzWqr4Dvv/9+ffTRRzp+/Lgk6fDhw8rIyNDgwYNdeh+Hw6GSkpJKGwAAQEPRqGYAZ8yYoeLiYnXr1k0eHh66fPmy5s2bpzFjxrj0PomJiZo7d65LrwkAAGCWRjUDuGbNGiUnJ+u9997TwYMHtXLlSv33f/+3Vq5c6dL7JCQkqLi42LmdPXvWpdcHAACoT41qBvD555/XzJkzNXr0aElSVFSUzpw5o8TERMXGxrrsPna7XXa73WXXAwAAMFOjmgH88ccf1aRJ5Y/k4eFx3WVgAAAArKZRzQAOGzZM8+bNU4cOHRQREaFDhw5p4cKFio+Pd7ZJSEhQXl6eVq1a5TyWmZkpSfr+++/19ddfKzMzU97e3rr99tvN/ggAAAD1zmYYhuHuIlyltLRUL730ktavX6/CwkIFBwdrzJgxevnll+Xt7S1JiouL0+nTp5Wenu7sZ7PZrrpWx44ddfr06Wrdt6SkRP7+/iouLpafn58rPgoAAMA11TV7NKoA6C4EQAAAYKa6Zo9G9QwgAAAAbowACAAAYDEEQAAAAIshAAIAAFgMARAAAMBiCIAAAAAWQwAEAACwGAIgAACAxRAAAQAALIYACAAAYDEEQAAAAIshAAIAAFgMARAAAMBiCIAAAAAWQwAEAACwGAIgAACAxRAAAQAALIYACAAAYDGNKgBeunRJL774osLCwuTj46NOnTrplVdeUUVFxXX77dixQ7169VLTpk3VqVMnvf322yZVDAAAYD5PM26yYcOGGvd5+OGH5ePjU6M+r776qt5++22tXLlSERER2r9/vyZMmCB/f39NnTq1yj6nTp3S4MGD9eSTTyo5OVkff/yxnn76abVp00a//OUva1w3AADAzc5mGIZR3zdp0qRmE402m00nTpxQp06datRv6NChatu2rZYvX+489stf/lK+vr7629/+VmWfGTNmaMOGDfrss8+cxyZOnKjDhw9rz5491bpvSUmJ/P39VVxcLD8/vxrVDAAAUFN1zR6mfQVcUFCgioqKam2+vr61usf999+vjz76SMePH5ckHT58WBkZGRo8ePA1++zZs0cDBgyodGzgwIHav3+/Ll68WKs6AAAAbmamfAUcGxtbo69zx44dW6s0O2PGDBUXF6tbt27y8PDQ5cuXNW/ePI0ZM+aafQoKCtS2bdtKx9q2batLly7p/PnzCgoKuqqPw+GQw+Fw7peUlNS4VgAAAHcxJQAmJSXVqP2SJUtqdZ81a9YoOTlZ7733niIiIpSZmalp06YpODhYsbGx1+xns9kq7V/5Vvw/j1+RmJiouXPn1qpGAAAAd2tUbwE///zzmjlzpkaPHq2oqCiNGzdO06dPV2Ji4jX7BAYGqqCgoNKxwsJCeXp6qlWrVlX2SUhIUHFxsXM7e/asSz8HAABAfTItAH766af64IMPKh1btWqVwsLCFBAQoN/97neVvlatjR9//PGqF048PDyuuwxMnz59tG3btkrHPvzwQ/Xu3VteXl5V9rHb7fLz86u0AQAANBSmBcA5c+YoKyvLuZ+dna3HH39cDz30kGbOnKmNGzded6auOoYNG6Z58+Zp06ZNOn36tNavX6+FCxfqF7/4hbNNQkKCxo8f79yfOHGizpw5o2effVafffaZVqxYoeXLl+u5556rUy0AAAA3K1OeAZSkzMxM/eEPf3Dur169Wvfcc4+WLl0qSQoJCdHs2bM1Z86cWt/jL3/5i1566SU9/fTTKiwsVHBwsJ566im9/PLLzjb5+fnKzc117oeFhWnz5s2aPn26Fi9erODgYL355pusAQgAABotU9YBlKSmTZvqxIkTCgkJkfTTki2PPPKIXnzxRUnS6dOnFRUVpdLSUjPKcSnWAQQAAGZqMOsAtm3bVqdOnZIklZeX6+DBg+rTp4/zfGlp6TWfuQMAAIDrmBYAH3nkEc2cOVO7du1SQkKCfH199fOf/9x5PisrS+Hh4WaVAwAAYFmmPQP4xz/+USNGjFB0dLSaN2+ud999V97e3s7zK1asuOoXOQAAAOB6pgXA4uJi7dq1S8XFxWrevLk8PDwqnX///ffVvHlzs8oBAACwLNMCYJcuXdSuXTvFxMTowQcfVL9+/RQaGuo837JlS7NKAQAAsDTTAuCOHTu0Y8cOpaena9KkSSorK1OHDh304IMPKiYmRjExMWrXrp1Z5QAAAFiWacvA/LuLFy9qz549Sk9PV3p6uj755BM5HA7ddtttOnbsmNnl1BnLwAAAADPVNXu4JQBeceHCBWVkZGjr1q1aunSpvv/+e12+fNld5dQaARAAAJiprtnDtK+AJamsrEy7d+/W9u3blZ6ern379iksLEzR0dFasmSJoqOjzSwHAADAkkwLgNHR0dq3b5/Cw8P1wAMP6JlnnlF0dLTatm1rVgkAAACQiQFw9+7dCgoKUkxMjPr166cHHnhArVu3Nuv2AAAA+P9M+yWQoqIivfPOO/L19dWrr76qdu3aKSoqSpMnT9batWv19ddfm1UKAACApbntJZDS0lJlZGQ4nwc8fPiwOnfurJycHHeUUye8BAIAAMxU1+xh2gzgf2rWrJlatmypli1b6tZbb5Wnp6c+++wzd5UDAABgGaY9A1hRUaH9+/crPT1d27dv18cff6wffvjB+esgixcvVkxMjFnlAAAAWJZpAbBFixb64YcfFBQUpH79+mnhwoWKiYlReHi4WSUAAABAJgbA1157TTExMerSpYtZtwQAAEAVTAuATz31lFm3AgAAwHWY+ksgV5SVlSkrK0uFhYWqqKiodG748OHuKAkAAMAyTA+AW7Zs0fjx43X+/Pmrztlstjr9FnBoaKjOnDlz1fGnn35aixcvrrLP4sWL9dZbb+n06dPq0KGDZs2apfHjx9e6BgAAgJud6cvATJ48WaNGjVJ+fr4qKioqbXUJf5K0b98+5efnO7dt27ZJkkaNGlVl+yVLlighIUFz5szRkSNHNHfuXE2aNEkbN26sUx0AAAA3M9MXgvbz89OhQ4dMeft32rRp+sc//qETJ07IZrNddb5v376677779Nprr1Xqs3//fmVkZFT7PiwEDQAAzNTgFoIeOXKk0tPT6/0+5eXlSk5OVnx8fJXhT5IcDoeaNm1a6ZiPj4/27t2rixcv1nuNAAAA7mD6M4BvvfWWRo0apV27dikqKkpeXl6Vzk+ZMsUl90lNTVVRUZHi4uKu2WbgwIFatmyZHnvsMd155506cOCAVqxYoYsXL+r8+fMKCgqqsp/D4ZDD4XDul5SUuKRmAAAAM5geAN977z1t3bpVPj4+Sk9PrzQ7Z7PZXBYAly9frkGDBik4OPiabV566SUVFBTo3nvvlWEYatu2reLi4rRgwQJ5eHhcs19iYqLmzp3rkjoBAADMZvozgIGBgZoyZYpmzpypJk3q5xvoM2fOqFOnTkpJSdGjjz56w/YXL17UV199paCgIL3zzjuaMWOGioqKrllfVTOAISEhPAMIAABMUddnAE2fASwvL9evf/3regt/kpSUlKSAgAANGTKkWu29vLzUvn17SdLq1as1dOjQ69Znt9tlt9tdUisAAIDZTH8JJDY2VmvWrKm361dUVCgpKUmxsbHy9KycbxMSEiqt8Xf8+HElJyfrxIkT2rt3r0aPHq2cnBz96U9/qrf6AAAA3M30GcDLly9rwYIF2rp1q3r06HHVSyALFy6s0/XT0tKUm5ur+Pj4q87l5+crNze3Ui1//vOfdezYMXl5eSkmJka7d+9WaGhonWoAAAC4mZn+DGBMTMw1z9lsNv3zn/80sRrXYB1AAABgpgb3DOD27dvNviUAAAD+jSnPAGZlZamioqLa7Y8cOaJLly7VY0UAAADWZUoAvOOOO/TNN99Uu32fPn0qPasHAAAA1zHlK2DDMPTSSy/J19e3Wu3Ly8vruSIAAADrMiUAPvDAAzp27Fi12/fp00c+Pj71WBEAAIB1mRIA09PTzbgNAAAAqsH0t4BRc5uz8/VF4ffuLgMAANTAI5GB6tL2FneXUSUCYAPwj6xz2pxd4O4yAABADYS2bkYARO3df1sbtfD1dncZAACgBjq2rN7Lr+5AAGwAfnNPB3eXAAAAGhFT1gEEAADAzcPtM4Dffvutli1bJi8vL02fPt3d5QAAADR6bp8BHDlypJo1a6Zly5ZJknJycjRr1iw3VwUAANB4uT0AlpaWatKkSfL2/uklh8jISG3evNnNVQEAADRebg+AAQEBOnfunGw2m/NYWVmZGysCAABo3Nz+DODrr7+u2NhYFRYWas2aNdqyZYu6devm7rIAAAAaLZthGIa7iygvL1dqaqqys7MVGBioCRMmyNf35l075z+VlJTI399fxcXF8vPzc3c5AACgkatr9nD7DODnn3+ujRs3qkWLFho0aJAiIyMbVPgDAABoaNz+DOCgQYNUXl6uoqIi/c///I/69eunrl27urssAACARsvtATAwMFCzZs3S888/r5UrV+rgwYM6evRora4VGhoqm8121TZp0qRr9vn73/+unj17ytfXV0FBQZowYYK++eab2n4cAACAm57bA+DAgQP1t7/9rdIxDw+PWl1r3759ys/Pd27btm2TJI0aNarK9hkZGRo/frwef/xxHTlyRO+//7727dunJ554olb3BwAAaAjc/gzg3r17tWLFCs2dO1d33323oqKiFBUVpaFDh9b4Wm3atKm0P3/+fIWHhys6OrrK9p988olCQ0M1ZcoUSVJYWJieeuopLViwoOYfBAAAoIFw+wzg5s2blZubq4MHD2ry5Mlq1aqV0tLS6nzd8vJyJScnKz4+vtIag/+ub9+++vLLL7V582YZhqGvvvpKa9eu1ZAhQ+p8fwAAgJuV25aBWbRokaZOnapjx46pc+fOatLEtVn0//7v//Sb3/xGubm5Cg4Ovma7tWvXasKECSorK9OlS5c0fPhwrV27Vl5eXtfs43A45HA4nPslJSUKCQlhGRgAAGCKui4D47YZwMjISEnS9OnT1bVrV915550aN26cXn31VW3atKnO11++fLkGDRp03fB39OhRTZkyRS+//LIOHDigLVu26NSpU5o4ceJ1r52YmCh/f3/nFhISUud6AQAAzHJTLAQt/ZRkc3JylJOTo6NHj+qNN96o9bXOnDmjTp06KSUlRY8++ug1240bN05lZWV6//33nccyMjL085//XOfOnVNQUFCV/ZgBBAAA7tRgF4J+9tln1aNHD/Xo0UMRERHy8/NT37591bdv3zpfOykpSQEBATd8lu/HH3+Up2flIbjyBvL1crHdbpfdbq9znQAAAO7gtgAYHR2trKwsbdq0SUeOHJGHh4ciIiKcobA2bwFLUkVFhZKSkhQbG3tVuEtISFBeXp5WrVolSRo2bJiefPJJLVmyRAMHDlR+fr6mTZumu++++7pfHQMAADRkbguAjz76aKWvZy9cuKCcnBxlZWUpLS2t1gEwLS1Nubm5io+Pv+pcfn6+cnNznftxcXEqLS3VW2+9pd///vdq0aKFHnzwQb366qu1ujcAAEBD4PZnAL/99lstW7ZM3t7emjZtmjtLqbW6fg8PAABQEw32LeArRo4cqWbNmmnp0qWSpJycHM2aNcvNVQEAADRebg+ApaWlmjRpkry9vSX9tDzM5s2b3VwVAABA4+X2ABgQEKBz585V+rWOsrIyN1YEAADQuLn9t4Bff/11xcbGqrCwUGvWrNGWLVvUrVs3d5cFAADQaLn9JRDpp9/tTU1NVXZ2tgIDAzVhwgT5+vq6u6xq4yUQAABgpga7EPQV2dnZeuONN/Tdd98pKipKw4cPb1DhDwAAoKFx+zOAI0eOVHR0tBISEhQcHKzhw4fro48+cndZAAAAjZbbZwD9/f01fvx4SdJdd92lESNG6KGHHtLhw4fdXBkAAEDj5PYZwE6dOmnhwoXO395t2bKlmjZt6uaqAAAAGi+3B0CHw6HFixerQ4cOeuSRRxQZGan+/fsrLy/P3aUBAAA0Sm57C3jRokWaOnWqjh07ps6dO+vChQvKysqqtJ07d07/+te/3FFejfAWMAAAMFODfQs4MjJSkjR9+nR98cUXat68uSIiIhQZGanBgwdr8eLF7ioNAACgUbsp1gGUfkqyOTk5ysnJ0dGjR/XGG2+4u6RqYwYQAACYqa7Z46YJgA0ZARAAAJipwX4FfEV2drZef/11FRUVKSoqSk888YRCQkLcXRYAAECj5fa3gEeOHKl+/fqxEDQAAIBJ3D4DyELQAAAA5nL7DCALQQMAAJjL7QGwrKyMhaABAABM5LYAmJeXp7y8PKWmpupf//qXPv/8c82ePVvTpk1TcXGxRo8erfDw8BpdMzQ0VDab7apt0qRJVbaPi4ursn1ERIQrPiIAAMBNyfRlYD7++GONHTtWubm5kqTWrVsrLi5Os2bNqvMSKl9//bUuX77s3M/JydHDDz+s7du3q1+/fle1Ly4u1oULF5z7ly5dUs+ePfXMM89ozpw51b4vy8AAAAAzNbhlYJ566ilFRERo3bp1stvtOnDggN58802lpKRoz549at26da2v3aZNm0r78+fPV3h4uKKjo6ts7+/vL39/f+d+amqqvvvuO02YMKHWNQAAANzsTJ8B9PHxUVZWljp37uw8ZhiGfvWrX8nLy0vvvfeeS+5TXl6u4OBgPfvss3rhhReq1WfYsGFyOBz68MMPr9vO4XDI4XA490tKShQSEsIMIAAAMEVdZwBNfwawe/fuKigoqHTMZrPplVde0caNG112n9TUVBUVFSkuLq5a7fPz8/XBBx/oiSeeuGHbxMRE5+yhv78/C1cDAIAGxfQAGBcXp9/97nfOZwCvKC4urvR1bF0tX75cgwYNUnBwcLXav/vuu2rRooUee+yxG7ZNSEhQcXGxczt79mwdqwUAADCP6c8ATps2TZLUpUsXjRgxQj/72c90+fJlJScn67XXXnPJPc6cOaO0tDSlpKRUq71hGFqxYoXGjRsnb2/vG7a32+2y2+11LRMAAMAtTA+ABQUFOnTokA4fPqzMzEy9++67OnHihGw2m+bPn69NmzapR48e6tGjhx555JFa3SMpKUkBAQEaMmRItdrv2LFDX3zxhR5//PFa3Q8AAKAhMf0lkKqUlZUpOztbmZmZzmCYk5OjoqKiGl+roqJCYWFhGjNmjObPn1/pXEJCgvLy8rRq1apKx8eNG6cTJ07ok08+qVX9LAMDAADM1OCWgalK06ZNddddd+muu+6q87XS0tKUm5ur+Pj4q87l5+dX+ezhunXrtGjRojrfGwAAoCG4KWYAGzpmAAEAgJka3DIwAAAAcC8CIAAAgMUQAAEAACyGAAgAAGAxBEAAAACLIQACAABYDAEQAADAYgiAAAAAFkMABAAAsBgCIAAAgMUQAAEAACyGAAgAAGAxBEAAAACLIQACAABYDAEQAADAYgiAAAAAFkMABAAAsBgCIAAAgMUQAAEAACymUQXA0NBQ2Wy2q7ZJkyZds4/D4dCsWbPUsWNH2e12hYeHa8WKFSZWDQAAYC5PdxfgSvv27dPly5ed+zk5OXr44Yc1atSoa/b51a9+pa+++krLly/XbbfdpsLCQl26dMmMcgEAANyiUQXANm3aVNqfP3++wsPDFR0dXWX7LVu2aMeOHTp58qRatmwp6adZRAAAgMasUX0F/O/Ky8uVnJys+Ph42Wy2Ktts2LBBvXv31oIFC9SuXTt16dJFzz33nC5cuHDdazscDpWUlFTaAAAAGopGNQP471JTU1VUVKS4uLhrtjl58qQyMjLUtGlTrV+/XufPn9fTTz+tb7/99rrPASYmJmru3Ln1UDUAAED9sxmGYbi7iPowcOBAeXt7a+PGjddsM2DAAO3atUsFBQXy9/eXJKWkpGjkyJH64Ycf5OPjU2U/h8Mhh8Ph3C8pKVFISIiKi4vl5+fn2g8CAADwH0pKSuTv71/r7NEoZwDPnDmjtLQ0paSkXLddUFCQ2rVr5wx/ktS9e3cZhqEvv/xSnTt3rrKf3W6X3W53ac0AAABmaZTPACYlJSkgIEBDhgy5brv77rtP586d0/fff+88dvz4cTVp0kTt27ev7zIBAADcotEFwIqKCiUlJSk2NlaenpUnOBMSEjR+/Hjn/m9+8xu1atVKEyZM0NGjR7Vz5049//zzio+Pv+bXvwAAAA1dowuAaWlpys3NVXx8/FXn8vPzlZub69xv3ry5tm3bpqKiIvXu3Vu//e1vNWzYML355ptmlgwAAGCqRvsSiJnq+iAmAABATdQ1ezS6GUAAAABcHwEQAADAYgiAAAAAFkMABAAAsBgCIAAAgMUQAAEAACyGAAgAAGAxBEAAAACLIQACAABYDAEQAADAYgiAAAAAFkMABAAAsBgCIAAAgMUQAAEAACyGAAgAAGAxBEAAAACLIQACAABYDAEQAADAYgiAAAAAFtOoAmBoaKhsNttV26RJk6psn56eXmX7zz//3OTKAQAAzOPp7gJcad++fbp8+bJzPycnRw8//LBGjRp13X7Hjh2Tn5+fc79Nmzb1ViMAAIC7NaoA+J/Bbf78+QoPD1d0dPR1+wUEBKhFixb1WBkAAMDNo1F9BfzvysvLlZycrPj4eNlstuu2veOOOxQUFKT+/ftr+/btN7y2w+FQSUlJpQ0AAKChaLQBMDU1VUVFRYqLi7tmm6CgIL3zzjtat26dUlJS1LVrV/Xv3187d+687rUTExPl7+/v3EJCQlxcPQAAQP2xGYZhuLuI+jBw4EB5e3tr48aNNeo3bNgw2Ww2bdiw4ZptHA6HHA6Hc7+kpEQhISEqLi6u9CwhAABAfSgpKZG/v3+ts0ejegbwijNnzigtLU0pKSk17nvvvfcqOTn5um3sdrvsdnttywMAAHCrRvkVcFJSkgICAjRkyJAa9z106JCCgoLqoSoAAICbQ6ObAayoqFBSUpJiY2Pl6Vn54yUkJCgvL0+rVq2SJL3xxhsKDQ1VRESE86WRdevWad26de4oHQAAwBSNLgCmpaUpNzdX8fHxV53Lz89Xbm6uc7+8vFzPPfec8vLy5OPjo4iICG3atEmDBw82s2QAAABTNdqXQMxU1wcxAQAAaqKu2aNRPgMIAACAayMAAgAAWAwBEAAAwGIIgAAAABZDAAQAALAYAiAAAIDFEAABAAAshgAIAABgMQRAAAAAiyEAAgAAWAwBEAAAwGIIgAAAABZDAAQAALAYAiAAAIDFEAABAAAsxtPdBTQGhmFIkkpKStxcCQAAsIIrmeNKBqkpAqALlJaWSpJCQkLcXAkAALCS0tJS+fv717ifzahtdIRTRUWFzp07p1tuuUU2m83l1y8pKVFISIjOnj0rPz8/l1/fShhL12AcXYNxdB3G0jUYR9cwYxwNw1BpaamCg4PVpEnNn+hjBtAFmjRpovbt29f7ffz8/PgX0kUYS9dgHF2DcXQdxtI1GEfXqO9xrM3M3xW8BAIAAGAxBEAAAACLIQA2AHa7XbNnz5bdbnd3KQ0eY+kajKNrMI6uw1i6BuPoGg1hHHkJBAAAwGKYAQQAALAYAiAAAIDFEAABAAAshgAIAABgMQTABuCvf/2rwsLC1LRpU/Xq1Uu7du1yd0k3lZ07d2rYsGEKDg6WzWZTampqpfOGYWjOnDkKDg6Wj4+P+vXrpyNHjlRq43A49Mwzz6h169Zq1qyZhg8fri+//NLET+F+iYmJuuuuu3TLLbcoICBAjz32mI4dO1apDWN5Y0uWLFGPHj2cC8D26dNHH3zwgfM8Y1g7iYmJstlsmjZtmvMYY1k9c+bMkc1mq7QFBgY6zzOO1ZeXl6exY8eqVatW8vX11c9+9jMdOHDAeb5BjaWBm9rq1asNLy8vY+nSpcbRo0eNqVOnGs2aNTPOnDnj7tJuGps3bzZmzZplrFu3zpBkrF+/vtL5+fPnG7fccouxbt06Izs72/j1r39tBAUFGSUlJc42EydONNq1a2ds27bNOHjwoBETE2P07NnTuHTpksmfxn0GDhxoJCUlGTk5OUZmZqYxZMgQo0OHDsb333/vbMNY3tiGDRuMTZs2GceOHTOOHTtmvPDCC4aXl5eRk5NjGAZjWBt79+41QkNDjR49ehhTp051Hmcsq2f27NlGRESEkZ+f79wKCwud5xnH6vn222+Njh07GnFxccann35qnDp1ykhLSzO++OILZ5uGNJYEwJvc3XffbUycOLHSsW7duhkzZ850U0U3t/8MgBUVFUZgYKAxf/5857GysjLD39/fePvttw3DMIyioiLDy8vLWL16tbNNXl6e0aRJE2PLli2m1X6zKSwsNCQZO3bsMAyDsayLW2+91Vi2bBljWAulpaVG586djW3bthnR0dHOAMhYVt/s2bONnj17VnmOcay+GTNmGPfff/81zze0seQr4JtYeXm5Dhw4oAEDBlQ6PmDAAO3evdtNVTUsp06dUkFBQaUxtNvtio6Odo7hgQMHdPHixUptgoODFRkZaelxLi4uliS1bNlSEmNZG5cvX9bq1av1ww8/qE+fPoxhLUyaNElDhgzRQw89VOk4Y1kzJ06cUHBwsMLCwjR69GidPHlSEuNYExs2bFDv3r01atQoBQQE6I477tDSpUud5xvaWBIAb2Lnz5/X5cuX1bZt20rH27Ztq4KCAjdV1bBcGafrjWFBQYG8vb116623XrON1RiGoWeffVb333+/IiMjJTGWNZGdna3mzZvLbrdr4sSJWr9+vW6//XbGsIZWr16tAwcOKDEx8apzjGX13XPPPVq1apW2bt2qpUuXqqCgQH379tU333zDONbAyZMntWTJEnXu3Flbt27VxIkTNWXKFK1atUpSw/sz6Wnq3VArNput0r5hGFcdw/XVZgytPM6TJ09WVlaWMjIyrjrHWN5Y165dlZmZqaKiIq1bt06xsbHasWOH8zxjeGNnz57V1KlT9eGHH6pp06bXbMdY3tigQYOc/xwVFaU+ffooPDxcK1eu1L333iuJcayOiooK9e7dW3/6058kSXfccYeOHDmiJUuWaPz48c52DWUsmQG8ibVu3VoeHh5X/a2gsLDwqr9hoGpX3nS73hgGBgaqvLxc33333TXbWMkzzzyjDRs2aPv27Wrfvr3zOGNZfd7e3rrtttvUu3dvJSYmqmfPnlq0aBFjWAMHDhxQYWGhevXqJU9PT3l6emrHjh1688035enp6RwLxrLmmjVrpqioKJ04cYI/kzUQFBSk22+/vdKx7t27Kzc3V1LD+28kAfAm5u3trV69emnbtm2Vjm/btk19+/Z1U1UNS1hYmAIDAyuNYXl5uXbs2OEcw169esnLy6tSm/z8fOXk5FhqnA3D0OTJk5WSkqJ//vOfCgsLq3Sesaw9wzDkcDgYwxro37+/srOzlZmZ6dx69+6t3/72t8rMzFSnTp0Yy1pyOBz67LPPFBQUxJ/JGrjvvvuuWhrr+PHj6tixo6QG+N9IU185QY1dWQZm+fLlxtGjR41p06YZzZo1M06fPu3u0m4apaWlxqFDh4xDhw4ZkoyFCxcahw4dci6VM3/+fMPf399ISUkxsrOzjTFjxlT5Wn779u2NtLQ04+DBg8aDDz5ouSUO/uu//svw9/c30tPTKy0X8eOPPzrbMJY3lpCQYOzcudM4deqUkZWVZbzwwgtGkyZNjA8//NAwDMawLv79LWDDYCyr6/e//72Rnp5unDx50vjkk0+MoUOHGrfccovz/yOMY/Xs3bvX8PT0NObNm2ecOHHC+Pvf/274+voaycnJzjYNaSwJgA3A4sWLjY4dOxre3t7GnXfe6VyWAz/Zvn27IemqLTY21jCMn17Nnz17thEYGGjY7XbjgQceMLKzsytd48KFC8bkyZONli1bGj4+PsbQoUON3NxcN3wa96lqDCUZSUlJzjaM5Y3Fx8c7/31t06aN0b9/f2f4MwzGsC7+MwAyltVzZS06Ly8vIzg42BgxYoRx5MgR53nGsfo2btxoREZGGna73ejWrZvxzjvvVDrfkMbSZhiGYe6cIwAAANyJZwABAAAshgAIAABgMQRAAAAAiyEAAgAAWAwBEAAAwGIIgAAAABZDAAQAALAYAiAAAIDFEAABAAAshgAIACbr16+fbDabbDabMjMzq9UnLi7O2Sc1NbVe6wPQ+BEAAcDFpk2bpscee+y6bZ588knl5+crMjKyWtdctGiR8vPzXVAdABAAAcDl9u3bp7vvvvu6bXx9fRUYGChPT89qXdPf31+BgYGuKA8ACIAA4CoXL16Ut7e3du/erVmzZslms+mee+6pdv+1a9cqKipKPj4+atWqlR566CH98MMP9VgxAKuq3l89AQA35OHhoYyMDN1zzz3KzMxU27Zt1bRp02r1zc/P15gxY7RgwQL94he/UGlpqXbt2iXDMOq5agBWRAAEABdp0qSJzp07p1atWqlnz5416pufn69Lly5pxIgR6tixoyQpKiqqPsoEAL4CBgBXOnToUI3DnyT17NlT/fv3V1RUlEaNGqWlS5fqu+++q4cKAYAACAAulZmZWasA6OHhoW3btumDDz7Q7bffrr/85S/q2rWrTp06VQ9VArA6AiAAuFB2drZ69OhRq742m0333Xef5s6dq0OHDsnb21vr1693cYUAwDOAAOBSFRUVysrK0rlz59SsWTP5+/tXq9+nn36qjz76SAMGDFBAQIA+/fRTff311+revXs9VwzAipgBBAAX+uMf/6g1a9aoXbt2euWVV6rdz8/PTzt37tTgwYPVpUsXvfjii/rzn/+sQYMG1WO1AKyKGUAAcKGxY8dq7NixNe7XvXt3bdmypR4qAoCrMQMIAG7w17/+Vc2bN1d2dna12k+cOFHNmzev56oAWIXNYJVRADBVXl6eLly4IEnq0KGDvL29b9insLBQJSUlkqSgoCA1a9asXmsE0LgRAAEAACyGr4ABAAAshgAIAABgMQRAAAAAiyEAAgAAWAwBEAAAwGIIgAAAABZDAAQAALAYAiAAAIDFEAABAAAshgAIAABgMQRAAAAAiyEAAgAAWAwBEAAAwGIIgAAAABZDAAQAALCY/wfK8VxqGCdOoQAAAABJRU5ErkJggg==", + "text/html": [ + "\n", + "
\n", + "
\n", + " Figure\n", + "
\n", + " \n", + "
\n", + " " + ], + "text/plain": [ + "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ - "%matplotlib qt5\n", + "%matplotlib widget\n", "fig1, (ax1, ax2, ax3) = plt.subplots(3, 1)\n", "fig1.set_figheight(10)\n", "fig1.suptitle('Ausgleichsbecken')\n", diff --git a/Druckrohrleitung/Druckrohrleitung_class_file.py b/Druckrohrleitung/Druckrohrleitung_class_file.py index 829ab55..c6ce967 100644 --- a/Druckrohrleitung/Druckrohrleitung_class_file.py +++ b/Druckrohrleitung/Druckrohrleitung_class_file.py @@ -9,14 +9,16 @@ sys.path.append(parent) from functions.pressure_conversion import pressure_conversion class Druckrohrleitung_class: -# units +# units + # make sure that units and display units are the same + # units are used to label graphs and disp units are used to have a bearable format when using pythons print() acceleration_unit = r'$\mathrm{m}/\mathrm{s}^2$' angle_unit = 'rad' area_unit = r'$\mathrm{m}^2$' density_unit = r'$\mathrm{kg}/\mathrm{m}^3$' flux_unit = r'$\mathrm{m}^3/\mathrm{s}$' length_unit = 'm' - pressure_unit = 'Pa' + pressure_unit = 'Pa' # DONT CHANGE needed for pressure conversion time_unit = 's' velocity_unit = r'$\mathrm{m}/\mathrm{s}$' # for flux and pressure propagation volume_unit = r'$\mathrm{m}^3$' @@ -27,33 +29,54 @@ class Druckrohrleitung_class: density_unit_disp = 'kg/m³' flux_unit_disp = 'm³/s' length_unit_disp = 'm' + # pressure_unit_disp will be set within the __init__() method time_unit_disp = 's' velocity_unit_disp = 'm/s' # for flux and pressure propagation volume_unit_disp = 'm³' - g = 9.81 + g = 9.81 # m/s² gravitational acceleration # init - def __init__(self,total_length,diameter,number_segments,pipeline_angle,Darcy_friction_factor,pw_vel,timestep,pressure_unit_disp,rho=1000): - self.length = total_length # total length of the pipeline - self.dia = diameter # diameter of the pipeline - self.n_seg = number_segments # number of segments for the method of characteristics - self.angle = pipeline_angle # angle of the pipeline - self.f_D = Darcy_friction_factor # = Rohrreibungszahl oder flow coefficient - self.c = pw_vel + def __init__(self,total_length,diameter,pipeline_head,number_segments,Darcy_friction_factor,pw_vel,timestep,pressure_unit_disp,rho=1000): + """ + Creates a reservoir with given attributes in this order: \n + Pipeline length [m] \n + Pipeline diameter [m] \n + Pipeline head [m] \n + Number of pipeline segments [1] \n + Darcy friction factor [1] \n + Pressure wave velocity [m/s] \n + Simulation timestep [s] \n + Pressure unit for displaying [string] \n + Density of the liquid [kg/m³] \n + """ + self.length = total_length # total length of the pipeline + self.dia = diameter # diameter of the pipeline + self.head = pipeline_head # hydraulic head of the pipeline + self.n_seg = number_segments # number of segments for the method of characteristics + self.f_D = Darcy_friction_factor # = Rohrreibungszahl oder flow coefficient + self.c = pw_vel # propagation velocity of pressure wave self.dt = timestep - self.density = rho # density of the liquid in the pipeline + self.density = rho # density of the liquid in the pipeline - self.A = (diameter/2)**2*np.pi + # derivatives of input attributes + self.angle = np.arcsin(self.head/self.length) # angle of the pipeline + self.A = (diameter/2)**2*np.pi # crossectional area of the pipeline + self.dx = total_length/number_segments # length of each segment + self.x_vec = np.arange(0,(number_segments+1),1)*self.dx # vector giving the distance from each node to the start of the pipeline + self.h_vec = np.arange(0,(number_segments+1),1)*self.head/self.n_seg # vector giving the height difference from each node to the start of the pipeline + + self.pressure_unit_disp = pressure_unit_disp # pressure unit for displaying - self.dx = total_length/number_segments # length of each segment - self.x_vec = np.arange(0,(number_segments+1),1)*self.dx # vector giving the distance from each node to the start of the pipeline - - self.pressure_unit_disp = pressure_unit_disp - -# setter - def set_initial_pressure(self,pressure): +# setter - set attributes + def set_initial_pressure(self,pressure,display_warning=True): # initialize the pressure distribution in the pipeline + if display_warning == True: + print('You are setting the pressure distribution in the pipeline manually. \n \ + This is not an intended use of this method. \n \ + Refer to the set_steady_state() method instead.') + + # make sure the vector has the right size if np.size(pressure) == 1: p0 = np.full_like(self.x_vec,pressure) elif np.size(pressure) == np.size(self.x_vec): @@ -64,11 +87,18 @@ class Druckrohrleitung_class: #initialize the vectors in which the old and new pressures are stored for the method of characteristics self.p_old = p0.copy() self.p = p0.copy() + # initialize the vectors in which the minimal and maximal value of the pressure at each node are stores self.p_min = p0.copy() self.p_max = p0.copy() - def set_initial_flow_velocity(self,velocity): + def set_initial_flow_velocity(self,velocity, display_warning=True): # initialize the velocity distribution in the pipeline + if display_warning == True: + print('You are setting the velocity distribution in the pipeline manually. \n \ + This is not an intended use of this method. \n \ + Refer to the set_steady_state() method instead.') + + # make sure the vector has the right size if np.size(velocity) == 1: v0 = np.full_like(self.x_vec,velocity) elif np.size(velocity) == np.size(self.x_vec): @@ -79,6 +109,7 @@ class Druckrohrleitung_class: #initialize the vectors in which the old and new velocities are stored for the method of characteristics self.v_old = v0.copy() self.v = v0.copy() + # initialize the vectors in which the minimal and maximal value of the velocity at each node are stores self.v_min = v0.copy() self.v_max = v0.copy() @@ -114,21 +145,19 @@ class Druckrohrleitung_class: self.p[0] = p_boundary_res self.p[-1] = p_boundary_tur - def set_steady_state(self,ss_flux,ss_level_reservoir,area_reservoir,x_vec,h_vec): + def set_steady_state(self,ss_flux,ss_pressure_res): # 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_like(self.x_vec,ss_flux/self.A) # the static pressure is given by static state pressure of the reservoir, corrected for the hydraulic head of the pipe and friction losses - ss_v_in_res = abs(ss_flux/area_reservoir) - ss_v_out_res = abs(ss_flux/self.A) - ss_pressure_res = self.density*self.g*(ss_level_reservoir)+self.density*ss_v_out_res*(ss_v_in_res-ss_v_out_res) - ss_pressure = ss_pressure_res+(self.density*self.g*h_vec)-(self.f_D*x_vec/self.dia*self.density/2*ss_v0**2) + ss_pressure = ss_pressure_res+(self.density*self.g*self.h_vec)-(self.f_D*self.x_vec/self.dia*self.density/2*ss_v0**2) - self.set_initial_flow_velocity(ss_v0) - self.set_initial_pressure(ss_pressure) + # set the initial conditions + self.set_initial_flow_velocity(ss_v0,display_warning=False) + self.set_initial_pressure(ss_pressure,display_warning=False) -# getter +# getter - return attributes def get_info(self): new_line = '\n' angle_deg = round(self.angle/np.pi*180,3) @@ -139,6 +168,7 @@ class Druckrohrleitung_class: f"----------------------------- {new_line}" f"Length = {self.length:<10} {self.length_unit_disp} {new_line}" f"Diameter = {self.dia:<10} {self.length_unit_disp} {new_line}" + f"Hydraulic head = {self.head:<10} {self.length_unit_disp} {new_line}" f"Number of segments = {self.n_seg:<10} {new_line}" f"Number of nodes = {self.n_seg+1:<10} {new_line}" f"Length per segments = {self.dx:<10} {self.length_unit_disp} {new_line}" @@ -148,17 +178,16 @@ class Druckrohrleitung_class: f"Density of liquid = {self.density:<10} {self.density_unit_disp} {new_line}" f"Pressure wave vel. = {self.c:<10} {self.velocity_unit_disp} {new_line}" f"Simulation timestep = {self.dt:<10} {self.time_unit_disp} {new_line}" - f"Number of timesteps = {self.nt:<10} {new_line}" - f"Total simulation time = {self.nt*self.dt:<10} {self.time_unit_disp} {new_line}" f"----------------------------- {new_line}" f"Velocity and pressure distribution are vectors and are accessible by the .v and .p attribute of the pipeline object") print(print_str) - def get_current_pressure_distribution(self,disp=False): - if disp == True: + def get_current_pressure_distribution(self,disp_flag=False): + # disp_flag if one wants to directly plot the return of this method + if disp_flag == True: # convert to pressure unit disp return pressure_conversion(self.p,self.pressure_unit,self.pressure_unit_disp) - elif disp == False: + elif disp_flag == False: # stay in Pa return self.p def get_current_velocity_distribution(self): @@ -167,16 +196,16 @@ class Druckrohrleitung_class: def get_current_flux_distribution(self): return self.v*self.A - def get_lowest_pressure_per_node(self,disp=False): - if disp == True: + def get_lowest_pressure_per_node(self,disp_flag=False): + if disp_flag == True: # convert to pressure unit disp return pressure_conversion(self.p_min,self.pressure_unit,self.pressure_unit_disp) - elif disp == False: + elif disp_flag == False: # stay in Pa return self.p_min - def get_highest_pressure_per_node(self,disp=False): - if disp == True: + def get_highest_pressure_per_node(self,disp_flag=False): + if disp_flag == True: # convert to pressure unit disp return pressure_conversion(self.p_max,self.pressure_unit,self.pressure_unit_disp) - elif disp == False: + elif disp_flag == False: # stay in Pa return self.p_max def get_lowest_velocity_per_node(self): @@ -194,14 +223,15 @@ class Druckrohrleitung_class: def timestep_characteristic_method(self): # use the method of characteristics to calculate the pressure and velocities at all nodes except the boundary ones - # they are set with the .set_boundary_conditions_next_timestep() method beforehand + # they are set with the .set_boundary_conditions_next_timestep() method beforehand + # constants for cleaner formula nn = self.n_seg+1 # number of nodes rho = self.density # density of liquid c = self.c # pressure propagation velocity f_D = self.f_D # Darcy friction coefficient dt = self.dt # timestep - D = self.dia # pipeline diametet + D = self.dia # pipeline diameter g = self.g # graviational acceleration alpha = self.angle # pipeline angle diff --git a/Druckrohrleitung/Druckrohrleitung_test_steady_state.ipynb b/Druckrohrleitung/Druckrohrleitung_test_steady_state.ipynb index 1477a53..5e73d33 100644 --- a/Druckrohrleitung/Druckrohrleitung_test_steady_state.ipynb +++ b/Druckrohrleitung/Druckrohrleitung_test_steady_state.ipynb @@ -34,36 +34,32 @@ "pUnit_calc = 'Pa' # [text] DO NOT CHANGE! for pressure conversion in print statements and plot labels \n", "pUnit_conv = 'mWS' # [text] for pressure conversion in print statements and plot labels\n", "\n", - "\n", " # for Turbine\n", - "Tur_Q_nenn = 0.85 # [m³/s] nominal flux of turbine \n", - "Tur_p_nenn = pressure_conversion(10.6,'bar',pUnit_calc) # [Pa] nominal pressure of turbine \n", - "Tur_closingTime = 90. # [s] closing time of turbine\n", - "\n", + "Tur_Q_nenn = 0.85 # [m³/s] nominal flux of turbine \n", + "Tur_p_nenn = pressure_conversion(10.6,'bar',pUnit_calc) # [Pa] nominal pressure of turbine \n", + "Tur_closingTime = 90. # [s] closing time of turbine\n", "\n", " # for PI controller\n", "Con_targetLevel = 8. # [m]\n", "Con_K_p = 0.1 # [-] proportional constant of PI controller\n", - "Con_T_i = 10. # [s] timespan in which a steady state error is corrected by the intergal term\n", + "Con_T_i = 10. # [s] timespan in which a steady state error is corrected by the intergal term\n", "Con_deadbandRange = 0.05 # [m] Deadband range around targetLevel for which the controller does NOT intervene\n", "\n", - "\n", " # for pipeline\n", - "Pip_length = (535.+478.) # [m] length of pipeline\n", - "Pip_dia = 0.9 # [m] diameter of pipeline\n", - "Pip_area = Pip_dia**2/4*np.pi # [m²] crossectional area of pipeline\n", - "Pip_head = 105. # [m] hydraulic head of pipeline without reservoir\n", - "Pip_angle = np.arcsin(Pip_head/Pip_length) # [rad] elevation angle of pipeline \n", - "Pip_n_seg = 50 # [-] number of pipe segments in discretization\n", - "Pip_f_D = 0.014 # [-] Darcy friction factor\n", - "Pip_pw_vel = 500. # [m/s] propagation velocity of the pressure wave (pw) in the given pipeline\n", + "Pip_length = (535.+478.) # [m] length of pipeline\n", + "Pip_dia = 0.9 # [m] diameter of pipeline\n", + "Pip_area = Pip_dia**2/4*np.pi # [m²] crossectional area of pipeline\n", + "Pip_head = 105. # [m] hydraulic head of pipeline without reservoir\n", + "Pip_angle = np.arcsin(Pip_head/Pip_length) # [rad] elevation angle of pipeline \n", + "Pip_n_seg = 50 # [-] number of pipe segments in discretization\n", + "Pip_f_D = 0.014 # [-] Darcy friction factor\n", + "Pip_pw_vel = 500. # [m/s] propagation velocity of the pressure wave (pw) in the given pipeline\n", " # derivatives of the pipeline constants\n", - "Pip_dx = Pip_length/Pip_n_seg # [m] length of each pipe segment\n", - "Pip_dt = Pip_dx/Pip_pw_vel # [s] timestep according to method of characteristics\n", - "Pip_nn = Pip_n_seg+1 # [1] number of nodes\n", - "Pip_x_vec = np.arange(0,Pip_nn,1)*Pip_dx # [m] vector holding the distance of each node from the upstream reservoir along the pipeline\n", - "Pip_h_vec = np.arange(0,Pip_nn,1)*Pip_head/Pip_n_seg # [m] vector holding the vertival distance of each node from the upstream reservoir\n", - "\n", + "Pip_dx = Pip_length/Pip_n_seg # [m] length of each pipe segment\n", + "Pip_dt = Pip_dx/Pip_pw_vel # [s] timestep according to method of characteristics\n", + "Pip_nn = Pip_n_seg+1 # [1] number of nodes\n", + "Pip_x_vec = np.arange(0,Pip_nn,1)*Pip_dx # [m] vector holding the distance of each node from the upstream reservoir along the pipeline\n", + "Pip_h_vec = np.arange(0,Pip_nn,1)*Pip_head/Pip_n_seg # [m] vector holding the vertival distance of each node from the upstream reservoir\n", "\n", " # for reservoir\n", "Res_area_base = 74. # [m²] total base are of the cuboid reservoir \n", @@ -75,38 +71,68 @@ "Res_dt = Pip_dt/Res_nt # [s] harmonised timestep of reservoir time evolution\n", "\n", " # for general simulation\n", - "flux_init = Tur_Q_nenn/1.1 # [m³/s] initial flux through whole system for steady state initialization \n", - "level_init = Con_targetLevel # [m] initial water level in upstream reservoir for steady state initialization\n", - "simTime_target = 600. # [s] target for total simulation time (will vary slightly to fit with Pip_dt)\n", - "nt = int(simTime_target//Pip_dt) # [1] Number of timesteps of the whole system\n", - "t_vec = np.arange(0,nt+1,1)*Pip_dt # [s] time vector. At each step of t_vec the system parameters are stored\n" + "flux_init = Tur_Q_nenn/1.1 # [m³/s] initial flux through whole system for steady state initialization \n", + "level_init = Con_targetLevel # [m] initial water level in upstream reservoir for steady state initialization\n", + "simTime_target = 600. # [s] target for total simulation time (will vary slightly to fit with Pip_dt)\n", + "nt = int(simTime_target//Pip_dt) # [1] Number of timesteps of the whole system\n", + "t_vec = np.arange(0,nt+1,1)*Pip_dt # [s] time vector. At each step of t_vec the system parameters are stored\n" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The pipeline has the following attributes: \n", + "----------------------------- \n", + "Length = 1013.0 m \n", + "Diameter = 0.9 m \n", + "Hydraulic head = 105.0 m \n", + "Number of segments = 50 \n", + "Number of nodes = 51 \n", + "Length per segments = 20.26 m \n", + "Pipeline angle = 0.104 rad \n", + "Pipeline angle = 5.95° \n", + "Darcy friction factor = 0.014 \n", + "Density of liquid = 1000.0 kg/m³ \n", + "Pressure wave vel. = 500.0 m/s \n", + "Simulation timestep = 0.04052 s \n", + "----------------------------- \n", + "Velocity and pressure distribution are vectors and are accessible by the .v and .p attribute of the pipeline object\n", + "The pipeline has the following attributes: \n", + "----------------------------- \n", + "Length = 1013.0 m \n", + "Diameter = 0.9 m \n", + "Hydraulic head = 105.0 m \n", + "Number of segments = 50 \n", + "Number of nodes = 51 \n", + "Length per segments = 20.26 m \n", + "Pipeline angle = 0.104 rad \n", + "Pipeline angle = 5.95° \n", + "Darcy friction factor = 0.014 \n", + "Density of liquid = 1000.0 kg/m³ \n", + "Pressure wave vel. = 500.0 m/s \n", + "Simulation timestep = 0.04052 s \n", + "----------------------------- \n", + "Velocity and pressure distribution are vectors and are accessible by the .v and .p attribute of the pipeline object\n" + ] + } + ], "source": [ "# create objects\n", "\n", "# Upstream reservoir\n", - "reservoir = Ausgleichsbecken_class(Res_area_base,Res_area_out,Res_dt,Res_level_crit_lo,Res_level_crit_hi,rho)\n", + "reservoir = Ausgleichsbecken_class(Res_area_base,Res_area_out,Res_dt,pUnit_conv,Res_level_crit_lo,Res_level_crit_hi,rho)\n", "reservoir.set_steady_state(flux_init,level_init)\n", "\n", "# pipeline\n", - "pipe = Druckrohrleitung_class(Pip_length,Pip_dia,Pip_n_seg,Pip_angle,Pip_f_D,Pip_pw_vel,Pip_dt,pUnit_conv,rho)\n", - "pipe.set_steady_state(flux_init,level_init,Res_area_base,Pip_x_vec,Pip_h_vec)\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "reservoir.get_info(full=True)\n", - "pipe.get_info(full=True)" + "pipe = Druckrohrleitung_class(Pip_length,Pip_dia,Pip_head,Pip_n_seg,Pip_f_D,Pip_pw_vel,Pip_dt,pUnit_conv,rho)\n", + "pipe.set_steady_state(flux_init,reservoir.get_current_pressure())\n", + "pipe.get_info()\n" ] }, { @@ -171,35 +197,41 @@ "metadata": {}, "outputs": [ { - "ename": "KeyboardInterrupt", - "evalue": "", - "output_type": "error", - "traceback": [ - "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[1;31mKeyboardInterrupt\u001b[0m Traceback (most recent call last)", - "\u001b[1;32my:\\KELAG\\KS\\KS-PW\\04 Digitalisierung\\KSPWDEV Server\\Digital Trainee Projekt\\DT_Slot_3_Project_Repo\\Druckrohrleitung\\Druckrohrleitung_test_steady_state.ipynb Cell 7\u001b[0m in \u001b[0;36m\u001b[1;34m()\u001b[0m\n\u001b[0;32m 35\u001b[0m lo_01, \u001b[39m=\u001b[39m axs1[\u001b[39m1\u001b[39m]\u001b[39m.\u001b[39mplot(pl_vec,v_old,marker\u001b[39m=\u001b[39m\u001b[39m'\u001b[39m\u001b[39m.\u001b[39m\u001b[39m'\u001b[39m,c\u001b[39m=\u001b[39m\u001b[39m'\u001b[39m\u001b[39mblue\u001b[39m\u001b[39m'\u001b[39m)\n\u001b[0;32m 37\u001b[0m fig1\u001b[39m.\u001b[39msuptitle(\u001b[39mstr\u001b[39m(\u001b[39mround\u001b[39m(t_vec[it_pipe],\u001b[39m2\u001b[39m)) \u001b[39m+\u001b[39m \u001b[39m'\u001b[39m\u001b[39m/\u001b[39m\u001b[39m'\u001b[39m \u001b[39m+\u001b[39m \u001b[39mstr\u001b[39m(\u001b[39mround\u001b[39m(t_vec[\u001b[39m-\u001b[39m\u001b[39m1\u001b[39m],\u001b[39m2\u001b[39m)))\n\u001b[1;32m---> 38\u001b[0m fig1\u001b[39m.\u001b[39;49mcanvas\u001b[39m.\u001b[39;49mdraw()\n\u001b[0;32m 39\u001b[0m fig1\u001b[39m.\u001b[39mtight_layout()\n\u001b[0;32m 40\u001b[0m plt\u001b[39m.\u001b[39mpause(\u001b[39m0.000001\u001b[39m)\n", - "File \u001b[1;32mc:\\ProgramData\\Anaconda3\\envs\\Georg_DT_Slot3\\lib\\site-packages\\matplotlib\\backends\\backend_agg.py:436\u001b[0m, in \u001b[0;36mFigureCanvasAgg.draw\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m 432\u001b[0m \u001b[39m# Acquire a lock on the shared font cache.\u001b[39;00m\n\u001b[0;32m 433\u001b[0m \u001b[39mwith\u001b[39;00m RendererAgg\u001b[39m.\u001b[39mlock, \\\n\u001b[0;32m 434\u001b[0m (\u001b[39mself\u001b[39m\u001b[39m.\u001b[39mtoolbar\u001b[39m.\u001b[39m_wait_cursor_for_draw_cm() \u001b[39mif\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mtoolbar\n\u001b[0;32m 435\u001b[0m \u001b[39melse\u001b[39;00m nullcontext()):\n\u001b[1;32m--> 436\u001b[0m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mfigure\u001b[39m.\u001b[39;49mdraw(\u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mrenderer)\n\u001b[0;32m 437\u001b[0m \u001b[39m# A GUI class may be need to update a window using this draw, so\u001b[39;00m\n\u001b[0;32m 438\u001b[0m \u001b[39m# don't forget to call the superclass.\u001b[39;00m\n\u001b[0;32m 439\u001b[0m \u001b[39msuper\u001b[39m()\u001b[39m.\u001b[39mdraw()\n", - "File \u001b[1;32mc:\\ProgramData\\Anaconda3\\envs\\Georg_DT_Slot3\\lib\\site-packages\\matplotlib\\artist.py:73\u001b[0m, in \u001b[0;36m_finalize_rasterization..draw_wrapper\u001b[1;34m(artist, renderer, *args, **kwargs)\u001b[0m\n\u001b[0;32m 71\u001b[0m \u001b[39m@wraps\u001b[39m(draw)\n\u001b[0;32m 72\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mdraw_wrapper\u001b[39m(artist, renderer, \u001b[39m*\u001b[39margs, \u001b[39m*\u001b[39m\u001b[39m*\u001b[39mkwargs):\n\u001b[1;32m---> 73\u001b[0m result \u001b[39m=\u001b[39m draw(artist, renderer, \u001b[39m*\u001b[39;49margs, \u001b[39m*\u001b[39;49m\u001b[39m*\u001b[39;49mkwargs)\n\u001b[0;32m 74\u001b[0m \u001b[39mif\u001b[39;00m renderer\u001b[39m.\u001b[39m_rasterizing:\n\u001b[0;32m 75\u001b[0m renderer\u001b[39m.\u001b[39mstop_rasterizing()\n", - "File \u001b[1;32mc:\\ProgramData\\Anaconda3\\envs\\Georg_DT_Slot3\\lib\\site-packages\\matplotlib\\artist.py:50\u001b[0m, in \u001b[0;36mallow_rasterization..draw_wrapper\u001b[1;34m(artist, renderer)\u001b[0m\n\u001b[0;32m 47\u001b[0m \u001b[39mif\u001b[39;00m artist\u001b[39m.\u001b[39mget_agg_filter() \u001b[39mis\u001b[39;00m \u001b[39mnot\u001b[39;00m \u001b[39mNone\u001b[39;00m:\n\u001b[0;32m 48\u001b[0m renderer\u001b[39m.\u001b[39mstart_filter()\n\u001b[1;32m---> 50\u001b[0m \u001b[39mreturn\u001b[39;00m draw(artist, renderer)\n\u001b[0;32m 51\u001b[0m \u001b[39mfinally\u001b[39;00m:\n\u001b[0;32m 52\u001b[0m \u001b[39mif\u001b[39;00m artist\u001b[39m.\u001b[39mget_agg_filter() \u001b[39mis\u001b[39;00m \u001b[39mnot\u001b[39;00m \u001b[39mNone\u001b[39;00m:\n", - "File \u001b[1;32mc:\\ProgramData\\Anaconda3\\envs\\Georg_DT_Slot3\\lib\\site-packages\\matplotlib\\figure.py:2810\u001b[0m, in \u001b[0;36mFigure.draw\u001b[1;34m(self, renderer)\u001b[0m\n\u001b[0;32m 2807\u001b[0m \u001b[39m# ValueError can occur when resizing a window.\u001b[39;00m\n\u001b[0;32m 2809\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mpatch\u001b[39m.\u001b[39mdraw(renderer)\n\u001b[1;32m-> 2810\u001b[0m mimage\u001b[39m.\u001b[39;49m_draw_list_compositing_images(\n\u001b[0;32m 2811\u001b[0m renderer, \u001b[39mself\u001b[39;49m, artists, \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49msuppressComposite)\n\u001b[0;32m 2813\u001b[0m \u001b[39mfor\u001b[39;00m sfig \u001b[39min\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39msubfigs:\n\u001b[0;32m 2814\u001b[0m sfig\u001b[39m.\u001b[39mdraw(renderer)\n", - "File \u001b[1;32mc:\\ProgramData\\Anaconda3\\envs\\Georg_DT_Slot3\\lib\\site-packages\\matplotlib\\image.py:132\u001b[0m, in \u001b[0;36m_draw_list_compositing_images\u001b[1;34m(renderer, parent, artists, suppress_composite)\u001b[0m\n\u001b[0;32m 130\u001b[0m \u001b[39mif\u001b[39;00m not_composite \u001b[39mor\u001b[39;00m \u001b[39mnot\u001b[39;00m has_images:\n\u001b[0;32m 131\u001b[0m \u001b[39mfor\u001b[39;00m a \u001b[39min\u001b[39;00m artists:\n\u001b[1;32m--> 132\u001b[0m a\u001b[39m.\u001b[39;49mdraw(renderer)\n\u001b[0;32m 133\u001b[0m \u001b[39melse\u001b[39;00m:\n\u001b[0;32m 134\u001b[0m \u001b[39m# Composite any adjacent images together\u001b[39;00m\n\u001b[0;32m 135\u001b[0m image_group \u001b[39m=\u001b[39m []\n", - "File \u001b[1;32mc:\\ProgramData\\Anaconda3\\envs\\Georg_DT_Slot3\\lib\\site-packages\\matplotlib\\artist.py:50\u001b[0m, in \u001b[0;36mallow_rasterization..draw_wrapper\u001b[1;34m(artist, renderer)\u001b[0m\n\u001b[0;32m 47\u001b[0m \u001b[39mif\u001b[39;00m artist\u001b[39m.\u001b[39mget_agg_filter() \u001b[39mis\u001b[39;00m \u001b[39mnot\u001b[39;00m \u001b[39mNone\u001b[39;00m:\n\u001b[0;32m 48\u001b[0m renderer\u001b[39m.\u001b[39mstart_filter()\n\u001b[1;32m---> 50\u001b[0m \u001b[39mreturn\u001b[39;00m draw(artist, renderer)\n\u001b[0;32m 51\u001b[0m \u001b[39mfinally\u001b[39;00m:\n\u001b[0;32m 52\u001b[0m \u001b[39mif\u001b[39;00m artist\u001b[39m.\u001b[39mget_agg_filter() \u001b[39mis\u001b[39;00m \u001b[39mnot\u001b[39;00m \u001b[39mNone\u001b[39;00m:\n", - "File \u001b[1;32mc:\\ProgramData\\Anaconda3\\envs\\Georg_DT_Slot3\\lib\\site-packages\\matplotlib\\axes\\_base.py:3082\u001b[0m, in \u001b[0;36m_AxesBase.draw\u001b[1;34m(self, renderer)\u001b[0m\n\u001b[0;32m 3079\u001b[0m a\u001b[39m.\u001b[39mdraw(renderer)\n\u001b[0;32m 3080\u001b[0m renderer\u001b[39m.\u001b[39mstop_rasterizing()\n\u001b[1;32m-> 3082\u001b[0m mimage\u001b[39m.\u001b[39;49m_draw_list_compositing_images(\n\u001b[0;32m 3083\u001b[0m renderer, \u001b[39mself\u001b[39;49m, artists, \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mfigure\u001b[39m.\u001b[39;49msuppressComposite)\n\u001b[0;32m 3085\u001b[0m renderer\u001b[39m.\u001b[39mclose_group(\u001b[39m'\u001b[39m\u001b[39maxes\u001b[39m\u001b[39m'\u001b[39m)\n\u001b[0;32m 3086\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mstale \u001b[39m=\u001b[39m \u001b[39mFalse\u001b[39;00m\n", - "File \u001b[1;32mc:\\ProgramData\\Anaconda3\\envs\\Georg_DT_Slot3\\lib\\site-packages\\matplotlib\\image.py:132\u001b[0m, in \u001b[0;36m_draw_list_compositing_images\u001b[1;34m(renderer, parent, artists, suppress_composite)\u001b[0m\n\u001b[0;32m 130\u001b[0m \u001b[39mif\u001b[39;00m not_composite \u001b[39mor\u001b[39;00m \u001b[39mnot\u001b[39;00m has_images:\n\u001b[0;32m 131\u001b[0m \u001b[39mfor\u001b[39;00m a \u001b[39min\u001b[39;00m artists:\n\u001b[1;32m--> 132\u001b[0m a\u001b[39m.\u001b[39;49mdraw(renderer)\n\u001b[0;32m 133\u001b[0m \u001b[39melse\u001b[39;00m:\n\u001b[0;32m 134\u001b[0m \u001b[39m# Composite any adjacent images together\u001b[39;00m\n\u001b[0;32m 135\u001b[0m image_group \u001b[39m=\u001b[39m []\n", - "File \u001b[1;32mc:\\ProgramData\\Anaconda3\\envs\\Georg_DT_Slot3\\lib\\site-packages\\matplotlib\\artist.py:50\u001b[0m, in \u001b[0;36mallow_rasterization..draw_wrapper\u001b[1;34m(artist, renderer)\u001b[0m\n\u001b[0;32m 47\u001b[0m \u001b[39mif\u001b[39;00m artist\u001b[39m.\u001b[39mget_agg_filter() \u001b[39mis\u001b[39;00m \u001b[39mnot\u001b[39;00m \u001b[39mNone\u001b[39;00m:\n\u001b[0;32m 48\u001b[0m renderer\u001b[39m.\u001b[39mstart_filter()\n\u001b[1;32m---> 50\u001b[0m \u001b[39mreturn\u001b[39;00m draw(artist, renderer)\n\u001b[0;32m 51\u001b[0m \u001b[39mfinally\u001b[39;00m:\n\u001b[0;32m 52\u001b[0m \u001b[39mif\u001b[39;00m artist\u001b[39m.\u001b[39mget_agg_filter() \u001b[39mis\u001b[39;00m \u001b[39mnot\u001b[39;00m \u001b[39mNone\u001b[39;00m:\n", - "File \u001b[1;32mc:\\ProgramData\\Anaconda3\\envs\\Georg_DT_Slot3\\lib\\site-packages\\matplotlib\\axis.py:1170\u001b[0m, in \u001b[0;36mAxis.draw\u001b[1;34m(self, renderer, *args, **kwargs)\u001b[0m\n\u001b[0;32m 1163\u001b[0m tick\u001b[39m.\u001b[39mdraw(renderer)\n\u001b[0;32m 1165\u001b[0m \u001b[39m# scale up the axis label box to also find the neighbors, not\u001b[39;00m\n\u001b[0;32m 1166\u001b[0m \u001b[39m# just the tick labels that actually overlap note we need a\u001b[39;00m\n\u001b[0;32m 1167\u001b[0m \u001b[39m# *copy* of the axis label box because we don't want to scale\u001b[39;00m\n\u001b[0;32m 1168\u001b[0m \u001b[39m# the actual bbox\u001b[39;00m\n\u001b[1;32m-> 1170\u001b[0m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49m_update_label_position(renderer)\n\u001b[0;32m 1172\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mlabel\u001b[39m.\u001b[39mdraw(renderer)\n\u001b[0;32m 1174\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_update_offset_text_position(ticklabelBoxes, ticklabelBoxes2)\n", - "File \u001b[1;32mc:\\ProgramData\\Anaconda3\\envs\\Georg_DT_Slot3\\lib\\site-packages\\matplotlib\\axis.py:2352\u001b[0m, in \u001b[0;36mYAxis._update_label_position\u001b[1;34m(self, renderer)\u001b[0m\n\u001b[0;32m 2348\u001b[0m \u001b[39mreturn\u001b[39;00m\n\u001b[0;32m 2350\u001b[0m \u001b[39m# get bounding boxes for this axis and any siblings\u001b[39;00m\n\u001b[0;32m 2351\u001b[0m \u001b[39m# that have been set by `fig.align_ylabels()`\u001b[39;00m\n\u001b[1;32m-> 2352\u001b[0m bboxes, bboxes2 \u001b[39m=\u001b[39m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49m_get_tick_boxes_siblings(renderer\u001b[39m=\u001b[39;49mrenderer)\n\u001b[0;32m 2354\u001b[0m x, y \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mlabel\u001b[39m.\u001b[39mget_position()\n\u001b[0;32m 2355\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mlabel_position \u001b[39m==\u001b[39m \u001b[39m'\u001b[39m\u001b[39mleft\u001b[39m\u001b[39m'\u001b[39m:\n", - "File \u001b[1;32mc:\\ProgramData\\Anaconda3\\envs\\Georg_DT_Slot3\\lib\\site-packages\\matplotlib\\axis.py:1880\u001b[0m, in \u001b[0;36mAxis._get_tick_boxes_siblings\u001b[1;34m(self, renderer)\u001b[0m\n\u001b[0;32m 1878\u001b[0m \u001b[39mfor\u001b[39;00m ax \u001b[39min\u001b[39;00m grouper\u001b[39m.\u001b[39mget_siblings(\u001b[39mself\u001b[39m\u001b[39m.\u001b[39maxes):\n\u001b[0;32m 1879\u001b[0m axis \u001b[39m=\u001b[39m \u001b[39mgetattr\u001b[39m(ax, \u001b[39mf\u001b[39m\u001b[39m\"\u001b[39m\u001b[39m{\u001b[39;00maxis_name\u001b[39m}\u001b[39;00m\u001b[39maxis\u001b[39m\u001b[39m\"\u001b[39m)\n\u001b[1;32m-> 1880\u001b[0m ticks_to_draw \u001b[39m=\u001b[39m axis\u001b[39m.\u001b[39;49m_update_ticks()\n\u001b[0;32m 1881\u001b[0m tlb, tlb2 \u001b[39m=\u001b[39m axis\u001b[39m.\u001b[39m_get_tick_bboxes(ticks_to_draw, renderer)\n\u001b[0;32m 1882\u001b[0m bboxes\u001b[39m.\u001b[39mextend(tlb)\n", - "File \u001b[1;32mc:\\ProgramData\\Anaconda3\\envs\\Georg_DT_Slot3\\lib\\site-packages\\matplotlib\\axis.py:1053\u001b[0m, in \u001b[0;36mAxis._update_ticks\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m 1051\u001b[0m tick\u001b[39m.\u001b[39mset_label1(label)\n\u001b[0;32m 1052\u001b[0m tick\u001b[39m.\u001b[39mset_label2(label)\n\u001b[1;32m-> 1053\u001b[0m minor_locs \u001b[39m=\u001b[39m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mget_minorticklocs()\n\u001b[0;32m 1054\u001b[0m minor_labels \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mminor\u001b[39m.\u001b[39mformatter\u001b[39m.\u001b[39mformat_ticks(minor_locs)\n\u001b[0;32m 1055\u001b[0m minor_ticks \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mget_minor_ticks(\u001b[39mlen\u001b[39m(minor_locs))\n", - "File \u001b[1;32mc:\\ProgramData\\Anaconda3\\envs\\Georg_DT_Slot3\\lib\\site-packages\\matplotlib\\axis.py:1282\u001b[0m, in \u001b[0;36mAxis.get_minorticklocs\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m 1280\u001b[0m \u001b[39m\"\"\"Return this Axis' minor tick locations in data coordinates.\"\"\"\u001b[39;00m\n\u001b[0;32m 1281\u001b[0m \u001b[39m# Remove minor ticks duplicating major ticks.\u001b[39;00m\n\u001b[1;32m-> 1282\u001b[0m major_locs \u001b[39m=\u001b[39m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mmajor\u001b[39m.\u001b[39;49mlocator()\n\u001b[0;32m 1283\u001b[0m minor_locs \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mminor\u001b[39m.\u001b[39mlocator()\n\u001b[0;32m 1284\u001b[0m transform \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_scale\u001b[39m.\u001b[39mget_transform()\n", - "File \u001b[1;32mc:\\ProgramData\\Anaconda3\\envs\\Georg_DT_Slot3\\lib\\site-packages\\matplotlib\\ticker.py:2114\u001b[0m, in \u001b[0;36mMaxNLocator.__call__\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m 2112\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39m__call__\u001b[39m(\u001b[39mself\u001b[39m):\n\u001b[0;32m 2113\u001b[0m vmin, vmax \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39maxis\u001b[39m.\u001b[39mget_view_interval()\n\u001b[1;32m-> 2114\u001b[0m \u001b[39mreturn\u001b[39;00m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mtick_values(vmin, vmax)\n", - "File \u001b[1;32mc:\\ProgramData\\Anaconda3\\envs\\Georg_DT_Slot3\\lib\\site-packages\\matplotlib\\ticker.py:2122\u001b[0m, in \u001b[0;36mMaxNLocator.tick_values\u001b[1;34m(self, vmin, vmax)\u001b[0m\n\u001b[0;32m 2119\u001b[0m vmin \u001b[39m=\u001b[39m \u001b[39m-\u001b[39mvmax\n\u001b[0;32m 2120\u001b[0m vmin, vmax \u001b[39m=\u001b[39m mtransforms\u001b[39m.\u001b[39mnonsingular(\n\u001b[0;32m 2121\u001b[0m vmin, vmax, expander\u001b[39m=\u001b[39m\u001b[39m1e-13\u001b[39m, tiny\u001b[39m=\u001b[39m\u001b[39m1e-14\u001b[39m)\n\u001b[1;32m-> 2122\u001b[0m locs \u001b[39m=\u001b[39m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49m_raw_ticks(vmin, vmax)\n\u001b[0;32m 2124\u001b[0m prune \u001b[39m=\u001b[39m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_prune\n\u001b[0;32m 2125\u001b[0m \u001b[39mif\u001b[39;00m prune \u001b[39m==\u001b[39m \u001b[39m'\u001b[39m\u001b[39mlower\u001b[39m\u001b[39m'\u001b[39m:\n", - "File \u001b[1;32mc:\\ProgramData\\Anaconda3\\envs\\Georg_DT_Slot3\\lib\\site-packages\\matplotlib\\ticker.py:2061\u001b[0m, in \u001b[0;36mMaxNLocator._raw_ticks\u001b[1;34m(self, vmin, vmax)\u001b[0m\n\u001b[0;32m 2059\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_nbins \u001b[39m==\u001b[39m \u001b[39m'\u001b[39m\u001b[39mauto\u001b[39m\u001b[39m'\u001b[39m:\n\u001b[0;32m 2060\u001b[0m \u001b[39mif\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39maxis \u001b[39mis\u001b[39;00m \u001b[39mnot\u001b[39;00m \u001b[39mNone\u001b[39;00m:\n\u001b[1;32m-> 2061\u001b[0m nbins \u001b[39m=\u001b[39m np\u001b[39m.\u001b[39mclip(\u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49maxis\u001b[39m.\u001b[39;49mget_tick_space(),\n\u001b[0;32m 2062\u001b[0m \u001b[39mmax\u001b[39m(\u001b[39m1\u001b[39m, \u001b[39mself\u001b[39m\u001b[39m.\u001b[39m_min_n_ticks \u001b[39m-\u001b[39m \u001b[39m1\u001b[39m), \u001b[39m9\u001b[39m)\n\u001b[0;32m 2063\u001b[0m \u001b[39melse\u001b[39;00m:\n\u001b[0;32m 2064\u001b[0m nbins \u001b[39m=\u001b[39m \u001b[39m9\u001b[39m\n", - "File \u001b[1;32mc:\\ProgramData\\Anaconda3\\envs\\Georg_DT_Slot3\\lib\\site-packages\\matplotlib\\axis.py:2524\u001b[0m, in \u001b[0;36mYAxis.get_tick_space\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m 2523\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mget_tick_space\u001b[39m(\u001b[39mself\u001b[39m):\n\u001b[1;32m-> 2524\u001b[0m ends \u001b[39m=\u001b[39m mtransforms\u001b[39m.\u001b[39;49mBbox\u001b[39m.\u001b[39;49mfrom_bounds(\u001b[39m0\u001b[39;49m, \u001b[39m0\u001b[39;49m, \u001b[39m1\u001b[39;49m, \u001b[39m1\u001b[39;49m)\n\u001b[0;32m 2525\u001b[0m ends \u001b[39m=\u001b[39m ends\u001b[39m.\u001b[39mtransformed(\u001b[39mself\u001b[39m\u001b[39m.\u001b[39maxes\u001b[39m.\u001b[39mtransAxes \u001b[39m-\u001b[39m\n\u001b[0;32m 2526\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mfigure\u001b[39m.\u001b[39mdpi_scale_trans)\n\u001b[0;32m 2527\u001b[0m length \u001b[39m=\u001b[39m ends\u001b[39m.\u001b[39mheight \u001b[39m*\u001b[39m \u001b[39m72\u001b[39m\n", - "File \u001b[1;32mc:\\ProgramData\\Anaconda3\\envs\\Georg_DT_Slot3\\lib\\site-packages\\matplotlib\\transforms.py:820\u001b[0m, in \u001b[0;36mBbox.from_bounds\u001b[1;34m(x0, y0, width, height)\u001b[0m\n\u001b[0;32m 813\u001b[0m \u001b[39m@staticmethod\u001b[39m\n\u001b[0;32m 814\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mfrom_bounds\u001b[39m(x0, y0, width, height):\n\u001b[0;32m 815\u001b[0m \u001b[39m\"\"\"\u001b[39;00m\n\u001b[0;32m 816\u001b[0m \u001b[39m Create a new `Bbox` from *x0*, *y0*, *width* and *height*.\u001b[39;00m\n\u001b[0;32m 817\u001b[0m \n\u001b[0;32m 818\u001b[0m \u001b[39m *width* and *height* may be negative.\u001b[39;00m\n\u001b[0;32m 819\u001b[0m \u001b[39m \"\"\"\u001b[39;00m\n\u001b[1;32m--> 820\u001b[0m \u001b[39mreturn\u001b[39;00m Bbox\u001b[39m.\u001b[39;49mfrom_extents(x0, y0, x0 \u001b[39m+\u001b[39;49m width, y0 \u001b[39m+\u001b[39;49m height)\n", - "File \u001b[1;32mc:\\ProgramData\\Anaconda3\\envs\\Georg_DT_Slot3\\lib\\site-packages\\matplotlib\\transforms.py:839\u001b[0m, in \u001b[0;36mBbox.from_extents\u001b[1;34m(minpos, *args)\u001b[0m\n\u001b[0;32m 822\u001b[0m \u001b[39m@staticmethod\u001b[39m\n\u001b[0;32m 823\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mfrom_extents\u001b[39m(\u001b[39m*\u001b[39margs, minpos\u001b[39m=\u001b[39m\u001b[39mNone\u001b[39;00m):\n\u001b[0;32m 824\u001b[0m \u001b[39m\"\"\"\u001b[39;00m\n\u001b[0;32m 825\u001b[0m \u001b[39m Create a new Bbox from *left*, *bottom*, *right* and *top*.\u001b[39;00m\n\u001b[0;32m 826\u001b[0m \n\u001b[1;32m (...)\u001b[0m\n\u001b[0;32m 837\u001b[0m \u001b[39m scales where negative bounds result in floating point errors.\u001b[39;00m\n\u001b[0;32m 838\u001b[0m \u001b[39m \"\"\"\u001b[39;00m\n\u001b[1;32m--> 839\u001b[0m bbox \u001b[39m=\u001b[39m Bbox(np\u001b[39m.\u001b[39;49mreshape(args, (\u001b[39m2\u001b[39;49m, \u001b[39m2\u001b[39;49m)))\n\u001b[0;32m 840\u001b[0m \u001b[39mif\u001b[39;00m minpos \u001b[39mis\u001b[39;00m \u001b[39mnot\u001b[39;00m \u001b[39mNone\u001b[39;00m:\n\u001b[0;32m 841\u001b[0m bbox\u001b[39m.\u001b[39m_minpos[:] \u001b[39m=\u001b[39m minpos\n", - "File \u001b[1;32mc:\\ProgramData\\Anaconda3\\envs\\Georg_DT_Slot3\\lib\\site-packages\\matplotlib\\transforms.py:775\u001b[0m, in \u001b[0;36mBbox.__init__\u001b[1;34m(self, points, **kwargs)\u001b[0m\n\u001b[0;32m 768\u001b[0m \u001b[39m\"\"\"\u001b[39;00m\n\u001b[0;32m 769\u001b[0m \u001b[39mParameters\u001b[39;00m\n\u001b[0;32m 770\u001b[0m \u001b[39m----------\u001b[39;00m\n\u001b[0;32m 771\u001b[0m \u001b[39mpoints : ndarray\u001b[39;00m\n\u001b[0;32m 772\u001b[0m \u001b[39m A 2x2 numpy array of the form ``[[x0, y0], [x1, y1]]``.\u001b[39;00m\n\u001b[0;32m 773\u001b[0m \u001b[39m\"\"\"\u001b[39;00m\n\u001b[0;32m 774\u001b[0m \u001b[39msuper\u001b[39m()\u001b[39m.\u001b[39m\u001b[39m__init__\u001b[39m(\u001b[39m*\u001b[39m\u001b[39m*\u001b[39mkwargs)\n\u001b[1;32m--> 775\u001b[0m points \u001b[39m=\u001b[39m np\u001b[39m.\u001b[39;49masarray(points, \u001b[39mfloat\u001b[39;49m)\n\u001b[0;32m 776\u001b[0m \u001b[39mif\u001b[39;00m points\u001b[39m.\u001b[39mshape \u001b[39m!=\u001b[39m (\u001b[39m2\u001b[39m, \u001b[39m2\u001b[39m):\n\u001b[0;32m 777\u001b[0m \u001b[39mraise\u001b[39;00m \u001b[39mValueError\u001b[39;00m(\u001b[39m'\u001b[39m\u001b[39mBbox points must be of the form \u001b[39m\u001b[39m'\u001b[39m\n\u001b[0;32m 778\u001b[0m \u001b[39m'\u001b[39m\u001b[39m\"\u001b[39m\u001b[39m[[x0, y0], [x1, y1]]\u001b[39m\u001b[39m\"\u001b[39m\u001b[39m.\u001b[39m\u001b[39m'\u001b[39m)\n", - "\u001b[1;31mKeyboardInterrupt\u001b[0m: " + "name": "stdout", + "output_type": "stream", + "text": [ + "The cuboid reservoir has the following attributes: \n", + "----------------------------- \n", + "Base area = 74.0 m² \n", + "Outflux area = 0.636 m² \n", + "Current level = 8.0 m\n", + "Critical level low = 0.0 m \n", + "Critical level high = inf m \n", + "Volume in reservoir = 592.0 m³ \n", + "Current influx = 0.773 m³/s \n", + "Current outflux = 0.773 m³/s \n", + "Current outflux vel = 1.215 m/s \n", + "Current pipe pressure = 7.854 mWS \n", + "Simulation timestep = 0.001013 s \n", + "Density of liquid = 1000.0 kg/m³ \n", + "----------------------------- \n", + "\n", + "The pipeline has the following attributes: \n", + "----------------------------- \n", + "Length = 1013.0 m \n", + "Diameter = 0.9 m \n", + "Hydraulic head = 105.0 m \n", + "Number of segments = 50 \n", + "Number of nodes = 51 \n", + "Length per segments = 20.26 m \n", + "Pipeline angle = 0.104 rad \n", + "Pipeline angle = 5.95° \n", + "Darcy friction factor = 0.014 \n", + "Density of liquid = 1000.0 kg/m³ \n", + "Pressure wave vel. = 500.0 m/s \n", + "Simulation timestep = 0.04052 s \n", + "----------------------------- \n", + "Velocity and pressure distribution are vectors and are accessible by the .v and .p attribute of the pipeline object\n" ] } ], @@ -246,12 +278,12 @@ " plt.pause(0.000001)\n", "\n", "reservoir.get_info(full=True)\n", - "pipe.get_info(full=True)" + "pipe.get_info()" ] }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 7, "metadata": {}, "outputs": [], "source": [ diff --git a/Regler/Pegelregler_test.ipynb b/Regler/Pegelregler_test.ipynb index 4d16b95..c913b6b 100644 --- a/Regler/Pegelregler_test.ipynb +++ b/Regler/Pegelregler_test.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 27, + "execution_count": 1, "metadata": {}, "outputs": [], "source": [ @@ -23,7 +23,7 @@ }, { "cell_type": "code", - "execution_count": 28, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ @@ -85,7 +85,7 @@ }, { "cell_type": "code", - "execution_count": 29, + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ @@ -93,7 +93,7 @@ "offset_pressure = pressure_conversion(Pip_head,'mws',pUnit_calc)\n", "\n", "# Upstream reservoir\n", - "reservoir = Ausgleichsbecken_class(Res_area_base,Res_area_out,Res_dt,Res_level_crit_lo,Res_level_crit_hi,rho)\n", + "reservoir = Ausgleichsbecken_class(Res_area_base,Res_area_out,Res_dt,pUnit_conv,Res_level_crit_lo,Res_level_crit_hi,rho)\n", "reservoir.set_steady_state(flux_init,level_init)\n", "\n", "# downstream turbine\n", @@ -108,7 +108,7 @@ }, { "cell_type": "code", - "execution_count": 30, + "execution_count": 4, "metadata": {}, "outputs": [], "source": [ @@ -124,7 +124,7 @@ }, { "cell_type": "code", - "execution_count": 31, + "execution_count": 5, "metadata": {}, "outputs": [ { @@ -181,7 +181,7 @@ }, { "cell_type": "code", - "execution_count": 32, + "execution_count": 6, "metadata": {}, "outputs": [], "source": [ diff --git a/Regler/Regler_class_file.py b/Regler/Regler_class_file.py index 88c298b..5b9135b 100644 --- a/Regler/Regler_class_file.py +++ b/Regler/Regler_class_file.py @@ -1,6 +1,7 @@ import numpy as np #based on https://en.wikipedia.org/wiki/PID_controller#Discrete_implementation +# performance parameters for controllers def trap_int(vec,timestep): # numerical integration via the trapeziod rule to calculate the performance parameters l = np.size(vec) @@ -41,6 +42,7 @@ def ITAE_fun(error_history,timestep): itae = trap_int(np.abs(e),dt) return itae +# P controller class P_controller_class: # def __init__(self,setpoint,proportionality_constant): # self.SP = setpoint @@ -72,14 +74,14 @@ class P_controller_class: def __init__(self): pass - +# PI controller class PI_controller_class: # 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 # integration time + self.Ti = Ti # ~integration time self.dt = timestep # use a list to be able to append more easily - will get converted to np.array when needed self.error_history = [0] @@ -88,15 +90,13 @@ class PI_controller_class: self.cv_upper_limit = upper_limit # limits for the controll variable # setter - - def set_setpoint(self,setpoint): self.SP = setpoint def set_control_variable(self,control_variable, display_warning=True): if display_warning == True: - print('WARNING! You are setting the control variable of the PI controller manually \ - and are not using the .update_controll_variable() method') + print('WARNING! You are setting the control variable of the PI controller manually! \ + Consider using the .update_controll_variable() method instead.') self.control_variable = control_variable # getter @@ -167,15 +167,14 @@ class PI_controller_class: # only if that is the case, change control variable if abs(self.error) > self.db: new_control = cv+Kp*(e0-e1)+dt/Ti*e0 + # 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 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 21f895e..d8fbf41 100644 --- a/Turbinen/Turbinen_class_file.py +++ b/Turbinen/Turbinen_class_file.py @@ -11,8 +11,8 @@ sys.path.append(parent) from functions.pressure_conversion import pressure_conversion class Francis_Turbine: - # units - # make sure that units and print units are the same +# units + # make sure that units and display units are the same # units are used to label graphs and disp 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}$' @@ -32,7 +32,16 @@ class Francis_Turbine: g = 9.81 # m/s² gravitational acceleration # init - def __init__(self, Q_nenn,p_nenn,t_closing,timestep,pressure_unit_disp): + def __init__(self,Q_nenn,p_nenn,t_closing,timestep,pressure_unit_disp): + """ + Creates a turbine with given attributes in this order: \n + Nominal flux [m³/s] \n + Nominal pressure [Pa] \n + Closing time [s] \n + Simulation timestep [s] \n + Pressure unit for displaying [string] \n + + """ self.Q_n = Q_nenn # nominal flux self.p_n = p_nenn # nominal pressure self.LA_n = 1. # 100% # nominal Leitapparatöffnung @@ -42,21 +51,21 @@ class Francis_Turbine: self.pressure_unit_disp = pressure_unit_disp - # initialize for get_info() - parameters will be converted to display -1 if not overwritten - self.p = pressure_conversion(-1,self.pressure_unit_disp,self.pressure_unit) - self.Q = -1. - self.LA = -0.01 + # initialize for get_info() + self.p = -np.inf + self.Q = -np.inf + self.LA = -np.inf -# setter +# setter - set attributes 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('You are setting the guide vane opening of the turbine manually. \n \ This is not an intended use of this method. \n \ Refer to the .update_LA() method instead.') + # set Leitapparatöffnung + self.LA = LA def set_pressure(self,pressure): # set pressure in front of the turbine @@ -70,7 +79,7 @@ class Francis_Turbine: raise Exception('LA out of range [0;1]') self.set_LA(ss_LA,display_warning=False) -#getter +#getter - get attributes 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 @@ -83,8 +92,11 @@ class Francis_Turbine: def get_current_LA(self): return self.LA - def get_current_pressure(self): - return pressure_conversion(self.p,self.pressure_unit,self.pressure_unit_disp) + def get_current_pressure(self,disp_flag=True): + if disp_flag == True: + return pressure_conversion(self.p,self.pressure_unit,self.pressure_unit_disp) + else: + return self.p def get_info(self, full = False): new_line = '\n' @@ -135,7 +147,9 @@ class Francis_Turbine: # methods def converge(self,convergence_parameters): # small numerical disturbances (~1e-12 m/s) in the velocity can get amplified at the turbine node, because the new velocity of the turbine and the - # new pressure from the forward characteristic are not compatible. + # new pressure from the forward characteristic are not perfectly compatible. + # Therefore, iterate the flux and the pressure so long, until they converge + eps = 1e-12 # convergence criterion: iteration change < eps iteration_change = 1. # change in Q from one iteration to the next i = 0 # safety variable. break loop if it exceeds 1e6 iterations @@ -150,20 +164,18 @@ class Francis_Turbine: rho = convergence_parameters[7] # density of the liquid dt = convergence_parameters[8] # timestep of the characteristic method - p_old = self.get_current_pressure() Q_old = self.get_current_Q() v_old = Q_old/area_pipe while iteration_change > eps: - self.set_pressure(p_old) + p_new = p-rho*c*(v_old-v)+rho*c*dt*g*np.sin(alpha)-f_D*rho*c*dt/(2*D)*abs(v)*v + self.set_pressure(p_new) Q_new = self.get_current_Q() v_new = Q_new/area_pipe - p_new = p-rho*c*(v_old-v)+rho*c*dt*g*np.sin(alpha)-f_D*rho*c*dt/(2*D)*abs(v)*v iteration_change = abs(Q_old-Q_new) Q_old = Q_new.copy() - p_old = p_new.copy() v_old = v_new.copy() i = i+1 if i == 1e6: diff --git a/Turbinen/Turbinen_test_steady_state.ipynb b/Turbinen/Turbinen_test_steady_state.ipynb index 4037691..2891587 100644 --- a/Turbinen/Turbinen_test_steady_state.ipynb +++ b/Turbinen/Turbinen_test_steady_state.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 8, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -23,7 +23,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -35,12 +35,10 @@ "pUnit_calc = 'Pa' # [text] DO NOT CHANGE! for pressure conversion in print statements and plot labels \n", "pUnit_conv = 'mWS' # [text] for pressure conversion in print statements and plot labels\n", "\n", - "\n", " # for Turbine\n", - "Tur_Q_nenn = 0.85 # [m³/s] nominal flux of turbine \n", - "Tur_p_nenn = pressure_conversion(10.6,'bar',pUnit_calc) # [Pa] nominal pressure of turbine \n", - "Tur_closingTime = 90. # [s] closing time of turbine\n", - "\n", + "Tur_Q_nenn = 0.85 # [m³/s] nominal flux of turbine \n", + "Tur_p_nenn = pressure_conversion(10.6,'bar',pUnit_calc) # [Pa] nominal pressure of turbine \n", + "Tur_closingTime = 90. # [s] closing time of turbine\n", "\n", " # for PI controller\n", "Con_targetLevel = 8. # [m]\n", @@ -48,23 +46,21 @@ "Con_T_i = 10. # [s] timespan in which a steady state error is corrected by the intergal term\n", "Con_deadbandRange = 0.05 # [m] Deadband range around targetLevel for which the controller does NOT intervene\n", "\n", - "\n", " # for pipeline\n", - "Pip_length = (535.+478.) # [m] length of pipeline\n", - "Pip_dia = 0.9 # [m] diameter of pipeline\n", - "Pip_area = Pip_dia**2/4*np.pi # [m²] crossectional area of pipeline\n", - "Pip_head = 105. # [m] hydraulic head of pipeline without reservoir\n", - "Pip_angle = np.arcsin(Pip_head/Pip_length) # [rad] elevation angle of pipeline \n", - "Pip_n_seg = 50 # [-] number of pipe segments in discretization\n", - "Pip_f_D = 0.014 # [-] Darcy friction factor\n", - "Pip_pw_vel = 500. # [m/s] propagation velocity of the pressure wave (pw) in the given pipeline\n", + "Pip_length = (535.+478.) # [m] length of pipeline\n", + "Pip_dia = 0.9 # [m] diameter of pipeline\n", + "Pip_area = Pip_dia**2/4*np.pi # [m²] crossectional area of pipeline\n", + "Pip_head = 105. # [m] hydraulic head of pipeline without reservoir\n", + "Pip_angle = np.arcsin(Pip_head/Pip_length) # [rad] elevation angle of pipeline \n", + "Pip_n_seg = 50 # [-] number of pipe segments in discretization\n", + "Pip_f_D = 0.014 # [-] Darcy friction factor\n", + "Pip_pw_vel = 500. # [m/s] propagation velocity of the pressure wave (pw) in the given pipeline\n", " # derivatives of the pipeline constants\n", - "Pip_dx = Pip_length/Pip_n_seg # [m] length of each pipe segment\n", - "Pip_dt = Pip_dx/Pip_pw_vel # [s] timestep according to method of characteristics\n", - "Pip_nn = Pip_n_seg+1 # [1] number of nodes\n", - "Pip_x_vec = np.arange(0,Pip_nn,1)*Pip_dx # [m] vector holding the distance of each node from the upstream reservoir along the pipeline\n", - "Pip_h_vec = np.arange(0,Pip_nn,1)*Pip_head/Pip_n_seg # [m] vector holding the vertival distance of each node from the upstream reservoir\n", - "\n", + "Pip_dx = Pip_length/Pip_n_seg # [m] length of each pipe segment\n", + "Pip_dt = Pip_dx/Pip_pw_vel # [s] timestep according to method of characteristics\n", + "Pip_nn = Pip_n_seg+1 # [1] number of nodes\n", + "Pip_x_vec = np.arange(0,Pip_nn,1)*Pip_dx # [m] vector holding the distance of each node from the upstream reservoir along the pipeline\n", + "Pip_h_vec = np.arange(0,Pip_nn,1)*Pip_head/Pip_n_seg # [m] vector holding the vertival distance of each node from the upstream reservoir\n", "\n", " # for reservoir\n", "Res_area_base = 74. # [m²] total base are of the cuboid reservoir \n", @@ -76,16 +72,16 @@ "Res_dt = Pip_dt/Res_nt # [s] harmonised timestep of reservoir time evolution\n", "\n", " # for general simulation\n", - "flux_init = Tur_Q_nenn/1.1 # [m³/s] initial flux through whole system for steady state initialization \n", - "level_init = Con_targetLevel # [m] initial water level in upstream reservoir for steady state initialization\n", - "simTime_target = 600. # [s] target for total simulation time (will vary slightly to fit with Pip_dt)\n", - "nt = int(simTime_target//Pip_dt) # [1] Number of timesteps of the whole system\n", - "t_vec = np.arange(0,nt+1,1)*Pip_dt # [s] time vector. At each step of t_vec the system parameters are stored\n" + "flux_init = Tur_Q_nenn/1.1 # [m³/s] initial flux through whole system for steady state initialization \n", + "level_init = Con_targetLevel # [m] initial water level in upstream reservoir for steady state initialization\n", + "simTime_target = 100. # [s] target for total simulation time (will vary slightly to fit with Pip_dt)\n", + "nt = int(simTime_target//Pip_dt) # [1] Number of timesteps of the whole system\n", + "t_vec = np.arange(0,nt+1,1)*Pip_dt # [s] time vector. At each step of t_vec the system parameters are stored\n" ] }, { "cell_type": "code", - "execution_count": 10, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -96,8 +92,8 @@ "reservoir.set_steady_state(flux_init,level_init)\n", "\n", "# pipeline\n", - "pipe = Druckrohrleitung_class(Pip_length,Pip_dia,Pip_n_seg,Pip_angle,Pip_f_D,Pip_pw_vel,Pip_dt,pUnit_conv,rho)\n", - "pipe.set_steady_state(flux_init,level_init,Res_area_base,Pip_x_vec,Pip_h_vec)\n", + "pipe = Druckrohrleitung_class(Pip_length,Pip_dia,Pip_head,Pip_n_seg,Pip_f_D,Pip_pw_vel,Pip_dt,pUnit_conv,rho)\n", + "pipe.set_steady_state(flux_init,reservoir.get_current_pressure())\n", "\n", "# downstream turbine\n", "turbine = Francis_Turbine(Tur_Q_nenn,Tur_p_nenn,Tur_closingTime,Pip_dt,pUnit_conv)\n", @@ -114,7 +110,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -133,22 +129,22 @@ "Q_in_vec = np.zeros_like(t_vec)\n", "Q_in_vec[0] = flux_init\n", "\n", - "v_boundary_res = np.zeros_like(t_vec)\n", - "v_boundary_tur = np.zeros_like(t_vec)\n", - "Q_boundary_res = np.zeros_like(t_vec)\n", - "Q_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", + "v_boundary_res = np.zeros_like(t_vec)\n", + "v_boundary_tur = np.zeros_like(t_vec)\n", + "Q_boundary_res = np.zeros_like(t_vec)\n", + "Q_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", - "level_vec = np.full_like(t_vec,level_init) # level at the end of each pipeline timestep\n", - "volume_vec = np.full_like(t_vec,reservoir.get_current_volume()) # volume at the end of each pipeline timestep\n", + "level_vec = np.full_like(t_vec,level_init) # level at the end of each pipeline timestep\n", + "volume_vec = np.full_like(t_vec,reservoir.get_current_volume()) # volume at the end of each pipeline timestep\n", "\n", - "v_boundary_res[0] = v_old[0]\n", - "v_boundary_tur[0] = v_old[-1] \n", - "Q_boundary_res[0] = Q_old[0]\n", - "Q_boundary_tur[0] = Q_old[-1]\n", - "p_boundary_res[0] = p_old[0]\n", - "p_boundary_tur[0] = p_old[-1]\n", + "v_boundary_res[0] = v_old[0]\n", + "v_boundary_tur[0] = v_old[-1] \n", + "Q_boundary_res[0] = Q_old[0]\n", + "Q_boundary_tur[0] = Q_old[-1]\n", + "p_boundary_res[0] = p_old[0]\n", + "p_boundary_tur[0] = p_old[-1]\n", "\n", "LA_soll_vec = np.full_like(t_vec,turbine.get_current_LA())\n", "LA_ist_vec = np.full_like(t_vec,turbine.get_current_LA())\n", @@ -158,7 +154,7 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -176,8 +172,8 @@ "axs1[1].set_ylabel(r'$Q$ [$\\mathrm{m}^3 / \\mathrm{s}$]')\n", "lo_p, = axs1[0].plot(Pip_x_vec,pressure_conversion(p_old,pUnit_calc, pUnit_conv),marker='.')\n", "lo_q, = axs1[1].plot(Pip_x_vec,Q_old,marker='.')\n", - "lo_pmin, = axs1[0].plot(Pip_x_vec,pipe.get_lowest_pressure_per_node(disp=True),c='red')\n", - "lo_pmax, = axs1[0].plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp=True),c='red')\n", + "lo_pmin, = axs1[0].plot(Pip_x_vec,pipe.get_lowest_pressure_per_node(disp_flag=True),c='red')\n", + "lo_pmax, = axs1[0].plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp_flag=True),c='red')\n", "lo_qmin, = axs1[1].plot(Pip_x_vec,pipe.get_lowest_flux_per_node(),c='red')\n", "lo_qmax, = axs1[1].plot(Pip_x_vec,pipe.get_highest_flux_per_node(),c='red')\n", "\n", @@ -191,7 +187,7 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -212,7 +208,7 @@ " # calculate the Con_T_ime evolution of the reservoir level within each pipeline timestep to avoid runaway numerical error\n", " for it_res in range(Res_nt):\n", " reservoir.timestep_reservoir_evolution() \n", - " level_vec[it_pipe] = reservoir.get_current_level() \n", + " level_vec[it_pipe] = reservoir.get_current_level() \n", " volume_vec[it_pipe] = reservoir.get_current_volume() \n", "\n", " # get the control variable\n", @@ -247,7 +243,6 @@ " v_old = pipe.get_current_velocity_distribution()\n", " Q_old = pipe.get_current_flux_distribution()\n", "\n", - "\n", " # plot some stuff\n", " # remove line-objects to autoscale axes (there is definetly a better way, but this works ¯\\_(ツ)_/¯ )\n", " lo_p.remove()\n", @@ -257,9 +252,9 @@ " lo_qmin.remove()\n", " lo_qmax.remove()\n", " # plot new pressure and velocity distribution in the pipeline\n", - " lo_p, = axs1[0].plot(Pip_x_vec,pipe.get_current_pressure_distribution(disp=True),marker='.',c='blue')\n", - " lo_pmin, = axs1[0].plot(Pip_x_vec,pipe.get_lowest_pressure_per_node(disp=True),c='red')\n", - " lo_pmax, = axs1[0].plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp=True),c='red')\n", + " lo_p, = axs1[0].plot(Pip_x_vec,pipe.get_current_pressure_distribution(disp_flag=True),marker='.',c='blue')\n", + " lo_pmin, = axs1[0].plot(Pip_x_vec,pipe.get_lowest_pressure_per_node(disp_flag=True),c='red')\n", + " lo_pmax, = axs1[0].plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp_flag=True),c='red')\n", " lo_q, = axs1[1].plot(Pip_x_vec,pipe.get_current_flux_distribution(),marker='.',c='blue')\n", " lo_qmin, = axs1[1].plot(Pip_x_vec,pipe.get_lowest_flux_per_node(),c='red')\n", " lo_qmax, = axs1[1].plot(Pip_x_vec,pipe.get_highest_flux_per_node(),c='red')\n", @@ -272,7 +267,7 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -316,8 +311,8 @@ "\n", "fig2,axs2 = plt.subplots(1,1)\n", "axs2.set_title('Min and Max Pressure')\n", - "axs2.plot(Pip_x_vec,pipe.get_lowest_pressure_per_node(disp=True),c='red')\n", - "axs2.plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp=True),c='red')\n", + "axs2.plot(Pip_x_vec,pipe.get_lowest_pressure_per_node(disp_flag=True),c='red')\n", + "axs2.plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp_flag=True),c='red')\n", "axs2.set_xlabel(r'$x$ [$\\mathrm{m}$]')\n", "axs2.set_ylabel(r'$p$ ['+pUnit_conv+']')\n", "\n", @@ -328,13 +323,6 @@ "axs2.set_xlabel(r'$x$ [$\\mathrm{m}$]')\n", "axs2.set_ylabel(r'$Q$ [$\\mathrm{m}^3/\\mathrm{s}$]')\n", "\n", - "# axs2[0,1].legend()\n", - "# axs2[1,0].legend()\n", - "# axs2[1,1].legend()\n", - "# # axs2[2,0].legend()\n", - "# # axs2[2,1].legend()\n", - "\n", - "\n", "fig2.tight_layout()\n", "plt.show()" ] diff --git a/Untertweng_mit_Pegelregler.ipynb b/Untertweng_mit_Pegelregler.ipynb index 93aa78a..64e9204 100644 --- a/Untertweng_mit_Pegelregler.ipynb +++ b/Untertweng_mit_Pegelregler.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 8, + "execution_count": 4, "metadata": {}, "outputs": [], "source": [ @@ -18,7 +18,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 5, "metadata": {}, "outputs": [], "source": [ @@ -30,12 +30,10 @@ "pUnit_calc = 'Pa' # [text] DO NOT CHANGE! for pressure conversion in print statements and plot labels \n", "pUnit_conv = 'mWS' # [text] for pressure conversion in print statements and plot labels\n", "\n", - "\n", " # for Turbine\n", - "Tur_Q_nenn = 0.85 # [m³/s] nominal flux of turbine \n", - "Tur_p_nenn = pressure_conversion(10.6,'bar',pUnit_calc) # [Pa] nominal pressure of turbine \n", - "Tur_closingTime = 90. # [s] closing time of turbine\n", - "\n", + "Tur_Q_nenn = 0.85 # [m³/s] nominal flux of turbine \n", + "Tur_p_nenn = pressure_conversion(10.6,'bar',pUnit_calc) # [Pa] nominal pressure of turbine \n", + "Tur_closingTime = 90. # [s] closing time of turbine\n", "\n", " # for PI controller\n", "Con_targetLevel = 8. # [m]\n", @@ -43,23 +41,21 @@ "Con_T_i = 1000. # [s] timespan in which a steady state error is corrected by the intergal term\n", "Con_deadbandRange = 0.05 # [m] Deadband range around targetLevel for which the controller does NOT intervene\n", "\n", - "\n", " # for pipeline\n", - "Pip_length = (535.+478.) # [m] length of pipeline\n", - "Pip_dia = 0.9 # [m] diameter of pipeline\n", - "Pip_area = Pip_dia**2/4*np.pi # [m²] crossectional area of pipeline\n", - "Pip_head = 105. # [m] hydraulic head of pipeline without reservoir\n", - "Pip_angle = np.arcsin(Pip_head/Pip_length) # [rad] elevation angle of pipeline \n", - "Pip_n_seg = 50 # [-] number of pipe segments in discretization\n", - "Pip_f_D = 0.014 # [-] Darcy friction factor\n", - "Pip_pw_vel = 500. # [m/s] propagation velocity of the pressure wave (pw) in the given pipeline\n", + "Pip_length = (535.+478.) # [m] length of pipeline\n", + "Pip_dia = 0.9 # [m] diameter of pipeline\n", + "Pip_area = Pip_dia**2/4*np.pi # [m²] crossectional area of pipeline\n", + "Pip_head = 105. # [m] hydraulic head of pipeline without reservoir\n", + "Pip_angle = np.arcsin(Pip_head/Pip_length) # [rad] elevation angle of pipeline \n", + "Pip_n_seg = 50 # [-] number of pipe segments in discretization\n", + "Pip_f_D = 0.014 # [-] Darcy friction factor\n", + "Pip_pw_vel = 500. # [m/s] propagation velocity of the pressure wave (pw) in the given pipeline\n", " # derivatives of the pipeline constants\n", - "Pip_dx = Pip_length/Pip_n_seg # [m] length of each pipe segment\n", - "Pip_dt = Pip_dx/Pip_pw_vel # [s] timestep according to method of characteristics\n", - "Pip_nn = Pip_n_seg+1 # [1] number of nodes\n", - "Pip_x_vec = np.arange(0,Pip_nn,1)*Pip_dx # [m] vector holding the distance of each node from the upstream reservoir along the pipeline\n", - "Pip_h_vec = np.arange(0,Pip_nn,1)*Pip_head/Pip_n_seg # [m] vector holding the vertival distance of each node from the upstream reservoir\n", - "\n", + "Pip_dx = Pip_length/Pip_n_seg # [m] length of each pipe segment\n", + "Pip_dt = Pip_dx/Pip_pw_vel # [s] timestep according to method of characteristics\n", + "Pip_nn = Pip_n_seg+1 # [1] number of nodes\n", + "Pip_x_vec = np.arange(0,Pip_nn,1)*Pip_dx # [m] vector holding the distance of each node from the upstream reservoir along the pipeline\n", + "Pip_h_vec = np.arange(0,Pip_nn,1)*Pip_head/Pip_n_seg # [m] vector holding the vertival distance of each node from the upstream reservoir\n", "\n", " # for reservoir\n", "Res_area_base = 74. # [m²] total base are of the cuboid reservoir \n", @@ -71,28 +67,28 @@ "Res_dt = Pip_dt/Res_nt # [s] harmonised timestep of reservoir time evolution\n", "\n", " # for general simulation\n", - "flux_init = Tur_Q_nenn/1.1 # [m³/s] initial flux through whole system for steady state initialization \n", - "level_init = Con_targetLevel # [m] initial water level in upstream reservoir for steady state initialization\n", - "simTime_target = 600. # [s] target for total simulation time (will vary slightly to fit with Pip_dt)\n", - "nt = int(simTime_target//Pip_dt) # [1] Number of timesteps of the whole system\n", - "t_vec = np.arange(0,nt+1,1)*Pip_dt # [s] time vector. At each step of t_vec the system parameters are stored\n" + "flux_init = Tur_Q_nenn/1.1 # [m³/s] initial flux through whole system for steady state initialization \n", + "level_init = Con_targetLevel # [m] initial water level in upstream reservoir for steady state initialization\n", + "simTime_target = 100. # [s] target for total simulation time (will vary slightly to fit with Pip_dt)\n", + "nt = int(simTime_target//Pip_dt) # [1] Number of timesteps of the whole system\n", + "t_vec = np.arange(0,nt+1,1)*Pip_dt # [s] time vector. At each step of t_vec the system parameters are stored\n" ] }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# create objects\n", "\n", "# Upstream reservoir\n", - "reservoir = Ausgleichsbecken_class(Res_area_base,Res_area_out,Res_dt,Res_level_crit_lo,Res_level_crit_hi,rho)\n", + "reservoir = Ausgleichsbecken_class(Res_area_base,Res_area_out,Res_dt,pUnit_conv,Res_level_crit_lo,Res_level_crit_hi,rho)\n", "reservoir.set_steady_state(flux_init,level_init)\n", "\n", "# pipeline\n", - "pipe = Druckrohrleitung_class(Pip_length,Pip_dia,Pip_n_seg,Pip_angle,Pip_f_D,Pip_pw_vel,Pip_dt,pUnit_conv,rho)\n", - "pipe.set_steady_state(flux_init,level_init,Res_area_base,Pip_x_vec,Pip_h_vec)\n", + "pipe = Druckrohrleitung_class(Pip_length,Pip_dia,Pip_head,Pip_n_seg,Pip_f_D,Pip_pw_vel,Pip_dt,pUnit_conv,rho)\n", + "pipe.set_steady_state(flux_init,reservoir.get_current_pressure())\n", "\n", "# downstream turbine\n", "turbine = Francis_Turbine(Tur_Q_nenn,Tur_p_nenn,Tur_closingTime,Pip_dt,pUnit_conv)\n", @@ -109,7 +105,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 7, "metadata": {}, "outputs": [], "source": [ @@ -158,7 +154,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 9, "metadata": {}, "outputs": [], "source": [ @@ -176,8 +172,8 @@ "axs1[1].set_ylabel(r'$Q$ [$\\mathrm{m}^3 / \\mathrm{s}$]')\n", "lo_p, = axs1[0].plot(Pip_x_vec,pressure_conversion(p_old,pUnit_calc, pUnit_conv),marker='.')\n", "lo_q, = axs1[1].plot(Pip_x_vec,Q_old,marker='.')\n", - "lo_pmin, = axs1[0].plot(Pip_x_vec,pipe.get_lowest_pressure_per_node(disp=True),c='red')\n", - "lo_pmax, = axs1[0].plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp=True),c='red')\n", + "lo_pmin, = axs1[0].plot(Pip_x_vec,pipe.get_lowest_pressure_per_node(disp_flag=True),c='red')\n", + "lo_pmax, = axs1[0].plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp_flag=True),c='red')\n", "lo_qmin, = axs1[1].plot(Pip_x_vec,pipe.get_lowest_flux_per_node(),c='red')\n", "lo_qmax, = axs1[1].plot(Pip_x_vec,pipe.get_highest_flux_per_node(),c='red')\n", "\n", @@ -191,7 +187,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 10, "metadata": {}, "outputs": [], "source": [ @@ -247,7 +243,6 @@ " v_old = pipe.get_current_velocity_distribution()\n", " Q_old = pipe.get_current_flux_distribution()\n", "\n", - "\n", " # plot some stuff\n", " # remove line-objects to autoscale axes (there is definetly a better way, but this works ¯\\_(ツ)_/¯ )\n", " lo_p.remove()\n", @@ -257,9 +252,9 @@ " lo_qmin.remove()\n", " lo_qmax.remove()\n", " # plot new pressure and velocity distribution in the pipeline\n", - " lo_p, = axs1[0].plot(Pip_x_vec,pipe.get_current_pressure_distribution(disp=True),marker='.',c='blue')\n", - " lo_pmin, = axs1[0].plot(Pip_x_vec,pipe.get_lowest_pressure_per_node(disp=True),c='red')\n", - " lo_pmax, = axs1[0].plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp=True),c='red')\n", + " lo_p, = axs1[0].plot(Pip_x_vec,pipe.get_current_pressure_distribution(disp_flag=True),marker='.',c='blue')\n", + " lo_pmin, = axs1[0].plot(Pip_x_vec,pipe.get_lowest_pressure_per_node(disp_flag=True),c='red')\n", + " lo_pmax, = axs1[0].plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp_flag=True),c='red')\n", " lo_q, = axs1[1].plot(Pip_x_vec,pipe.get_current_flux_distribution(),marker='.',c='blue')\n", " lo_qmin, = axs1[1].plot(Pip_x_vec,pipe.get_lowest_flux_per_node(),c='red')\n", " lo_qmax, = axs1[1].plot(Pip_x_vec,pipe.get_highest_flux_per_node(),c='red')\n", @@ -272,7 +267,7 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 11, "metadata": {}, "outputs": [], "source": [ @@ -316,8 +311,8 @@ "\n", "fig2,axs2 = plt.subplots(1,1)\n", "axs2.set_title('Min and Max Pressure')\n", - "axs2.plot(Pip_x_vec,pipe.get_lowest_pressure_per_node(disp=True),c='red')\n", - "axs2.plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp=True),c='red')\n", + "axs2.plot(Pip_x_vec,pipe.get_lowest_pressure_per_node(disp_flag=True),c='red')\n", + "axs2.plot(Pip_x_vec,pipe.get_highest_pressure_per_node(disp_flag=True),c='red')\n", "axs2.set_xlabel(r'$x$ [$\\mathrm{m}$]')\n", "axs2.set_ylabel(r'$p$ ['+pUnit_conv+']')\n", "\n", diff --git a/functions/pressure_conversion.py b/functions/pressure_conversion.py index d4a91f9..12c9048 100644 --- a/functions/pressure_conversion.py +++ b/functions/pressure_conversion.py @@ -27,8 +27,6 @@ def pa_to_torr(p): def pa_to_atm(p): return p*1/(101.325*1e3) - # converstion function - def pa_to_psi(p): return p/6894.8