{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "51266afe",
   "metadata": {},
   "source": [
    "# <center><b> Master thesis : Analysis of the saturated absorption of the HD molecule in an optical cavity </b></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20aab70c",
   "metadata": {},
   "source": [
    "# <center><b> Lamb Dip </b></center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "314ee040",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import integrate\n",
    "from qutip import Qobj, Bloch, about, basis, mesolve, fock, qutip, Options\n",
    "\n",
    "class Atom(object):\n",
    "    '''\n",
    "    This class creates an atom object containing all the informations needed to compute the spectrum.\n",
    "    The atom must be instanciated with properties expressed in atomic units.\n",
    "\n",
    "    Atomic units:\n",
    "\n",
    "    hbar = 1\n",
    "    m_e = 1\n",
    "    e = 1\n",
    "    4 pi epsilon_0 = 1\n",
    "\n",
    "    More precisely:\n",
    "    - Transition frequency f_0 in E_H/hbar (1 E_H = 1 Hartree, hbar = constante de Planck réduite)\n",
    "    - Dipole constant D12 in e*a_0 (e = elementary electric charge, a_0 = Bohr)\n",
    "    '''\n",
    "\n",
    "    def __init__(self, default_mass = 1, default_transition_frequency = 5, default_dipole_moment = 1,\n",
    "                 default_quantum_state = [[1, 0], [0, 0]]):\n",
    "\n",
    "        #Universal constants in atomic units\n",
    "        self.h = 2*np.pi\n",
    "        self.c = 137\n",
    "        self.epsilon_0 = 1 / (4*np.pi)\n",
    "\n",
    "        #Constants specific to the atom\n",
    "        self.mass = default_mass\n",
    "        self.f_0 = default_transition_frequency\n",
    "        self.D12 = default_dipole_moment\n",
    "\n",
    "        # Basis\n",
    "        self.psi_1 = fock(2, 0)\n",
    "        self.psi_2 = fock(2, 1)\n",
    "\n",
    "        #Quantum state\n",
    "        self.Qstate = Qobj(default_quantum_state)\n",
    "\n",
    "    def get_h(self):\n",
    "        return self.h\n",
    "\n",
    "    def get_c(self):\n",
    "        return self.c\n",
    "\n",
    "    def get_epsilon_0(self):\n",
    "        return self.epsilon_0\n",
    "\n",
    "    def getTransitionFrequency(self):\n",
    "        return self.f_0\n",
    "\n",
    "    def getDipoleMoment(self):\n",
    "        return self.D12\n",
    "\n",
    "    def getDipoleMoment_X12(self):\n",
    "        X12 = np.sqrt(self.D12**2 /3)\n",
    "        return X12\n",
    "\n",
    "    def getAngularFrequency(self):\n",
    "        omega_0 = 2*np.pi*self.f_0\n",
    "        return omega_0\n",
    "\n",
    "    def getWaveLength(self):\n",
    "        lambda_0 = self.c/self.f_0\n",
    "        return lambda_0\n",
    "\n",
    "    def getEinsteinCoefficient_A21(self):\n",
    "        omega_0 = 2*np.pi*self.f_0\n",
    "        A21 = (4/3) * (1/self.c**3) * omega_0**3 * self.D12**2\n",
    "        return A21\n",
    "\n",
    "    def getSaturationIntensity(self):\n",
    "        omega_0 = 2*np.pi*self.f_0\n",
    "        A21 = (4/3) * (1/self.c**3) * omega_0**3 * self.D12**2\n",
    "        lambda_0 = self.c/self.f_0\n",
    "        I_sat = (np.pi*self.h*self.c*A21)/(3*(lambda_0**3))\n",
    "        return I_sat\n",
    "\n",
    "    def getQuantumState(self):\n",
    "        return self.Qstate\n",
    "\n",
    "    def updateQuantumState(self, new_Qstate):\n",
    "        self.Qstate = new_Qstate\n",
    "\n",
    "    def getBasis(self):\n",
    "        return self.psi_1, self.psi_2\n",
    "\n",
    "    def getMass(self):\n",
    "        return self.mass\n",
    "\n",
    "class Laser(object):\n",
    "    '''\n",
    "    This class creates a laser object containing all the informations needed to compute the spectrum.\n",
    "    The laser must be instanciated with properties expressed in atomic units units.\n",
    "\n",
    "    Atomic units:\n",
    "\n",
    "    hbar = 1\n",
    "    m_e = 1\n",
    "    e = 1\n",
    "    4 pi epsilon_0 = 1\n",
    "\n",
    "    '''\n",
    "\n",
    "    def __init__(self, default_frequency = 5.01, default_intensity = 1, default_phi_time = 0.0, default_phi_space = 0.0):\n",
    "\n",
    "        #Universal constants in atomic units\n",
    "        self.h = 2*np.pi\n",
    "        self.c = 137\n",
    "        self.epsilon_0 = 4*np.pi\n",
    "\n",
    "        #Constants specific to the laser\n",
    "        self.f = default_frequency\n",
    "        self.I = default_intensity\n",
    "        self.phi_time = default_phi_time\n",
    "        self.phi_space = default_phi_space\n",
    "\n",
    "    def getFrequency(self):\n",
    "        return self.f\n",
    "\n",
    "    def getIntensity(self):\n",
    "        return self.I\n",
    "\n",
    "    def getPhaseTime(self):\n",
    "        return self.phi_time\n",
    "\n",
    "    def getPhaseSpace(self):\n",
    "        return self.phi_space\n",
    "\n",
    "    def getAngularFrequency(self):\n",
    "        omega = 2*np.pi*self.f\n",
    "        return omega\n",
    "\n",
    "    def getElectricFieldAmplitude(self):\n",
    "        E = np.sqrt(2*self.I/(self.c*self.epsilon_0))\n",
    "        return E\n",
    "\n",
    "class GazeousMedium(object):\n",
    "    '''\n",
    "    This class creates a gazeous medium. This medium takes the properties of an Atom() or Molecule() object\n",
    "    and compute the velocity classes according to the Maxwell-boltzmann distribution.\n",
    "    '''\n",
    "    def __init__(self, default_quantum_system = Atom(), default_temperature = 300,\n",
    "                 default_pressure = 3.44e-9, default_density = 1):\n",
    "        self.Qsyst = default_quantum_system\n",
    "\n",
    "        self.T = default_temperature\n",
    "        self.P = default_pressure\n",
    "        self.D = default_density\n",
    "\n",
    "        #Boltzmann constant\n",
    "        self.k_b = 3.167*1e-6\n",
    "\n",
    "    def getQuantumSystem(self):\n",
    "        return self.Qsyst\n",
    "\n",
    "    def distMaxwellBoltzmann(self):\n",
    "        'Compute velocity classes'\n",
    "        M = self.Qsyst.getMass()\n",
    "        v_p = np.sqrt(2*self.k_b*self.T/M)\n",
    "\n",
    "        n = 1/2\n",
    "        scale = 3\n",
    "        subdivision = 40\n",
    "        v = np.zeros(2*scale*subdivision)\n",
    "        for index in range(scale):\n",
    "            start = subdivision*index\n",
    "            v[start:start + subdivision] = np.linspace(-n*v_p/(10**(index)),\n",
    "                                                       -n*v_p/(10**(index+1)),\n",
    "                                                       subdivision)\n",
    "            v[len(v) - subdivision - start:len(v)-start] = np.linspace(n*v_p/(10**(index+1)),\n",
    "                                                                       n*v_p/(10**(index)),\n",
    "                                                                       subdivision)\n",
    "\n",
    "        D_v = np.zeros(len(v))\n",
    "        for index in range(len(v)):\n",
    "            D_v[index] = self.D*(M/(2*np.pi*self.k_b*self.T))**(1/2)*np.exp(-(M*v[index]**2)/(2*self.k_b*self.T))\n",
    "\n",
    "        return D_v, v\n",
    "\n",
    "    def distMB_AdaptativeSampling(self, v_shift):\n",
    "        'Compute velocity classes'\n",
    "        M = self.Qsyst.getMass()\n",
    "        v_p = np.sqrt(2*self.k_b*self.T/M)\n",
    "\n",
    "        n = 2\n",
    "        scale = 3\n",
    "        subdivision = 20\n",
    "        v_plus = np.zeros(2*(scale*subdivision-(scale-1)))\n",
    "        v_minus = np.zeros(2*(scale*subdivision-(scale-1)))\n",
    "        start = 0\n",
    "        for index in range(scale):\n",
    "            v_plus[start:start + subdivision] = np.linspace(-n*v_p/(10**(index)) + v_shift,\n",
    "                                                            -n*v_p/(10**(index+1)) + v_shift,\n",
    "                                                            subdivision)\n",
    "            v_plus[len(v_plus)-subdivision-start:len(v_plus)-start] = np.linspace(n*v_p/(10**(index+1)) + v_shift,\n",
    "                                                                                  n*v_p/(10**(index)) + v_shift,\n",
    "                                                                                  subdivision)\n",
    "            \n",
    "            v_minus[start:start+subdivision] = np.linspace(-n*v_p/(10**(index)) - v_shift,\n",
    "                                                           -n*v_p/(10**(index+1)) - v_shift,\n",
    "                                                           subdivision)\n",
    "            v_minus[len(v_minus)-subdivision-start:len(v_minus)-start] = np.linspace(n*v_p/(10**(index+1)) - v_shift,\n",
    "                                                                                     n*v_p/(10**(index)) - v_shift,\n",
    "                                                                                     subdivision)\n",
    "            start = start + subdivision - 1\n",
    "        \n",
    "        if v_minus[-1] < 0.0:\n",
    "            length_minus = len(v_minus)\n",
    "        else:\n",
    "            length_minus = np.argwhere(v_minus >= 0.0)[0][0]\n",
    "\n",
    "        if v_plus[0] > 0.0:\n",
    "            start_v_plus = 0\n",
    "            length_plus = len(v_plus)\n",
    "        else:\n",
    "            start_v_plus = np.argwhere(v_plus > 0.0)[0][0]\n",
    "            length_plus = len(v_plus) - start_v_plus\n",
    "\n",
    "        v = np.zeros(length_plus+length_minus)\n",
    "        for index_1 in range(length_minus):\n",
    "            v[index_1] = v_minus[index_1]\n",
    "        for index_2 in range(length_plus):\n",
    "            v[index_1+1+index_2] = v_plus[start_v_plus + index_2]\n",
    "\n",
    "        D_v = np.zeros(len(v))\n",
    "        for index in range(len(v)):\n",
    "            D_v[index] = self.D*(M/(2*np.pi*self.k_b*self.T))**(1/2)*np.exp(-(M*v[index]**2)/(2*self.k_b*self.T))\n",
    "\n",
    "        ### Debugging ###\n",
    "        #D_v_plus = np.zeros(len(v_plus))\n",
    "        #for index in range(len(v_plus)):\n",
    "        #    D_v_plus[index] = self.D*(M/(2*np.pi*self.k_b*self.T))**(1/2)*np.exp(-(M*v_plus[index]**2)/(2*self.k_b*self.T))\n",
    "        #\n",
    "        #D_v_minus = np.zeros(len(v_minus))\n",
    "        #for index in range(len(v_minus)):\n",
    "        #    D_v_minus[index] = self.D*(M/(2*np.pi*self.k_b*self.T))**(1/2)*np.exp(-(M*v_minus[index]**2)/(2*self.k_b*self.T))\n",
    "        #return D_v, v, D_v_plus, v_plus, D_v_minus, v_minus\n",
    "        ### Debugging ###\n",
    "\n",
    "        return D_v, v\n",
    "\n",
    "class AbsorptionSpectrum(object):\n",
    "    '''\n",
    "    This class creates a python object that can compute the absorption spectrum\n",
    "    of a given Atom() or Molecule() and a given Laser().\n",
    "    '''\n",
    "\n",
    "    def __init__(self, default_medium = GazeousMedium(), default_laser = Laser()):\n",
    "\n",
    "        #Medium\n",
    "        self.medium = default_medium\n",
    "\n",
    "        #Laser\n",
    "        self.laser = default_laser\n",
    "\n",
    "    def spectrum_v(self, spectral_window, velocity):\n",
    "\n",
    "        #Quantum system\n",
    "        Qsyst = self.medium.getQuantumSystem()\n",
    "        psi_1, psi_2 = Qsyst.getBasis()\n",
    "        omega_0 = Qsyst.getAngularFrequency()\n",
    "        X12 = Qsyst.getDipoleMoment_X12()\n",
    "        A21 = Qsyst.getEinsteinCoefficient_A21()\n",
    "        v = velocity\n",
    "\n",
    "        #Universal constants\n",
    "        h = Qsyst.get_h()\n",
    "        c = Qsyst.get_c()\n",
    "        epsilon_0 = Qsyst.get_epsilon_0()\n",
    "\n",
    "        #Laser\n",
    "        E = self.laser.getElectricFieldAmplitude()\n",
    "        I = self.laser.getIntensity()\n",
    "        phi_time = self.laser.getPhaseTime()\n",
    "        phi_space = self.laser.getPhaseSpace()\n",
    "\n",
    "        phi_1 = phi_time + phi_space\n",
    "        phi_2 = phi_time - phi_space\n",
    "\n",
    "        #Collapse operators\n",
    "        c_ops = []\n",
    "        c_ops.append(np.sqrt(A21) * psi_1 * psi_2.dag())\n",
    "\n",
    "        length_sw = len(spectral_window)\n",
    "        alpha = np.zeros(length_sw)\n",
    "\n",
    "        if v!=0.0:\n",
    "\n",
    "            #Rabi frequency\n",
    "            Omega_R = -X12*(E/2)\n",
    "\n",
    "            #Hamiltonian (without time dependance)\n",
    "            H1 = Qobj(np.array([[0, Omega_R/2], [0, 0]]))\n",
    "            H2 = Qobj(np.array([[0, 0], [np.conjugate(Omega_R/2), 0]]))\n",
    "\n",
    "            for i in range(length_sw):\n",
    "\n",
    "                # Initial density matrix\n",
    "                Qstate_initial = Qobj([[1, 0], [0, 0]])\n",
    "\n",
    "                omega = spectral_window[i]\n",
    "\n",
    "                # Doppler\n",
    "                k = omega / c\n",
    "                delta_blue = omega + k*v - omega_0\n",
    "                delta_red = omega - k*v - omega_0\n",
    "\n",
    "                #Time\n",
    "                T_kv = abs(2*np.pi/(k*v))\n",
    "                n_time_steady_state = 10\n",
    "                n_time_oscillations = 4\n",
    "                time_end = n_time_steady_state*(1/A21) + n_time_oscillations*(T_kv/2)\n",
    "                T_Omega_R_demi = abs(2*(2*np.pi)/Omega_R)\n",
    "                time_subdivision = max(int((time_end/T_Omega_R_demi)*20000), int((time_end/T_kv)*2000))\n",
    "                tlist = np.linspace(0, time_end, time_subdivision)\n",
    "\n",
    "                #Solution of the Lindblad equation\n",
    "                H_time = [[H1, np.exp(1j*(delta_blue*tlist+phi_1)) + np.exp(1j*(delta_red*tlist+phi_2))],\n",
    "                          [H2, np.exp(-1j*(delta_blue*tlist+phi_1)) + np.exp(-1j*(delta_red*tlist+phi_2))]]\n",
    "                rho = mesolve(H_time, Qstate_initial, tlist, c_ops, [])\n",
    "\n",
    "                rho22_LB = np.zeros(len(rho.states))\n",
    "                for index in range(0,len(rho.states)):\n",
    "                    rho22_LB[index] = rho.states[index].full()[1][1].real\n",
    "\n",
    "                n = 4\n",
    "                t_steady_state = np.abs(tlist - n*(1/A21)).argmin()\n",
    "\n",
    "                # number density\n",
    "                N = 1\n",
    "\n",
    "                # Loss\n",
    "                I_loss = omega*N*A21/(tlist[-1] - tlist[t_steady_state]) * \\\n",
    "                        integrate.simpson(rho22_LB[t_steady_state:len(tlist)], tlist[t_steady_state:len(tlist)])\n",
    "\n",
    "                alpha[i] = I_loss/I\n",
    "        else:\n",
    "\n",
    "            for i in range(length_sw):\n",
    "\n",
    "                omega = spectral_window[i]\n",
    "                delta = omega - omega_0\n",
    "                k = omega / c\n",
    "\n",
    "                lambda_demi_SW = np.pi/k\n",
    "                subdivision_x = 11\n",
    "                positions = np.linspace(0,lambda_demi_SW, subdivision_x)\n",
    "                alpha_0 = np.zeros(len(positions))\n",
    "\n",
    "                for index_loop in range(len(positions)):\n",
    "\n",
    "                    E_x = E*np.sin(k*positions[index_loop])\n",
    "\n",
    "                    #Rabi frequency\n",
    "                    Omega_R = -X12*E_x\n",
    "                    Omega_R_max = -X12*E\n",
    "\n",
    "                    #Hamiltonian (without time dependance)\n",
    "                    H1 = Qobj(np.array([[0, Omega_R/2], [0, 0]]))\n",
    "                    H2 = Qobj(np.array([[0, 0], [np.conjugate(Omega_R/2), 0]]))\n",
    "\n",
    "                    #Time\n",
    "                    n_time = 20\n",
    "                    time_end = n_time*(1/A21)\n",
    "                    T_Omega_R_max_demi = abs(2*(2*np.pi)/Omega_R_max)\n",
    "                    time_subdivision = int((time_end/T_Omega_R_max_demi)*20000)\n",
    "                    tlist = np.linspace(0, time_end, time_subdivision)\n",
    "\n",
    "                    # Initial density matrix\n",
    "                    Qstate_initial = Qobj([[1, 0], [0, 0]])\n",
    "\n",
    "                    #Solution of the Lindblad equation\n",
    "                    H_time = [[H1, np.exp(1j*(delta*tlist+phi_time))], [H2, np.exp(-1j*(delta*tlist+phi_time))]]\n",
    "                    rho = mesolve(H_time, Qstate_initial, tlist, c_ops, [])\n",
    "\n",
    "                    rho22_LB = np.zeros(len(rho.states))\n",
    "                    for index in range(0,len(rho.states)):\n",
    "                        rho22_LB[index] = rho.states[index].full()[1][1].real\n",
    "\n",
    "                    n = 10\n",
    "                    t_steady_state = np.abs(tlist - n*(1/A21)).argmin()\n",
    "\n",
    "                    # number density\n",
    "                    N = 1\n",
    "\n",
    "                    # Loss\n",
    "                    I_loss = omega*N*A21/(tlist[-1] - tlist[t_steady_state]) * \\\n",
    "                            integrate.simpson(rho22_LB[t_steady_state:len(tlist)], tlist[t_steady_state:len(tlist)])\n",
    "\n",
    "                    alpha_0[index_loop] = I_loss/I\n",
    "\n",
    "                alpha[i] = integrate.simpson(alpha_0/lambda_demi_SW, positions)\n",
    "\n",
    "        return alpha\n",
    "\n",
    "    def spectrum(self, omega):\n",
    "\n",
    "        #Quantum system\n",
    "        Qsyst = self.medium.getQuantumSystem()\n",
    "        psi_1, psi_2 = Qsyst.getBasis()\n",
    "        omega_0 = Qsyst.getAngularFrequency()\n",
    "        X12 = Qsyst.getDipoleMoment_X12()\n",
    "        A21 = Qsyst.getEinsteinCoefficient_A21()\n",
    "\n",
    "        #Universal constants\n",
    "        h = Qsyst.get_h()\n",
    "        c = Qsyst.get_c()\n",
    "        epsilon_0 = Qsyst.get_epsilon_0()\n",
    "\n",
    "        #Laser\n",
    "        E = self.laser.getElectricFieldAmplitude()\n",
    "        I = self.laser.getIntensity()\n",
    "        phi_time = self.laser.getPhaseTime()\n",
    "        phi_space = self.laser.getPhaseSpace()\n",
    "\n",
    "        phi_1 = phi_time + phi_space\n",
    "        phi_2 = phi_time - phi_space\n",
    "\n",
    "        #Collapse operators\n",
    "        c_ops = []\n",
    "        c_ops.append(np.sqrt(A21) * psi_1 * psi_2.dag())\n",
    "\n",
    "        #Velocityn sampling\n",
    "        k = omega / c\n",
    "        v_shift = abs(omega_0 - omega)/k\n",
    "        D_v, v = self.medium.distMB_AdaptativeSampling(v_shift)\n",
    "        alpha_v = np.zeros(len(v))\n",
    "\n",
    "        #Rabi frequency\n",
    "        Omega_R = -X12*(E/2)\n",
    "\n",
    "        #Hamiltonian (without time dependance)\n",
    "        H1 = Qobj(np.array([[0, Omega_R/2], [0, 0]]))\n",
    "        H2 = Qobj(np.array([[0, 0], [np.conjugate(Omega_R/2), 0]]))\n",
    "\n",
    "        for i in range(len(v)):\n",
    "\n",
    "            # Initial density matrix\n",
    "            Qstate_initial = Qobj([[1, 0], [0, 0]])\n",
    "\n",
    "            # Doppler\n",
    "            delta_blue = omega + k*v[i] - omega_0\n",
    "            delta_red = omega - k*v[i] - omega_0\n",
    "\n",
    "            #Time\n",
    "            T_kv = abs(2*np.pi/(k*v[i]))\n",
    "            n_time_steady_state = 10\n",
    "            n_time_oscillations = 4\n",
    "            time_end = n_time_steady_state*(1/A21) + n_time_oscillations*(T_kv/2)\n",
    "            T_Omega_R_demi = abs(2*(2*np.pi)/Omega_R)\n",
    "            time_subdivision = max(int((time_end/T_Omega_R_demi)*20000), int((time_end/T_kv)*2000))\n",
    "            tlist = np.linspace(0, time_end, time_subdivision)\n",
    "\n",
    "            #Solution of the Lindblad equation\n",
    "            H_time = [[H1, np.exp(1j*(delta_blue*tlist+phi_1)) + np.exp(1j*(delta_red*tlist+phi_2))],\n",
    "                      [H2, np.exp(-1j*(delta_blue*tlist+phi_1)) + np.exp(-1j*(delta_red*tlist+phi_2))]]\n",
    "            rho = mesolve(H_time, Qstate_initial, tlist, c_ops, [])\n",
    "\n",
    "            rho22_LB = np.zeros(len(rho.states))\n",
    "            for index in range(0,len(rho.states)):\n",
    "                rho22_LB[index] = rho.states[index].full()[1][1].real\n",
    "\n",
    "            n = 4\n",
    "            t_steady_state = np.abs(tlist - n*(1/A21)).argmin()\n",
    "\n",
    "            # number density\n",
    "            N = 1\n",
    "\n",
    "            # Loss\n",
    "            I_loss = omega*N*A21/(tlist[-1] - tlist[t_steady_state]) * \\\n",
    "                    integrate.simpson(rho22_LB[t_steady_state:len(tlist)], tlist[t_steady_state:len(tlist)])\n",
    "\n",
    "            alpha_v[i] = I_loss/I\n",
    "\n",
    "        alpha = integrate.simpson(D_v*alpha_v, v)\n",
    "\n",
    "        return alpha\n",
    "        ### Debugging ###\n",
    "        #return alpha_v\n",
    "\n",
    "    def spectrum_broadened_MB(self, spectral_window):\n",
    "\n",
    "        length_sw = len(spectral_window)\n",
    "\n",
    "        palier = 0.0\n",
    "        step = 5\n",
    "        alpha = np.zeros(length_sw)\n",
    "        for i in range(length_sw):\n",
    "            alpha[i] = self.spectrum(spectral_window[i])\n",
    "\n",
    "            if (i/length_sw*100) >= palier:\n",
    "                print(i/length_sw*100, \"%\")\n",
    "                palier += step\n",
    "\n",
    "        return alpha\n",
    "\n",
    "#######################\n",
    "## Instanciations ##\n",
    "\n",
    "# Gaseous medium\n",
    "dummy_atom = Atom(default_mass=0.001)\n",
    "I_sat = dummy_atom.getSaturationIntensity()\n",
    "dummy_medium = GazeousMedium(dummy_atom)\n",
    "\n",
    "# Spectral window\n",
    "omega_0 = dummy_atom.getAngularFrequency()\n",
    "width = 0.05\n",
    "frequency_window = np.linspace(omega_0-width, omega_0+width, 21)\n",
    "np.savetxt('frequency_window_21pt.txt', frequency_window)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c797d407",
   "metadata": {},
   "outputs": [],
   "source": [
    "#######################\n",
    "## RUNS ##\n",
    "# Varying LASER intensity\n",
    "dummy_laser_0p5 = Laser(default_intensity= 0.5*I_sat)\n",
    "dummy_spectrum_0p5 = AbsorptionSpectrum(dummy_medium, dummy_laser_0p5)\n",
    "alpha_0p5 = dummy_spectrum_0p5.spectrum_broadened_MB(frequency_window)\n",
    "np.savetxt('alpha_gas_0p5Isat.txt', alpha_0p5)\n",
    "\n",
    "dummy_laser_1 = Laser(default_intensity= 1*I_sat)\n",
    "dummy_spectrum_1 = AbsorptionSpectrum(dummy_medium, dummy_laser_1)\n",
    "alpha_1 = dummy_spectrum_1.spectrum_broadened_MB(frequency_window)\n",
    "np.savetxt('alpha_gas_1Isat.txt', alpha_1)\n",
    "\n",
    "dummy_laser_50 = Laser(default_intensity= 50*I_sat)\n",
    "dummy_spectrum_50 = AbsorptionSpectrum(dummy_medium, dummy_laser_50)\n",
    "alpha_50 = dummy_spectrum_50.spectrum_broadened_MB(frequency_window)\n",
    "np.savetxt('alpha_gas_50Isat.txt', alpha_50)\n",
    "\n",
    "dummy_laser_100 = Laser(default_intensity= 100*I_sat)\n",
    "dummy_spectrum_100 = AbsorptionSpectrum(dummy_medium, dummy_laser_100)\n",
    "alpha_100 = dummy_spectrum_100.spectrum_broadened_MB(frequency_window)\n",
    "np.savetxt('alpha_gas_100Isat.txt', alpha_100)\n",
    "\n",
    "dummy_laser_500 = Laser(default_intensity= 500*I_sat)\n",
    "dummy_spectrum_500 = AbsorptionSpectrum(dummy_medium, dummy_laser_500)\n",
    "alpha_500 = dummy_spectrum_500.spectrum_broadened_MB(frequency_window)\n",
    "np.savetxt('alpha_gas_500Isat.txt', alpha_500)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7b8fe6c",
   "metadata": {},
   "source": [
    "## Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "cff0210e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.ticker import MaxNLocator\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "f03638b0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmIAAAG9CAYAAAC21hqBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAADaBElEQVR4nOzdd1hTZxsH4F/C3kMQRAEFB+4B4kTEPaq1bquto666tWqrVgH3aN174qyj1m2tE3ECggMUFRd7yt4Zz/fHkSAfKCMJSeC9ryuXcnJy3udkPuedPCIiMAzDMAzDMBWOr+gAGIZhGIZhqiqWiDEMwzAMwygIS8QYhmEYhmEUhCViDMMwDMMwCsISMYZhGIZhGAVhiRjDMAzDMIyCsESMYRiGYRhGQVgixjAMwzAMoyDqig6gKhGLxYiOjoaBgQF4PJ6iw2EYhmEYphSICOnp6bCysgKfL9s6LJaIVaDo6GhYW1srOgyGYRiGYcohIiICtWrVkukxWSJWgQwMDABwL6ShoaGCo2EYhmEYpjTS0tJgbW0t+R2XJZaIVaD85khDQ0OWiDEMwzCMipFHtyLWWZ9hGIZhGEZBWCLGMAzDMAyjICwRYxiGYRiGURCWiDEMwzAMwygIS8QYhmEYhmEUhCViDMMwDMMwCsISMYZhGIZhGAVhiRjDMAzDMIyCsESMYRiGYRhGQVgixjAMwzAMoyBKn4hFRUVh48aN6NGjB2xsbKCpqQlLS0sMGjQIvr6+ZT5eeno63N3d0aRJE+jq6sLY2BitWrWCp6dnsfv7+/ujT58+MDExgZ6eHpydnXHs2DFpT4thGIZhGAY8IiJFB/E1v/32G9asWQN7e3u4urqievXqCA0NxdmzZ0FE+OuvvzB06NBSHSs8PBxdunTBu3fv0K1bN7Rs2RK5ubl48+YNwsPD8ezZs0L7e3t7o2fPntDU1MTw4cNhZGSEf/75B+/fv8eKFSuwcOHCMp1LWloajIyMkJqaytaaZBiGYRgVIc/fb6VPxP755x+Ym5vDxcWl0PY7d+6ga9euMDAwQHR0NLS0tL56HJFIhHbt2iE4OBiXLl2Cm5tbofuFQiHU1dUL/e3g4IDIyEg8ePAALVu2BMDVqLVr1w6vXr3CixcvUK9evVKfC0vEGKYSEgqBoCDg6VPAxgZo3RowMFB0VAzDyJA8f7/VS95FsQYOHFjsdhcXF7i5ueHq1asICgqCk5PTV4/z999/w9/fH4sXLy6ShAEolIQBwM2bN/H27VuMHTtWkoQBgIGBARYvXozhw4fjwIEDWLlyZTnOimEqQE4OcP06cm7fhqaaGvhaWoCmpmxvWloAX+l7OMhWfDzw8CHw4AF38/cHsrIkdxMAXtOmQNu2BTcHh8r5PAmFQFISkJAAJCYWulFyMkQCAdRFIm4/gQAfMzKQmZsLUz4f+gAgECAzJwfBKSlQE4vhpKMDCASAUIiHqamIyMlBK3V12PN4gFCIpNxcHMvMRP3ateE2ahQ0Bg0C6tdX9LPAMFJR+kTsazQ0NAAUTaKKc+LECQDAkCFDEBERgUuXLiElJQX29vbo3bs39PX1C+3v7e0NAOjRo0eRY+Vvu337tjThM4zsJScDly4BZ88i5d9/0S0rCwEAQgA4fNplP4B5AAYA2PfZQ78BkATgAIAGn7bdBrALgCOAXz7bdyOAVAA/8HiwMzYGTEyQbGCACG1tVDc3h2WNGoCJScHt0z6F/jY2Bkrx2VUogYCr7cpPuh48AN69AwD4AlgCoBqAY0ZGQIsWQFgYBn34gIigIKwLCkLnPXu44xgaAm3aFCRmbdoA1aop6KS+QCwGUlKKJFT5t7y4OMRERSH340fUz8zkticnYwuApwCmAsi/ZL0GoAe491zIZ0V8D+AqgEMAfvi0LQRAWwDWAMI/23ctgDMAdgCw/7QtDsB0AJovXiBx4UJoLFwIODiAvv0WvAEDAGfnypnwMpWakn8Lfll4eDiuX78OS0tLNG3atMT9Hz16BAC4e/cuZs+ejdzcXMl95ubmOHnyJDp37izZFhoaCgDFNj2amJjAzMxMss+X5ObmFionLS2txDgZpswiIiA6cwb3Dx5E8pMn6C8WAwCMAKSrqwNCITTHjgV0dIC8PKQ/fYokf39k1qrFJQ95eUBeHvzv3UO8QIC8evUANTUgLw9vkpPxV3Iy0tXV8YumJrevUIgtAN4B6EEEu+RkIDkZ1wEMBeACwOez8IYBiAKwAUDrT9veA7gMwE5HB73NzSXJWbq+PnTNzKBWrRpgYQFYWhbcatTgkhd5/tDGxxdOuh49ArKycOhTvFMBuPB4QKNG4Nnb4+r586hmbAxKTARPTQ1EhLvm5kj4+BE6P/4IhIUB/v7wTkvD0mvX0O/aNczOL6tePaBdu4LkrGlT+SSmmZlAdHSxt7SICDyPjARSU9EuLY1LxgAsBhAALtFs++kwtwD0AtAMXOKV7yyAmwDcatRASzs7wMwMpgBw7hwExsbA9OnceWloQPvIEWi9fg0aPhzo0AHQ0IBuXBxqb9wIKxMTYPVqbl91dTQ+ehQfg4JgMWwY0LUroKEB4+RkfLdyJYyysmBgaAjcvAm8fInBL18iY80arKhWDU6DBgEDBgBdunA1tgyj7EgF5eXlUadOnQgAHTp0qFSP0dLSIgCkpqZGv/76K0VERFBCQgJt3ryZNDU1ycjIiKKjoyX7d+/enQBQaGhoscezs7MjTU3Nr5bp7u5O4FoqCt1SU1NLf7IM8//EYqJnz4iWLSNydCQC6Myn95YdQOLGjYl+/53o0SPy8/WlyMhIEolEkocnJSXR8+fP6cOHD4UOe+3aNTpz5kyh92dQUBBt2LCB/vnnn4IdRSJyX7SIJo8bR5EBAUTPnxPdu0fH58+n6oaGNKh5c6KlS4lmzyYaM4bq6OoSALpvb09kbU2kr08nP8XrAhB9dnP8tP3fz7a9AmgBQEcAIjU1IisrolatSNy7N9G4cUQLFxJt3kx08iTRnTtEoaFE6eklP495eUSPHhFt2UL0/fdEdepQGkD7Afr987iMjGhkjRoEgDx/+IEoOZmIiHJycmj79u0UEBBAYrH400sjpvfv39OJEycoOzubK0cgoGVTphAAGlGnDlGDBpJjTwZoGUAJAJGuLlGnTkTz5xP98w/RZ99HxcrOJnr3jujuXe7cN27kHjtqFOV27kxxdesSGRpKyloO0LcA+X12bv99er6bfX6+hobUVUeHANDh5s2JfvyRaM4cCpg6lTTV1cnR3p7Ix4foxQui+Hjy2rePli1bRs+ePfvsqc2juLg4SklJKfl1kEZyMmUeOEDaamoEgJ5+dh4fdHUpqEcPEh8+TJSUJN84mEovNTVVbr/fKpeIiUQiGjVqFAGgCRMmlPpxGhoaBIC+/fbbIvf9+uuvBICWLVsm2SaLRCwnJ4dSU1Mlt4iICJaIMeUjFHI/fnPm0GVLSxoI0KH8Hx0ej9Lbt6cahob0w4ABlJGRoehoC/Hx8aFTp05R0mc/hrdv3qRB33xDv0+ZQuTnR/Tff0THj1OdatW4pO3777kEoEcPOm5tXWzS1g4gK4BufbYtDCAvgO4BRHp6RPb2RB06EA0aRDR1Kpe8zp9P5OJCsdradAagB589PvGzC6akzZu5ZEMkoosXL9KKFSsoMDCwXM/B27dvae/evXTt2jVuw8ePlHzqlKSsBAMDSQzXANoIUDBAZGNDNHQo0a+/Eo0ZQ9SjB1GTJkSmpkQAXQJoLUAvPzuHa5+O2fzz50tPj7p9Soi92rYlmjuXaP16erp6NdW2tKQ+Li5c4pebS0RE586do3379tGbN28k5yAWiyUJp7J5+fIlbdmwgcT//ks0eTKRlRXN//Q8zACI1NWJunXjku7wcEWHy6ggloh9IhaLady4cQSARo0aVegqvyRmZmYEgPbt21fkvrt37xZJ0gYPHkwA6NGjR188nrm5eZnil+cLyVRCWVlE58/Th6FDSWhmVqhmAwD1rV6daN8+org4IqIyfR6UVX5NSu6nhICIyM/Pj6ZNm0Z/rl1LFBVFFBBAdOkS1TQxIQDkN2QI0eDBRB070nELi2KTtqEAuQEU8Nm2xZ+exx+trIg8PYmuXiVKSaERI0bQvHnzKO7T8yovKSkptGHDBpo+fTqRSEQUEkK0fz+NqV+f8Cm+/FgDAXIGqPv/nVc3Ho9LrurV4xK2WbMoaM4cAkBW1aoRvXxJlJZGRER///03bd++nV69eiXX81IKIhFNHTyYtNTU6ETNmpLnKwGgMQCdtbMjsacn0dOnXA0zw5RAnr/fKtNHTCwWY/z48Thw4ABGjBgBLy8v8MvQV6RBgwZITEyEsbFxkfvyt2VnZ0u25fcNCw0NhaOjY6H9k5OTkZiYiPbt25f9RBjma5KSJJ3tceUKumZl4SaAuwA6mJgA33yDIc7OyIuKwnfDhnF9vD4py+dBWWloaKB69eqFtrVu3RqtW7cu2GBlBQDwCw5GdHQ0GjVqBOjqAgCM//sP3f/8Ey0aNQKmTQNiY4HYWNwdPx7RqakQ9+/P9Tdr0wbtBAI03bYNdQYOBJYskRy+oiZsNjIywqxZswo2ODgADg5oLxTi44UL6DR+PDcNxoMH0Hj+HH7HjsFUTw84c4Z7Dqys0H33blgEBcF63DiuTxQAB6EQ8b/9BjMzM4DHkxx+0KBBFXJeSoHPx9ZTp7A6I4MbzBURAZw7h4u7d8MrNBRP373Dt+7ugLs7UKcOcvr2hfbgwVy/NWUfQMJUPjJP7eRAJBLR2LFjCQANGzaMhEJhmY+xePHiIs2P+U6ePEkAaOLEiZJtV65cIQA0duzYIvsfP36cANCCBQvKFAOrEWOKJRCQYOdOutWyJW3k8QrVeIzU1SU+j0dbp0/n+jQx5eLt7U2HDx9W2c9eVlYWnTlz5os19EzpPHnyhGZOmEDbRo4k6tePSFubRADV+lSLGmZiQvTnn5ImWobJV6WbJkUiEY0ZM4YA0JAhQ0ggEHx1/4SEBAoJCaGEhIRC29+9e0daWlpUvXp1ioyMlGxPS0ujFi1aEAC6fv26ZLtAICA7OzvS0tKix48fF9q/cePGpK6uXuYqfpaIMUXcukXUpAmFf2om4wMU17Ah0eLFRAEB9P7duyLvZYZhZCQjg578+ScBIEMej3LzL4Lq1ye6cIE1WzIS8vz9VvqZ9T08PODp6Ql9fX3MnDmz2DnDBgwYgBafmmjy93d3d4eHh0eh/bZs2YIZM2agWrVq+O6776ClpYVLly7hw4cPmDhxInbt2lVo/1u3bqFnz57Q0tLCiBEjYGhoKFniaPny5Vi0aFGZzoXNrM9IREQgY+ZM6J85w/1taopeZmawbNIEHn/+idq1ays0PIapSsLCwvAiKAi9Y2OBRYuA+HgsAdC9dWu4eHkBjRopOkRGweT6+y3z1E7GRo8eXewUEJ/fDhw4INk/f8oId3f3Yo93/vx5cnFxIX19fdLW1iZHR0favXv3F8v39fWlXr16kZGREeno6JCTkxMdOXKkXOfCasQYys6mtN9/p5/V1KgWQGk8HtGUKUQfPyo6MoZhiIhSUshnxAhJDfUHPp9o2jSixERFR8YoUJWuEatMWI1YFUYEnD8PzJ6N3Pfv0RjAWwCHV6zAqDIuHs8wjHzFx8dj8YwZUPfzw7b377mNJiYgd3fwpkwBPq3qwlQdVXrR78qEJWJV1KtXeDZuHJrevw8eANSsCe+xY4EuXdC5mHVPGYZRDkQE3q1bwKxZSAgKgiuAX2vUwA9794Lfp4+iw2MqkDx/v1V/vDvDKKu0NNDcufipYUM0v38f59XVgQULgJcv0XnZMpaEMYyS4/F43LQggYHY2KsXQgBsiokB9e0L9OkDvHyp6BCZSoAlYgwja2IxcPgw0KABeH/+CUsi8AA8/vlnYOVK4P8WmGcYRsmpq2PJ2bNY4+GBzcOGQU1dHfj3X4ibNEHkTz8BycmKjpBRYaxpsgKxpskqIDAQl0aORJOXL2ELAHXrImP1aryqXbvIxMAMw6io16+BuXPhdeECfgbgqaOD+evWAZMmsQlhKynWNMkwyi4xEZg0CYsdHfHNy5eYq6YGrFoFBAdDf9AgloQxTGVSvz5w/jz+69wZOQD42dncSg4tWgDXrik6OkbFsESMYaQhFALbtnFfzLt3YwgAbTU11JkwAeL58wEtLYWGJxKLEJsRC1bxzai6tNw03Hp/C3/c/wO7Hu2Cf5Q/coQ5Co3p2M2bOH/mDKZv2ACYmgLPn+Nljx546OLC1ZoxTCmwpskKxJomKxfy9sZfP/6I3IgIjAWA5s2BLVvwsVEjVKtWrcLjEYlFeJn4EgExAQiMCURATAAexzxGpiATJtomcK7pDOeazmhTsw3a1GoDM12zCo+RYUpDIBIgKD4IflF+8I3yhV+UH0ISQkAo/HOlzldHY/PGaFWjleTW3KI59DT1Kj7opCSQhwd6bNmC6wA28/mYPmsWsHgxUMwax4xqYdNXVBIsEaskIiOBefNw5vhxDARgyOPh9apVsJg7F1BTq5AQhGIhQhJCCiVdT2KfIEuQVepj2JnYcUnZp8SspWVLaKkrtgaPqXqICO+S30kSLr8oPzyOfVxsbZetkS2crJyQnpeOwJhAJGYlFtmHz+PDwcyBS8wsW8HRyhEtLFvAUEv+37m5ubn4ecQI/HXuHJ6LxbADADMzYPlyYPz4Cvt+YGSPJWKVBEvEVFx2NrBhA7BiBZCVBRGPhy6Wlug5bhzm/P47tLW15VKsQCTAi4QXCIgJQEB0AAJjA/E09imyhdlF9tXT0EPLGi3hWMMRjjUc0apGK9ib2uN5/HP4Rvlyt0hfvPr4qshjNfgaaGHZQpKYtanZBnVN63JD+BlGRhIyEyQJl180929SdlKR/Yy1jblaXCtnSW2uhb6F5H4iQmRapORCJDAmEIExgYjJiCm23Hqm9dCqRivJ56JljZYw1TGVyznGxcXBIjAQmDMHePkSGwEYWVjg+99+g9ZPPwEGBnIpl5EflohVEiwRU1Hh4cjYtAmbd+yAT3Y2LgPgd+wIbN4MatFCpolKnigPz+OfF0m6ckW5RfY10DQoknTVr1YfavySr7qTs5PhH+0P30hfSYJWXO2CqY5pQXNmzTZwrumMaroV3+zKqKYsQRYCYwIliZdvlC8+pHwosp+mmiZaWraUvNecazqX+yIgJj1GkpQFxnL/hqeGF7tvbePaks9OfpJmrmde5jK/SCDAh+XL4bB0KXIB3AHQ0dAQGDcOmDoVqFtXdmUxcsUSsUqCJWIqhAjw8QG2bAHOnEGWWIyaAFIAHJ82DcM2bwZklIDFpMdgw8MNuPXhFp7FPUOeKK/IPoZahpIfCscajnC0ckRd07rg82Qz3oaI8D7lfaHE7HHM42ITwLqmdSWJWVe7rmhkzhZEZgokZiVih/8OnA45jeD4YIhIVGQfBzOHQklXM4tm0FTTlGtMkuTs0+1t8tsi+/HAw7Amw7C402KZva+zs7Ox9Y8/cOfkSZzLzQUvNBQAcBGAabt2aOfuDl6PHjL7PmHkgyVilQRLxFRAVhYEhw7h9KpVCAgPx7r87V26YHvdujBo3x4jRo6EugzmCorPjMeau2uw/dH2Qv1hjLWNiyRddiZ2Mku6SitPlIensU8LNWmGJoUW2a+bXTfMaTsHPev2rPAYGeXxJukNNjzYgANPDhRqNq+hXwNtarWRNDE6WTnBSNtIgZFyUnJS8DjmcaGas5eJ3Ez5PPAwpPEQLO60GE2qN5FdoWIxcPUqRJs2oe6VK/gA4BSAwQ4O3PQXP/7Imi2VlNIkYnZ2dlIXOGvWLMyYMUPq46gilogpsQ8fgO3bgb178T45GfYACEDI0KFwWLwYaCK7L+OPWR/xx/0/sMVvCzIFmQCAdrXaYUabGXCu6Yw6xnWUtl9WUnYS18QU6Yv7kfdx/d11iEkMAGho1hCz287GqGajoKOho+BImYryMPIh1t1fhzMhZySjGltatsSstrPQpU4X1DSoqbTv5//3NPYplvksw+mQ05JtgxsNxuJOi9HMopnMyklJScEv48fjvytXEMrjQScjAwDwQk8PRiNGoOavv7JmSyWjNIkYn8+HkZERjMs5FDc8PBzu7u5YsmRJuR6v6lgipmSIgFu34Ld0KV75+OCH/I9CnTqYaGGBmq6umDp3LszMZDPNQ0pOCtY/WI+NDzciPS8dAOBk5YRlbsvQ076nyvxYfe5Dygds9t2MvYF7JedkpmuGKU5TMKX1lEKdq5nKQyQW4cLrC/jj/h+4F3FPsr1PvT6Y224uOtfurJLv53xBcUFY5rMMp16ckmz7zuE7LHFdghaWLWRWTl5eHjRzcoCDB4GtW9H99Wt4A/ACMLJvX2DGDKB7d9ZsqQSUKhHz8PAodyIl7eNVHUvElERmJnDkCLBlC+4/f44OAPQBRLm5wXD2bG4xXxkOM0/PTccm303488GfSMlJAQA0t2iOpW5L0a9+P5X+wcqXlpuGfYH7sMl3E8JSwwBwHbBHNh2J2W1no6lFUwVHyMhCtiAbB58exPoH6yXN1Bp8DYxqNgq/tPsFjas3VnCEshUcH8wlZM9PSWr7vm3wLdxd3dGyRkuZlpWTlYWebdrgbnAw3gHcEmkAkurVg+6UKdBmoy0ViiVilQRLxBTs3TtEr12LiKNH0eZTUwDp6qKljg6adeiA1Tt2wMrKSmbFZeZlYqvfVqy7vw4fsz8CABqbN4ZnZ0981/C7StmfSigW4kzIGax/uB4PIx9Ktne364457eaobM1fVZeQmYBt/tuwzX+bZHStsbYxfnb6GdOdp6OGQQ0FRyhfLxJeYJnPMpwIPiFJyPrV7wd3V3c4Wsl2+bIPHz6gdl4et2LHgQOYkp6OUwA26+hgxKRJbLSlgihNIvb27VuYmprCxMSkXIVJ+3hVxxIxBSACbtwANm/Gfxcu4BsA9gBe2NmBP306MHYshHp6Mul8ny9bkI2dj3Zi9b3ViM+MBwDUr1YfHq4eGNp4aKmml6gMHkQ8wIaHG3A65LSkH1kj80aY3XY2RjYdyfqRqYDXH19j/YP1OPj0oGRASW3j2pjddjbGtRwHfU19BUdYsUISQrD8znIcDz4ueU/3rdcX7q7uaF2ztczLE6ekoImDA0Li4nATgBsA8HgQ9uoFtZkz2WjLCqQ0iRgjHZaIVaCMDOTs24ekrVth9eYNACAdQC11dTSrXx+nb9xAdUtLmRaZK8zFnsA9WHlnpWRSSTsTO7i7uuP7pt9DnS+7ZE+VFNePzFzXHD87/cz6kSkhIsL9iPv448EfOPfynKQGyMnKCfPaz8PAhgOr7Hs536vEV1h+ZzmOBR2TJGS96/aGu6s72tRqI9OyhEIhrl65gt5qauBt2QL8+y82AjgEwKNmTfRfuBAYOxbQYRc28sQSsUqCJWIV4ONHYNUqXNq+HWOzs9EGwAV9fe6LaupUxBgaokYN2TajCEQCHHhyAMt9liMiLQIAYGNkg8WdFmN089HQUNOQaXmq6kv9yEY1HYXZ7WbLdpoApsxEYhHOvjyLPx78UahZuV/9fpjbfi5cbFxYs/L/ef3xNVbcWYGjz45K5kvrad8T7q7uaGfdTk6Fvkaz9u0R9PEjdgKYBAC2tsDatcCQIayGTE5YIlZJsERMjgQCbvoJDw8gJQWvATQAYGtqiufBwdCTcfIFcP2hjjw7gqW3l+J9ynsAQE2Dmljksgg/tfpJrhNUqrIv9SPrYd8Dc9rOQQ/7HuwHvwJl5mXC64kX1j9cj3fJ7wBwCfKPzX7EnHZz0NC8oYIjVH5vkt5gxZ0VOPz0sCQh627XHe6u7uhg00Hm5X38+BEHduzAzzo60Nu4EYiMxBMAkY0aoe/Bg+A5Ocm8zKpOpRIxgUCAc+fO4eHDh4iLiwMAWFhYoF27dujXrx80NavujxNLxOSACKILF7B34kTkxcVhOgA0bw6sXIk7enpo16GDTPt/AVzNwfHg4/C87SkZOWahZ4EFHRdgktMkaKvLZ83Jyuhr/ch+bP4jS2blKEeYg9V3V2OL3xbJWo+mOqaY4jQFU52nwlJftk33VcHbpLdYeWclDj49KEnIutbpCndXd7jYusin0MxM0Lp16LRsGe6KxVgDYP6YMdyauDIcfFTVqUwiFhoail69eiEmJgZt2rRB9erVAQDx8fHw9fWFlZUV/v33X9SrV09WRaoUlojJ2PPnwJw5OH/1Kr4FoAvg9dq1qDlnjkynn8gnJjH+fvE3PLw9EJIYAoCbM+vXDr9iSusp0NXQlXmZVcX75PfY4relUD+yNjXb4PTQ06hpWFPB0VU+kWmRGHRyEPyi/ABwfRnntJ2DMS3GQE9TT8HRfZlYDNy/D2RlARYWQPXqgLk5IONrLam9T36PlXdWwuupF4RiIQCgg3UHjG0xFkMaD4Ghlmy///Py8rBk9mzs3rsXz/LyUAsA9PSAhQuB2bNZ/zEZUJlErGvXrjAxMYGXlxf09QuPpsnIyMDYsWORnJyM69evy6pIlcISMRlJTIRw8WKo79kDiEQgDQ30t7FBtwkTMGXOHGhoyLZPVujHUBx8ehCHnx2WLB5som2Cue3nYrrzdBhoKWZun8xMIDYWiIsruBX3d1YW0Ls3MH480K6dcnchSc1Jxb7H+7DcZzmSc5JhoWeBU0NOya82oQq6G34Xg08ORlxmHEy0TbC973YMaTREqUfzhoUBBw5wt/Bi1u+uVq0gMfv83+K26Vbg9dKHlA9YdWcVDjw5AIFYAADQUdfBoEaDMLr5aLjVdpPp856RkQH94GBg1izA1xeeAAxNTDB12zZoDh+u3B9+JacyiZiuri4ePXqERo2KXyw1ODgYzs7OyMrKklWRKoUlYlLKy0Pan3/C08MDt/Ly4AdAfeBArpOqvb1Mi0rJScGJ4BM4+PQgHkQ+kGw31jbGzDYzMbvtbLmsl5eRUXJilf//zMyyH79RIy4h++EHQEYLBsjFu+R3+O7Ed3gW9wzqfHWs77Ee05ynsb5jUiAi7ArYhen/TodQLETT6k1xdvhZ2JlIv3SdPOTmAufOAfv2AdeucTPRAICxMWBtzX0GEhO5WrKy0NMrPkGrXh2wsQF69JB9BVJUWhSOPDsCr6dekvUsAcDa0Bo/NPsBo1uMRv1q9WVXoFiMd5s2wWHOHAgAXAPQzcUF2LgRaNVKduVUISqTiNnY2GDdunUYNmxYsfefOHECc+fORUREhKyKVCksESsnIuDSJeCXX5D8+jXqAkgCcHbZMnz7++8yK0YoFuLa22s4+PQgzr48i1xRLgCAz+Ojp31PjG4+Gv0b9Jfp/Fe5ucDx49w4g+fPy55c6egAlpYFV//5t8+35eUBhw8DJ09ytWMAoKEBDBjAJWXdugF8JZxbNjMvExMuTMBfwX8BAH5s/iN29t3J5h8rh1xhLqZdnoa9j/cCAIY2Hor9/fcrZTNkUBCXfB05wg2CztelC/DTT8B33xUkSiIRt098PJeYlfRvTk7J5Zuacp+Ln38GateW7bkREfyi/HDw6UH8FfyXZKUNgFtvdkyLMRjaeCiMtY2lLkskEsFr1y7c2b0bXq9fA9nZAI+H7FGjoLNmDSCHAUyVmcokYhs2bMCiRYswbdo0uLm5FeojduvWLWzfvh0rVqzAzJkzZVWkSmGJWDkEB+P5xIlo/OBTrVT16jgxYACM+vdHr759ZVNEfDAOPjmII0FHEJsRK9ne2LwxxrQYg5FNR8p85vC4OGDHDu4WH1/4Pl3d4hOq4rbp65e+tSE1lUv69u4FHj0q2G5rC4wbx83wYW0tu3OUBSLCxocbMe/aPIhIhFY1WuGfof/A1ti25AczAIDo9GgMOjkIDyMfggceVnVdhfkd5itV7WJaGvfe3LcP8PMr2F6zJjBmDPf+tJOy4o6ooMb5S4mavz/XDApwFyf9+gHTpgFdu8q+VS9HmIPzr87j4NODuPLmimSwira6NgY4DMCY5mPQza6bbJouIyKABQuQffQomgLooa6OVQsXwmjBAkCbDS4qDbn+fpOMHT16lFq3bk3q6urE4/GIx+ORuro6tW7dmo4ePSrr4lRKamoqAaDU1FRFh6L84uNJMGkS9QMIAPmpqxP9+iuRjJ67hMwE2vxwM7Xa1YrgAcmt2ppqNP3ydHoU9YjEYrFMyvpcYCDR6NFEmppE3E8DUa1aRKtXE716RZSeLvMii/X4MdG0aUTGxgVx8PlEvXsTnT5NlJdXMXGU1s13N8lsrZnkNbrx7oaiQ1IJ98Pvk+UflgQPkPFqY7oSekXRIUmIxUR37hCNGUOkq1vwPlRXJxo4kOjSJSKhsGJjEgqJzp0j6t69IB6AyMGBaOtWorQ0+ZQbnRZN6+6to8bbGhf6PrL604p+vfYrvYh/IZNyTq1YQQCoFkCZAFHt2kSnTnEvBvNV8vz9lnkili8vL4+io6MpOjqa8pTtW11BWCJWCrm5RH/+SWRkRATQDwCp83i0zdNT+kMLc+lMyBkacHwAaSzVkHzZqS9VpwHHB9CZkDOUK8yV/hz+j1BI9M8/RJ06Ff5yb9uW6PhxxSY9WVlER44Qde5cOLbq1YnmzSN6+VJxsf2/sJQwctzlSPAA8T35tO7eOrkky5XF7ke7Je/zxtsaU+jHUEWHREREMTFEa9YQNWhQ+D3XsCHRH38QxcUpOkJOSAjR9OlEBgYFMRoYcBcwISHyKVMsFpN/lD9NuzSNTNeYFkrKnPc40za/bfQx66NUZXjfvEmX584lqllTcmIPW7Qg8aNHMjqLykklEzGmKJaIfYVYTIIzZ2iHuTkl5X/rtWhB0X//TSFSfOuJxWIKiA6gGZdnSGpU8m+tdrWiTQ83UXxGvAxPpEBKCtH69dxF5+dX+yNGED18KJcipfL6NdFvvxFZWhb+gXRxITp4kCgzU9EREmXlZdGYs2Mkr+GwU8MoIzdD0WEplVxhLk26MEnyHA06MYjScyuoqvULBAKiCxeIBgwgUlMreG/p6RGNG0d0757yVsqkpnK1YQ4OhT8X3btztWfyqrXLEeTQ6Renqd+xfqTmqSZ5PTWXadKQk0Po4quLJBAJyl9ARgaRuztd19QkANQbIOGYMVymzBSh1InYixcvaMmSJfTTTz+Rp6cnnTt3jsLDw2URW6XDErEvePaMqFs3GvSpGXK2jg7R3r1SfcPlV/U32d6kUPJl+Yclzf1vLgXFBcnwBAp7/Zq7ktbXL/jSNjUlWrCAKCJCbsXKTF4e0dmzRN98wzVX5p+DoSHRzz8TBQQoNj6xWEzb/LaR+lJ1ggeo6fam9ObjG8UGpSRi0mOow74OBA8Qz4NHK3xWKLTWMDSUe99bWRVOYtq14z7i8mrqkwexmOjaNaJvvy38ubC15Wr4EhPlV3Zseiytv7+emu1oVuz3WXhK+X9zty1fTpp8Pk3PPyF9faKVK4mys2V4BqpPaROxe/fuka6urqQvGJ/Pl9yqVatG3bp1o3nz5tGxY8dkFa9KY4nY/4mPJ5o8WfKt9p+6OlXT0aHdmzaV+5D+Uf7U71g/4nvyJV9WWsu0aNipYXT59WXpriC/Qiwmun6dqF8/Ih6v4Eu6USOi3buVozapPCIjiZYvJ6pTp/APacuWRNu2ESUnKy62O2F3CvV/uvz6suKCUQIPIx6S1Z9WBA+Q0SojuvT6ksJiuXWraHO3mRnRnDlEz58rLCyZef+eaP587gIr//y0tbnavcBA+Zb9OOYxzfx3ZqEafq1lWvTLf79QQmZCuY759u1b+njlCpGzMxFA0QCtNzGh3JMnZRy96lLaRKxHjx6ko6NDFy9epKioKOLxeOTm5kZdu3YlNTU14vP5kgSNYYmYhFhMSRs30gxNTTqV/y02aBDR27eUXs7e6iEJITT45OBCV4vt9rajnf47KSkrScYnUCAri7uyb9Kk8I9O375EV68qb3NLWYlEXKI5fHjhgQbVqnHnqShRaVHUbm87SQ3Q8tvLSSQWKS4gBdkXuI80l2kSPEANtzak14mvFRKHWMwNPMmvMcofAPL331z3z8omK4to/37uwuTzz3/79kTHjsn3nHOFuXQ25Cx1OtBJ8p1nsNKAPL09KS2nnFWNIhHR4cM0TleXAFAfgGjqVKKcHNkGr4KUNhGzsLCgYcOGSf7m8Xjk+alT9ZMnT6hevXo0ZMgQ2rVrl3RRVhIsESOuunvcOJr/qRnSVkODcq9dK/fhwlLCaNzZcZIaMJ4Hj0b9M4pCEuTUm/aTqCiiRYu4q/zP+7tMncqNfqzMEhOJNm4kql+/4Md2xQruO1wRcgQ5hfpEDTg+gFJzqsZnLE+YR1MvTS107uX+EZZSWhp3PZX/eRgzhqiq9FIRi7l+biNGcP1A858DS0sid3fu+0J+ZYvp39B/qcXOFpL3gflac9r4YCPlCMqXQB3YuZNsDA3pUf6JODsThYXJOHLVorSJmK6uLi1cuFDyN5/PpyVLlkj+fv36Nenp6dGtW7ekKabSqPKJWFgYkZMTEUDZPB4Na9aMrv/3X7kOFZ8RT7P+nSWpBYAHqP9f/elZ7DMZB12Ynx/RyJGFv2xtbbnRXopsplOE7Gyin34qeB7691fsc7A3YK/k/eCw1YFeJijRkE85iE2PJZf9LpL3v6e3p8JqA1++5EY9AkQaGkS7dlWe2uCyio4m8vAoPOhFXZ0bqHDkCDeIRx5EYhEdDzpO9TbXk7wnbDbY0IHHB0goKnt/W0H+CItPc9w8MzKinPPn5RC5alDaRMzOzo4mT54s+dvQ0JBmzJhRaJ8hQ4ZQz549pSmm0qjKiVjulSv0V/44cFPTcrdnpeakkvstd9JfqS/5snE94Er3w+/LOGKOWMzNubV4MVHjxkVHE/79NzcirCrbs4dIS4t7TurW5cZeKIpvpC/V/LOmpJnmbMhZxQUjR36RflRrfS3JeZ57eU5hsZw5UzDFg5UV0YMHCgtFqeTmctPTdOhQ+HtDQ4Nrrt2zh+smK2t5wjza9WiXpL8gPECNtjWif178U76BG+/eUXiTJlQNoNYAxfzyS8VP8KYElDYRGzhwIHXt2lXyd8uWLcnNza3QPr/++isZGxtLU0ylUSUTMbGYBKtXU4dPTZH7bW25nq5llC3Ipj/v/0nV1lQrNP3EldArMh8VJhJxzQy//FK0k7qGBtEPPxCxKXcK8/cnsrHhniMdHe7KX1Fi02ML9ZtZfHNxpeo35vXYi7SWaRE8QA22NJB7M/yXCIVECxcWfDY6dSKKjVVIKErv6VOi338vqDXMv/H5RK6uRJs2yb4ZNysvi9beXUsmq00KzUVWnsmQfa5dI1MtLXIEKBsg6tGDKKF8AwNUldImYrt37yYNDQ1K/tQesWDBAlJXV6dnn10Su7i4ULVq1aQKsrKocolYejrRkCFEAHkCZKypSZf++adMhxCIBLQnYI/k6h8eoPpb6tPJ4JMy/XHNy+Mq6SZPLjqPlo4O0XffER0+XPWaH8siIYH7fs5/3qZNU1wH7TxhHs38d6bkPdPnaB9Kzk5WTDAykifMoxmXZ0jO6Ztj31BKtpzauUrw8SNRz54Fr/WsWcq3GoOyevGC61Pp6Fj4eya/K9bq1dwUOLKSnJ1Mi24sIt0VupL3TvdD3ck/yr9Mx/nw4QO9//NP7gsRILK2JtF9+bREKCOlTcSIiMLCwiSBxcXFkYmJCRkbG9P3339Pzs7OxOfzaejQoVIHWhlUpURM/PIlZefPgKiuTqItWyiiDJd8IrGITgafpPpb6ku+PGqtr0V7AvbIbAqKrCxuvqwffyQyMSn8hWhkRDRqFDcjvqpOPaEIQiF35Z//PLZrx02BoSiHnx4m7eXaBA9Q3c115Tp/nDzFZ8ST6wFXyWdhyc0lCqvle/y4oKZYR4eoiq9cJ5X377lJnzt2LDztDUDUtCnX0f/pU9n0t4tJj6Fpl6YVWlVk0IlBZa9RffaMqF492g9QNx6P4lesqBIdApU6Eft/vr6+VK9ePcncYm3btqXo6GhZF6OSqkoilnnyJI3S0KA+AIksLbl2vlISi8V0JfRKoTUgq62pRn/e/5OyBdJPMJiSwg0rHzy48Pp2ALesz8SJRFeuVM6h9hXp/HnJKlVUvTo3r5SiBEYHku0GW4IHSG+FHh16cqhcnZcVQSwW0413N8hmgw3BA6S/Up/OhJxRWDyHD3PzZQFEdnZcksDIRkwM0c6d3Iz9nw8Gyu97OX8+tyKHtKOT3yW9ox/++YF4HjzJcmHjzo6jsJTSj4pMj4ois08z8v8JcPPaVNRCuQqiVInY/v37Kb4UPQxDQ0PZDPv/p9InYkIh0eLF9AwgbYDUALpfhlE2DyIeUGevzpIETH+lPrnfcpd6KoL4eK5jbO/eXB+vz7/gbGy4ZhUfnyrZ/1SuQkO5q3qAW9bmjz8Ud+GckJlA3Q51k7y3avxRg+Zdnae0NWQfkj/QUu+lZL/JXhJzvc316Hm8YmZDzcvjVovI/9z07k2UJL/p+aq8jx+5ZcW+/bYg8c2/1azJNfvfvCndQKGguCD69q9vC00KO/vK7FIv+fY8OJimdexIovw1qxwcKsdsvV+gVIkYj8cjdXV16tixI/3xxx/0WpaN2ZVcpU7EkpK4b+dP3xbHe/Yk71LODxYUF0T9/+pfaC21Wf/OkmoNyPBwrgOsq2vh5UgArsPsokXcUj1VoEZdoTIzuSbe/Od+8GDFLWsjEAnI45ZHoQEfFbHmaGml56aT12MvcvNyKxSf3go9mnB+gsL6uEVHc01n+a/h4sXsoqUipacTnTzJVTp9vmwaPq1W8NNPXHNxed0Pv1+o2dtgpQF53PIo/Xx0d+8SWVmRCKDVGhqUtm9f+YNRYkqViD148IB+++03atSokWTWfAcHB1qwYAE9YOOWv6qyJmKiwEBaa2pK7/FpnY9Dh0r1uLdJb2nUP6OkqiIvFIeIG0rfrl3RTrCOjlwH2RcvynVoRgpiMbdocn5tpIODYl+HXGEunQk5QwOODyjUX0Z9qTr1/6s//f3873JPhFlWIrGIbr67SaPPjCa9FXqSWHgePOpysAsdenJIoYua37tHVKMG97oZGnJNzoziZGcTXbzILaVUrVrh77h+/Yh8fct33OK6hJitNaMtvltK1xcxLo6W2tkRAGoJkPDnnyvdbPxKlYh97s2bN7Ru3TpycXGRLGlkaWlJEydOpEuXLlFOJXshpFUpE7GjR2mBujr3AdTUpLxSfBMIRAJaeH1hoR/BwScHl3sYvlDIzdfz+TJDPB43z9eGDUQfPpTrsIyM3b9fsPizvj7RqVOKjohrstziu4Va725dqBbKdI0pTbk4hR5GPJTLotmhH0Np8c3Fkr5r+be6m+vS8tvLy30xIitiMbeWaH7y3KhR5V8xQtUIBEQ3bnA1ZZ/X+vfoQXTnTvmOWdwgqV5HelFcRlyJj31w9y7VMjAgr/xAKtls/EqbiH0uMTGR9u/fT99++y3p6ekRn88nfX19GjhwIB06dIg+fvwoq6JUVqVKxPLyiGbOJAIoDKBa2tq0txSLdUemRhaaDbw8w6jzCQRcP4oGDQq+hAwMuLmNYmLKdUhGzmJjCy8G/csvyjMh7vP45/TrtV8lE8Lm3xpsaUArfFZQeIp0fV5TslNoT8Ae6ri/Y6HjG64ypInnJ9K98HtySfrKKiuLaPTogtdoyJBK3w9b5b16xb1m+d21AO5zduNG+bpfCEQC2uK7RTLi2PIPS7r+9nqJj0tLS+Oq7D4NQ482MaHsSlKNqhKJ2Oeys7Pp3LlzNG7cOLKwsJD0K3N1dZVHcSqj0iRiMTEU7uxc8IlftIiySvFN/d+b/8h8rbmkH8LxoOPlKj43l2j3bm7UVn4IJiZEnp6sA7EqEAiI5s0reO1cXZVrIlChSEhX31ylUf+MIp3lOoWaC7se7EoHnxyk9NzSZSZCkZD+e/MffX/6e8mPWn4zfM/DPemvoL8oKy9LzmdUeu/fFyxgzecTrVvH+lGqkrdviSZMKDwoqX17on//Ld/rGBQXRI22NZK8/xfdWFS66YPev6fsli2pFUCOAH2YMUPlOxaqXCL2ObFYTPfu3aN58+ZRgwYN5F2cUqsMiZjwzh2ar69P2gD56+pynbJKeoxISL/f+F3SF6z5jub0OrHsgzyys7m+RtbWBV8y5ubcBIiK6gDOlN/ffxd0PrayKtMsJxUmLSeN9gfuLzSaN78D/egzo+nGuxvF9qEJSQih3679VqR2reHWhrTm7hqKSpPjKtDldPVqQb8jMzOuNoVRTWFhRFOnFiw9BnDL/J47V/aELDMvkyacnyB5D3fY16FUTedP/fzIVEuLzAAKz28zVeHZ+FU6EWMKqHQiJhYT7dxJInV16g9uuaLVv/xS4sNi0mMKjQKbeH5imWsAMjKI/vyz8Iz3NWpw/b8yFNePmZGBkJCCZV/U1Ym2bFHeGpj3ye9p2e1lVHdz3ULJlfV6a1p4fSEFRAfQdr/t1GZPm0L3m6w2oamXppJfpJ9SND3+P7GYaNWqgn5GTk6VqmtPlRYdTTR7dsFk+ABRs2bcKMyyzkd2POg4Ga4yJHiAjFcb0z8vSl4lJSwsjO7+/ntBALVqqexipCwRqyRUNhHLzuaG6Xz6JKf070/njpfcrHjj3Q2yWGchqUE4+qxsU3CnphKtXMldned/idjYEG3fzoXEVA5paZKVsAggGjlSuRNssVhM98Lv0aQLk8holVGhpCv/puapRv2O9avQEZjlkZZGNHBgwXM/diz7bFVGcXFEv/1WePqLhg25NWHL0kfzbdJbct7jLHmfT7k4pXQTbQcFEdWvT/4A9eHxKE4FZ+NX6UTs7Nmz5OnpKe9iVIJKJmJhYbTf1paW53caWb26xA+QUCQkT29P4nvyCR6gJtublGlEZFISt7SHsXHBl4a9PdHevWzG+8pKLOZqPfM7GzdpQvTunaKjKlm2IJtOBp+kvkf7kvpSdWq2oxmtv7+eYtOVqNPbF8TGFtRGamhws7qr2G8jU0YfPxItWVKw6gU+zdq/f3/p1wrNFebSvKvzJMlYsx3NSvX9LkpOpsaGhgSAJgBE/ftzVXYqQqUTsTFjxhCfz5d3MSpB5RKxU6fo4acPDgC6t2FDiQ+Jy4grNIP5uLPjKDOvdIs1xsdzV20GBgVfEg4O3LIqyjKyjpGv27eJLCy41755c1Y7Iy8iUcEC7VZWKttaxJRTSgrR8uWF5yKrXZtLxks769S/of9KBl/prtCl/YH7S2x6fx4cTN80akTJ+Ws4GRtzQ99V4AqAJWKVhMokYnFxhdqKZpibk+fs2SQqoVPB7Q+3qcYfNQgeIJ3lOuT12KtUxUVHE82ZU3jtx/x+DCo+0IYph/DwguboGTMUHU3l9Oef3POrrV2pV6VhSpCeTrR2LbcebP53b61aRJs3c9OYlCQ6LZq6HuwqufD+/vT3pZuRPyiIm2UboKUAXXd2JoqMlP6E5Eiev988IiKUwaFDh8qyO/bs2YP79+9DJBKV6XGVUVpaGoyMjJCamgpDQ0NFh1MUEbIOH8bayZMxJzsbhmpqwIIFoEWLwNPW/uLDxCTGmrtr8Put3yEmMRqaNcSpIafQuHrjrxYXEQGsWQPs3Qvk5nLbnJyAxYuBb74B+HxZnhyjSi5fBvr25f5/4QL3fmBk4/FjoE0bQCAAduwAJk9WdESMomVlAXv2AGvXAtHR3DYLC2DhQmDqVEBN7cuPFYlFWHNvDZbcWgIRiVDXtC6ODzoORyvHrxcqFMJn6lS47t4NHoAQfX002LwZGDMG4PFkdWoyI9ff7zJnbp+WNSrtLX9/RslrxGJiiL77jnp+aoacZGrKLcZYgoTMBOp9pLfkiuiHf34ocUmW3Fyi+fMLz3XToQPRlSsqUUPNVJBP8wWTmRlRlPLN9qCSMjIKJkAeMIB93pjCcnKIduwgsrUtPM9faVYnuRt2l2w22BA8QBpLNWj9/fUlNlWmpaXRlOHD6efPq+R69uSqxZWMUtWIaWtrw8rKCpMmTSrV/qdOncLjx49ZjRiUtEaMCDh2DJgxA0hKgjefjx/19bHz0CH0+fbbrz70Xvg9DD89HJFpkdBW18bW3lsxruU48L5yNfPqFfD990BgIPd3ly7A778DnTsr5UUQo0C5uUDbtsCTJ0DXrsDVq6yWVFoTJnA10FZWwLNnQLVqio6IUUYCAfc+mTcPyMwEDA2BbduAkSO//j2dnJ2Mn87/hDMvzwAA+tbriwPfHoC5nvlXyyOBALyNG4HFi5Gam4ulGhpYvG4djGfMUJofBqWqEXNycqLq1auXen/WR6yA0tWIRUfTrfbt6Vr+lUiLFkRPnpS4RqhYLKZ199aRmqcawQNUf0t9ehr7tITHEO3bV9APrFo1orNnZXkyTGUUElLwnlm9WtHRqLZTpwrWYb15U9HRMKrgzRuitm0LKquGDSt59RKxWEzb/LaR1jItggfI6k8ruvX+VukKDAmhn6pXJwDkBhB17640iwUrVWf9SZMmEZ/Pp/BSVh2yRKyA0iRiYjGRlxf9o6tLAMgKoKSFC0s1fvlj1kfqd6yfpClyxN8jSuycmZRENHRowYe5Sxel75fJKJG9ewsmfC3FmvJMMcLCCqaDWbBA0dEwqkQgIFq6tGBqmZo1ia6XvOwkPYl5Qg22NJAsj7Tk5pJSLY/kc+sW1Tc3pzuamlyB+vpce2lZZ6CVMaVKxI4cOUK1a9emG6Vc/2Lv3r00ZsyYMgdWGSlFIhYRQdSnDxFAWQA5aGvTpCFDuMVaS/Aw4iHZbrAleIC0lmnRTv+dJfYBuHOHm4Q1/4d09WqFf54YFSMWFyTydnbcRL9M6QmFRC4u3PPn7Fz6+aIY5nO+vkT16hVcUM+eXfL0Mhm5GTTu7DjJhbvLfhcKTym5EkcgEBC9fk3UsSMRQJcAOt+0qUInF1SqRKyiRUZG0oYNG6h79+5kbW1NGhoaZGFhQQMHDqSHDx+W+ji3bt2SzIdV3O1BMRPp2NrafnH/SZMmlflcFJqIicWUunUr7dbW5j5FmppEq1ZRenJyKR4qpg0PNpDGUg2CB8h+kz0FRgd+9TECATcpa/6yKfb2RH5+sjkVpupJTi7oQDxypKKjUS1LlxZULLx5o+hoGFWWkUE0eXJBMtakCdGTJyU/7uizo2Sw0oDgATJdY0rnXp4rXYEiESWuXEkWn353j2tpceugKeBqXqk661e03377DWvWrIG9vT1cXV1RvXp1hIaG4uzZsyAi/PXXXxg6dGiJx/H29oabmxtcXV3RuXPnIvePHz8etWrVKrStdu3aSElJwaxZs4rs7+TkhG/KOKZeYZ31w8OR89NPaHT9Ot4DOFW3LgafOwc0alTiQ1NyUjDu3DhJ58vBjQZjb7+9MNI2+uJjwsK4Tp337nF///gjsHUrYGAgi5Nhqqr794FOnQCRCDh0CPjhB0VHpPzYc8bIw8WLwE8/AfHxgKYmsGIFMGfO1wfTvEl6g+F/D0dATAAAYHzL8VjbfS1MdEy+WlZ2djY8Zs3Cv0ePwi8zE9oA96bevx+wt5fdSZVAqTrrV7TTp0+Tj49Pke0+Pj6koaFBpqamJXYuJyqoEXN3dy912ba2tmRra1uGaL+uwmvExGKiXbskU9UvUlMjO1NTuuPtXaqH33x3k+psrEPwAGku06QtvltKbIo8frxg+QxDQ6KjZVtekmG+6vPandevFR2NcktJYbWIjPzExXGrFOXXjnXuXPJi8bnCXJp9ZbakqdJinQUdDzpe4u8KEVFOVhbRtm1EenokBmiDhgYlLF9eYbVjStM0OWLECDp9+nS5C5P28f+vR48eBID8/f1L3LfKJWLv39OlFi0oOv9T0q4d5Tx9ShmlWE05MTORxpwdI/mw1NlYh/yjvv4cp6dzCwbnF9e2rWqsFcioFqGQqFMn7j3m6MjWHv0SsZho+HDueapTh/WrY+RDLCbavZtIT497rxkZle7i2+eDDzlsdZD8xvQ92pc+JJdydOS7d3S8cWMCQDUBymzblujVK6nOozTk+ftdpll5jh8/juDg4HLXvkn7+P+noaEBAFBXVy/1Y0JDQ7F582asXr0af/31FxITE7+6f25uLg4ePIiVK1dix44dePr0qVQxy51YDOzYgRUNGqDvkyeYzOeD/vwTuHMHWs2aQU9P74sPJSIcfXYUDtsc4PXECzzwMMVpCh5PegwnK6cvPi4gAGjVCjhwgJvy5fffAR8foE4deZwgU5WpqQFHjgAmJtz77vffFR2Rcjp0CDh+nHu+jh3j5oFiGFnj8bi56Z484VZrSE3luqWMGAEkJ3/5cS62Lngy6Qk8XD2gqaaJS6GX0Hh7Y2x8uBEicQlzjtapg7peXmhiZYUJmprQffgQaN4c+PNPrg1eBZWpjxifz8eAAQMwYMCAchU2ZswYeHh4YMmSJeV6/OfCw8NRv359mJiYIDIyEmpfW4MBBX3E/p+Ojg48PT0xb968IvfVrl0bYWFhRbb36tULhw8fhpmZWZlilnsfsXfvuIZ7b28EAXDm8zF1zBis3rWrxGT1XfI7/HzpZ1x9exUA0Ni8Mfb024N21u2++BixmHvvL1rETQBYqxb3I+nqKsuTYpiizpwBBg7k/v/ff0CPHoqNR5mEhgItW3ITcS5fzn0+GUbehEKur9iyZVw+VKsWcPAgN2n314QkhGDSxUm4E34HAOBYwxF7+u1Byxotv/q43Nxc8CMioDFlCnDtGqIA3K1bF0PPnwevYUMZnVUBpekjxuPxyrzE0edLHfF4PPL09JS6Gi8vL486depEAOjQoUOlekxwcDCtW7eOQkJCKDMzk6KioujIkSNUs2ZNAkA7d+4s8hhPT0/y9vamhIQESktLo4cPH1Lv3r0JALVr167Edu2cnBxKTU2V3CIiIuRTtSkSUeKqVXRTS4urH9bVJdq8mWJKsS6MQCSgtXfXks5yHcm0FCt8VlCu8OttPtHRRN26FTRFDhxI9PGjrE6IYUr288/ce8/CguuvwnBNtU5OBUvTCIWKjoipah4+JKpbt+C3Yc6ckqe5EIlFtPvRbjJaZUTwAKl5qtHc/+aWuFweERGJxSTes4f6qqsTANqopsatlydjSjNq8uDBg1Infi1atEDz5s3L/XixWIzRo0fjyJEjmDBhAnbv3i1VPMHBwXB0dISJiQmio6PBL2ENFbFYDFdXV9y9excXL15E3/yViYvh4eEBT0/PIttlmlGnpCC0e3d0fPQI2QCC27aFzdGjgJ1diQ/1j/LHxIsT8ST2CQDArbYbdn2zC/Wq1fvq4y5eBMaOBRITAR0dYNMmYPx4pVmJgqkisrOB1q2B58+B3r2592VVXwLpt9+ANWu4ptunTwFra0VHxFRFGRnAL78A+T/PTZsCR49y/35NTHoMZv03CyefnwQA1DaujZ19d6Jn3Z5ffZxIJMLy+fOxd/t2PDE1RbWQEJm3xytNjZiiicViGjduHAGgUaNGkUhGoyVcXFwIAL0qZYe/ffv2EQBaUMIU1RVSIyYWk9DNjdrx+dTI0pKCnn59qSEiovTcdJr17yzie/Il87oceHygxBq+7Gyi6dMLrnRatOCWoGEYRQkKIsqfGm/DBkVHo1jXr3PLFwFEMhwTxTDldv48kbl5wdSVf/xRukGOF15dIOv11pLO/N+f/p7iMkqu9s7KzJTbsi1KM2pSkUQiEY0dO5YA0IgRI0gowzr37777jgDQ48ePS7X/uXPnCADNnDmzTOXI7YUMC6Oohw9LNY3H/7/BR54eWao3eHAwUdOmBUnYrFlEpSiOYeRu2zbuPamhQRQQoOhoFCMhgahGDe55mDhR0dEwTIHYWKJvvin47XBzIyrNConFVRjsD9xfqqku5KHKJ2KfJ2HDhg2TaRImEAjI1taWeDwefSxlJ6eFCxcSANpQxktwRc6sH50WTUNODik0JcWV0JLb0cViou3bC2odqlcnuny5AgJmmFISi4m+/ZZ7f9avz02lUpWIxUT9+nHn7+BAlJmp6IgYprD8KS11dbn3qbEx0b59pasd84v0o+Y7mkt+u9y83OhVovynq/h/VToRE4lENGbMGAJAQ4YM4dag+oqEhAQKCQmhhISEQtvv379fJJMWCAQ0a9YsAkC9evUqdN/z588puZjlf+7cuUPa2tqkpaVFYSXNXvd/FJGIicQi2vVoV6FOkPOvzqfMvJK/raOiCn7gAKKePbmrG4ZRNomJ3GLEANG4cYqOpmLl1whqapZuuRmGUZRXr7j1Tj+b3pJK0xCVJ8wrMqhs+e3lJQ4qkyWl6ayvCPkd3vX19TFz5sxip2EYMGAAWrRoUWh/d3d3eHh4SPapXbs2eDwe2rdvj5o1ayIlJQU+Pj549eoVbGxs4OPjA1tb20Llrl27Fl27dkXt2rWhpaWF4OBgXL16FXw+Hzt37sT48ePLdC4VvcTRi4QXmHhhIu5FcGsNOVk5YU+/PWhh2eKrjxMIuCWJ3N2B9HRAQ4PrADxzJusMzSgvb29uqDwRN4fWsGGKjkj+goO5AQs5OcDGjdxnlGGUmUDADfDy8OCmWOHzgalTuWkvjL68ch6A8k2zJCtVurP+6NGjv7pYNwA6cOCAZH93d/diZ9BfvXo1de7cmaysrEhTU5N0dXWpWbNmtGjRIkpKSipSrre3Nw0dOpTq1q1LBgYGpKGhQbVq1aLhw4eTr69vuc6lomrEsgXZtOTmEski3Xor9Gjjg40kFJXcpOvjU7gvWJs2pbtiYRhl8PvvBctrVfaVHbKyuEWXAaLevbnmH4ZRFRERREOHFvzWWFgQHT5c8vtYLBbTkadHyHytOcEDxPPg0ZSLUyglO0Wu8VbpGrHKpCJqxG5/uI1JFyfh1cdXAIC+9fpie9/tsDGy+erj4uKA+fO5GbkBoFo1YPVqYNw4VgvGqA6hkFsP+MEDoF07boWHMiy8oVKmTQO2bQMsLIBnz4Dq1RUdEcOU3bVr3Hv59Wvu706duPd1kyZff9zHrI+Ye20uvJ54AQCsDKywtfdWfNfwO7nEKc/fb6l+Yv/9919ZxcFIKTk7GRPOT0Dng53x6uMrWOhZ4OTgk7gw4sJXkzCRiGuGbNCAS8J4PGDiRODVK25uMJaEMapEXb1gSZ8HD4BipvGrFC5c4H6sAG72cpaEMaqqe3fuQmLlSm5eSh8foEULYO5crmvMl1TTrYYD3x7AjR9voK5pXUSnR2PgyYH47sR3iEyLrLD4ZUGqn9m+ffuiXbt2uHr16hf3yc7OlqYIpgREhBPBJ9BwW0PsfbwXADCx1USETA3BkMZDwPvKLKsPH3L9S6ZP59YIc3Tktu3axdWIMYwqql27YCLJFSu4vmOVSXQ0N6EyAMyZA/T8+lyXDKP0tLSABQuAkBBgwACuguDPPwEHB+DECa7x8ku61OmCZ5OfYZHLIqjz1XHu5TlEpEZUWOyyIFUidu3aNejo6KB3797o0KEDrl+/XmSflStXwsTERJpimK9IyUnBlMtTEJcZh4ZmDXFn7B3s6rcLJjpffs4TE7narnbtgMePAWNjYPt2wNcXcHauuNgZRl6GDeOa1YmAUaOAjx8VHZFsiMXAjz9y59OyJVeLwDCVha0tt47spUuAvT130TF8OLeW7KtXX36cjoYOlndZjseTHmNDzw0V0nlflqTuI/bkyROsXbsWx48fl4xK7Nu3L9TV1REfH4+9e/dCQ0MDcXFxsopZZcmrjfnIsyN4l/wOv3b4FVrqWl/cTywG9uzhrjySk7ltY8ZwIyJZ0wZT2WRmcrW8r14B337LfcGr+jJca9cCv/4K6OoCgYFclwKGqYxycrjfplWrgNxcbvT+3LncIvZ6ehUfjzz7iEmViO3duxeTJ0+GWCwueuBP33g6OjrYtm0bRo8eXf4oK4mKnr7ic48eAVOmAP7+3N/NmnG1YB06VGgYDFOhHj8G2rYF8vK4PlVTpig6ovJ79IirxRYKuQuqMs6ewzAq6e1bYMYM4PJl7m8bG26qlgEDKvbCSmk7669Zswbm5ua4du0akpOTkZmZiYyMDJw4cQK1a9cGEWHBggUsCVOgpCTg55+5Jkd/f64T86ZNQEAAS8KYyq9lS+6qGuD6UwUFKTae8kpPB0aM4JKwwYOBn35SdEQMUzHs7YGLF4GzZ7mmy/BwYOBAoG9fLkmrFKSZ+0JbW5t++eWXYu/LycmhyZMnE5/Pp61bt0pTTKVRkTPri0RE+/cTmZkVzNMyahRRdLTci2YYpSIWc/NsAUSNG6vmEkBjxnDxW1sTFTPtIcNUCZmZRAsXcuvKAkRaWkRLlnBz6smbPH+/paoRs7W1/WLfLy0tLezYsQOurq5Yu3atNMUwZfTkCeDiwnVWTkwEGjfmRo4dPgzUqKHo6BimYvF4gJcXN9/W8+dczZgqOXaMi5/PB44eBdjYJ6aq0tXlRkIHB3PTXuTmAkuXcnOOXbqk6OjKT6pEbPjw4Th58iQufeUZaNasGeuoX0FSU7klThwdgfv3AX194I8/uH4yrq6Kjo5hFKd69YLJinft4vqKCQSKjak0tm/nRkkCXCdlFxfFxsMwyqB+feC//4CTJ4GaNYF374BvvuH6jX34oOjoyk6qRGzevHmoU6cO+vfvj5EjR8LX17fQ/REREThz5gyqsUmp5IoIOHKEG0G1eTM3OnLoUODlS+CXX7jRJgxT1fXoAaxfz9WQ7djBzb+lrNNaCIXcbONTp3JzKo0aBSxZouioGEZ58HjAkCHc79y8edxkzufOAY0aAV+Z2lQpST19RVRUFAYOHAh/f3/weDyYmJigQYMGUFNTQ2BgILKzszF16lRs3rxZVjGrLHmMukhJ4a4Cbt/m/m7QgJspv1s3mRyeYSqdc+eAkSO56S3s7blZ6hs2VHRUBZKTuQup69e5H5uVK7kpK1R96g2Gkafnz7mLl5cvuVtJC4iXldJOX5GPiHDu3DkcO3YMd+7ckTRFamtrY+TIkdi8eTN0dHSkDlbVyeOFJAK6duVmxF+8mOv/ovXlqcQYhgE3erJ/f64Zw9AQOH4c6N1b0VFxc5716weEhnJzJR09ys2BxjBMyYi4UZW2trI/ttInYv8vLS0NWVlZMDc3h5qamqwPr7Lk9UK+fcs1P9p8fV1vhmE+k5AADBoE3LnDdYRftw6YPVtxNU/XrnE1YSkp3Gf5/HmgeXPFxMIwTGFKO4/YlxgaGsLS0pIlYRXE3p4lYQxTVubmXPPfTz9x/Sp/+YX7f25uxceybRtXI5eSArRvD/j5sSSMYaoKuSRiDMMwqkBTk5ulfuNGrlbswAGuqT8+vmLKFwi4EZzTpnGd8n/8Ebh5k5tqg2GYqoElYgzDVGk8Hjfty+XLXAffe/eA1q2BZ8/kW25SEtCrFzeCk8fjVgDw8mJ9PBmmqmGJGMMwDLjpLB4+BOrW5Tr8tm/PLasiDy9fAm3acLVf+vpcOfPns5GRDFMVsUSMYRjmEwcHwNeXa57MzAS++46bPkKWQ5quXuUWIn/zBqhdm5t8uX9/2R2fYRjVwhIxhmGYz5iaAv/+y02mCnAz2o8aBWRnS3dcIm7C5d69uVUwOnTgkr6mTaWPmWEY1SVVIhYeHo60tLSv7pOeno7w8HBpimEYhqlQGhrcxMjbtwNqatx6j507AzEx5TueQABMnsz1RROLgTFjgBs3uKWXGIap2qRKxOrUqYNNmzZ9dZ/t27ejTp060hTDMAyjED//zDUlmphwU0q0bg0EBJTtGB8/cssr7d7N9QH74w9g/37WKZ9hGI5UiRgRoaT5YOUwXyzDMEyF6dKFS8IaNgSioriFt0+eLN1jQ0K4Tvne3oCBAbec0i+/sE75DMMUUJd3AZGRkTAwMJB3MQzDMHJTty7w4AEwYgTXf2zYMG5tO3d3bv6x4ly5wu2XlgbUqcMlYY0bV2zcjGIQEUhMEAvFIBFBLBID+XUSPIDH4xX6Fyi67Yv3MZVOmROxpUuXFvrb29u72P1EIhEiIyNx/PhxtGnTplzBMUxZifJEiHsWh+zkbMkXHxEBVMZ/gWLv+6piviO/+MVZzGbL5paoVr9aqc+VqVhGRlwy9euvwJ9/AkuXcsnYwYPcupD5iIBNm7iaL7GYq0E7fZqbyZ9RTsJcIQJ2BSDknxAIc4SS5OnzRKo0/yfRp+RLLOeWoM8SNDUtNdRoWQPWHa1h09EG1u2toVtNV77lMzJV5rUm+Z9d/vF4vBKbHq2srHDmzBm0bt26fBFWIvJcq6qqyk7ORuSDSITfDUfEvQhE+UVBmCNUdFjlwwNajG6Bzp6dYWRjpOhomK84cACYNInrhN+iBbcupLU1kJfHjbbcu5fb76efuA7/mpoKDZf5AhITgo8H4+bvN5HyPkXR4ciMWUMz2HS0kdyM6xiz2jQpKdWi37dv3wbA1RB06dIFY8aMwejRo4vsp6amBlNTUzg4OBRK3qoylohJh4iQ8j4F4XfDEX6PS7wSnicU2U/HVAeGtT49v8VU9Rep5i/jPl8OsOT4v0SYLUSUXxQAQE1LDc7TneGywAU6pjpfPyijMHfvAgMHcouHW1gA+/ZxC4ffvs01V/7xBzBrFusPpqzeXn2L679eR+yTWACAfg19dPytI4xrG4OnxgNfnQ++Gr/0/1fjg69e/P95ajzw+Lwv1rJ/qXa+NPflpuUi8mEkIu5FIPxuOBJDEoucq34Nfdh0sJHUmlk2twRfnf0ul4VSJWKf8/T0hJubGzp16iTLmCotloiVjUggQuzjWEnSFXEvAhmxGUX2M61nylXJd7CGTQcbVGtQTSWv/qL8onD91+v44P0BAKBlpIWOCzqizYw20NDRUGxwTLE+fOAmYw0KKthmYAAcPw706aOwsJiviA6Ixo3fbuDd9XcAAC1DLXT4rQPazmwLDV3V/5xlJWYh4j6XlIXfDUf0o2iIBeJC+2joacC6nTWXmHWwQa22taCpz6ptv0ZpEzGmbFgi9nU5KTmIeBAhubKL8ouCMLtwMyNfgw8rJytJ0mXd3hp61fW+cETVQ0R4+x93pR73LA4AYFDTAJ09O6PF6BbsKlYJZWRwE76eOwfY2XH9yBo1UnRUzP9LepuEW7/fQvDxYACAmqYaWk9tDZdFLpW6T5UgW4DoR9Fc9427EQi/F47c1NxC+/DUeLBsYSlpyrTuYA2DGmyQ3eeUOhHLy8vD2bNn4e/vj5SUFIhEoqKF8HjYt2+fNMVUCiwRK0BESPmQIkm6Iu5FIP55fJHmPR1THVi3t+YSr442sHKygrq23Af7KhyJCUHHgnDz95tIDUsFwPX76LqqKxr0b6CSNX6VmVjMLVXUrBlQxT/aSiczPhM+y33waOcjrmaIBzQb2Qxuy9xgXNtY0eFVOBITEl4kSGrMwu+GS75jPmdiZ4K6veui/dz2VfJ5+n9Km4iFhYWhe/fuePv27Vf7v/B4vGITtKqGJWJcAvb+xnvcWnILkQ8ii9xvWtdUknRZd7CGWQMzrm9FFSXMFcJ/uz/uLL+D7CRujR3r9tbotqYbbDraKDg6hlFeeRl5eLD+Ae6vu4+8jDwAQN1eddF1VVdYtrBUcHTKJTUilbsovsfVmsU+jZVcFPPV+Wg+ujlcFrrAxM5EsYEqkNImYgMHDsTZs2fxww8/YNy4cahVqxbU1YuvrbC1tS13kJVFVU/Ewu6E4dbiWwi7HQaA+4DnNzNad7CGdXtr6FvoKzhK5ZSTmoN7a+/h4YaHkubaBv0boMvKLqjemK2TwzD5RAIRAvcE4rbnbWTGZwIAajjWQPe13VGnC1vlpTRyUnMQfjccvht9JX3peGo8NP+hOTou7Ihq9areNDtKm4gZGxujdevWuHbtmixjqrSqaiIW6RuJW4tv4d017gOtpqkGx8mOcFngAn1LlniVRXp0Orw9vfF432OQiMDj89B8THN09ugMI2s25QVTdRERXvz9AjcX3kTSmyQAgIm9Cbqu7IpGgxtV6Zp1aUTcj4DPMh+8ufIGAMDj89D0+6ZwWeQCMwczBUdXcZQ2ETM0NMTkyZOxdu1aWcZUaVW1RCzmcQy8l3jj9cXXALgasJbjW6LTok4F00sw5ZL4MhE3F91EyD8hAAB1bXU4z3BGx986QseETXnBVC3vb73H9V+vI9o/GgCgV10PnZZ0guMER6hpqik4usohyi8KPst8JN/n4AFNhjdBp987wbxR5Z+tWGkTsZ49e0JTUxMXLlyQZUyVVlVJxOKD4+Ht7i1JEnhqPDT/sTk6Le4EkzpVt4+BPEQ+jMS1+dcQficcAKBtrI2OCzvCeZozm/KCqfRin8bixm83JLU1GnoaaD+vPdrNaQctA7aqujxEB0TDZ5kPXp17xW3gAY0GN0KnxZ1g0dRCscHJkdImYo8fP4aLiwu8vLwwePBgWcZVKVX2ROzj64/w9vDmhocTAB7QdERTuLq7sqV75IiIEHo5FDd+u4H44HgAgGEtQ3T27Izmo5uDr8amvGAql5QPKbi1+BaeHX0GEFfb7jjJEZ0Wd2L9TCtI7JNY+Cz3QcjpEMm2hgMbotPiTpVyMITSJmJLly6Fv78/Ll++DFdXV7Rs2RJGRkX7qfB4PCxevFiqQCuDypqIJb9Phs9SHzw99FSyxlqjwY3g6uHKOpJXILFIjGdHnuHW4ltIi0gDAJg3MkfXVV1Rv199NuUFo/LEIjFuLLwB342+EOVxI/EbD2uMLsu7wLSuqYKjq5riguJwZ8UdPD/5XDLSskH/Bui0uBOsnKwUG5wMKW0iVtqli9j0FZzKloilRqTCZ7kPnux/ArGQm7m5fr/6cFvqVimviFSFMEcIv21+uLPiDnKScwAAbsvc0Ol3tgIGo9puud+Cz1IfAECdLnXQbU23SvVjr8oSXiTgzoo7CD4eLLkgr9enHjot6YRabWopODrpKW0ilr/uZGm4urqWt5hKo7IkYukx6bi76i4CdgVIrkrte9jDbZkbajrXVHB0TL6clBzcXnYbD9c/BHjAjzd+RB03NnyfUU3vrr/D4R6HAQL67emHlj+1ZLW8SijxVSLurryLZ0eeSRIy+572cF3iCuv21gqOrvyUNhFjykbVE7HMhEzcW3sP/tv8JXNZ2braosvyLmxyUSV27qdzeLL/CfQt9THpySTWh4ZRORmxGdjZYicy4zLRcnxL9N/TX9EhMSVIepOEOyvvcF1WRFyaUadrHbi6u8LWRfXmFWWJWCWhqolYdnI27v9xH76bfCHIFAAAarWrBbdlbqjTpQ67KlVygiwB9jjvQcLzBNh1t8OoK6PYnEqMyhCLxDjS4wje33yP6k2qY7zv+EqxOHdVkfwuGXdW3cFTr6eSLiy1O9dG+/ntUbdnXZX5LpLn77fUw6mEQiE2bNgAZ2dnGBoaFppZ/8mTJ5gyZQpev34tbTGMAmQmZMLb0xub6mzC3ZV3IcgUoIZjDXx/+XuMuzcOdl3tWBKmAjR0NTDk5BBo6Grg3bV3uLPqjqJDYphS81nug/c330NDTwNDTg1hSZiKMbEzQf89/TH9zXQ4TnYEX4OPD94fcKzPMWx12Arfzb7ITcst+UCVmFQ1YtnZ2ejRowfu378PMzMzaGhoICYmRtIxPzU1FZaWlvjll1+wfPlymQWtqlSlRiwmMAZ+W/wQ9FcQRLnca1m9aXW4LXVDg2/ZgtOq6onXE5wbew48Pg+jb42GbSfVax5gqpb3N9/jULdDAAEDDg1A8x+aKzokRkqpEal4uPEhHu97jNxULgHT1NdEi7Et4DzNWWmnOlLaGrGVK1fi3r17WLVqFWJjYzF+/PhC9xsZGcHV1RX//fefVEEy8icSiBB8Ihj7O+7HbsfdeOL1BKJcEaycrDDo+CBMfjIZDgMcWBKmwlqMaYHmPzYHiQmnR5xGZkKmokNimC/KiMvAPyP/AQhoMbYFS8IqCSNrI/T8syfmRM5B3x19YdbQDHkZefDb4oetDbbiaO+jCP03VNLRvyoofoXuUjpx4gQ6d+6M+fPnA0CxP9J2dnZ4/PixNMUwcpQZn4mA3QF4tOMR0qPTAXCTIzYa0ghtZrRBzTY1WfJVifTZ1gdRflFIfJmIsz+exfeXvleZPhpM1SEWiXFm1BlkxGbAvLE5+mzto+iQGBnT1NeE02QnOE5yxPsb7+G3xQ+vLrzCmytv8ObKG5jWM4XzNGe0GNMCWoaVe5UEqRKx8PBwfPfdd1/dx9DQEKmpqdIUw8hB9KNo+G3xQ/DxYMkUFHoWepIPhkENAwVHyMiDpr4mBp8cjL3Oe/HmyhvcW3sPHX/rqOiwGKaQu6vu4t31d4X6NzKVE4/Hg103O9h1s0Pyu2T4bfPD432PkRSahCszr+DmoptoPqY5nKc5w6xB5VxkXKpEzMDAAAkJCV/d5+3btzA3r/wLgqoCUZ4IL06/gN8WP0Q+iJRsr+lcE84znNFocCOoa0n1lmBUgEVTC/Te2hsXxl/Azd9vwqajDZt+hFEaH25/gLe7NwCuBrcqLCjNcEzsTNDzz55w83TDsyPP4LvZF4khifDf6g//rf6o26sunKc7o24v1RltWRpS/eq2bdsWFy5cQGpqarFLG0VGRuLy5csYMGCANMUwUsqIy0DArgA82vkIGTEZAAC+Bh+NhzaG83TnSjHrMVM2Lce1xIdbHxB0NAinR5zGpCeToFtNV9FhMVVcZnwmTo84DRITmo9ujhZjWig6JEYBCjVb3nwPv83/12xZ1xTO0ytPs6VUoyZ9fHzg5uaGVq1aYdOmTfj333+xcuVKpKen48GDB5g+fTrevHmDBw8ewNHRUZZxq6SKHjUZ5RfFNT+eCIZYwM3fom+pD8fJjnCa5AR9SzaxZ1WWm56LPU578PH1R9TrWw8jzo+oVFeZjGohMeFo76N4e/UtzBqaYYL/BGjqaSo6LEZJfN5s+floy4pqtlTqCV137tyJGTNmFLuWpJqaGrZv315kNGVVVRGJmChPhOennsNvix+ifKMk22u1rQXn6Vzzo5qmmlzKZlRP7NNY7G2zF6JcEbqv6472c9srOiSmirqz8g5uLroJdR11TPCbgOpNqis6JEYJ5WXkFWq2zGff0x5tZrSRW7OlUidiABASEoKdO3fC19cXSUlJMDQ0RJs2bTBlyhQ0btxYFnFWCvJ8IdNj0hGwKwABuwKQEcs1P6ppqqHxMK75sWZrtgYkU7xHux7h0uRL4KvzMfbOWNRqy5qqmYoVdicMBzsfBIkJ/ff1R8txLRUdEqPkiKhQsyU+ZTKmdU3x3eHvZP49pvSJGFM68ngh8zLzcHHiRTw/9byg+bGGPpx+doLjREe2riBTIiJuXrHnJ57DyMYIkx5Pgo6pjqLDYqqIzIRM7Gq5C+lR6Wg2qhkGHBrApsxhyiT5XTL8t/sjcG8gBFkCzAqbJfOR//JMxNgQORWnoauBhBcJEAvEsG5vDecZzmg4sCHUNFjzI1M6PB4P/Xb3Q0xADJLeJOHc2HMYdnYY+zFk5I7EhLM/nkV6VDqqNaiGvjv6svcdU2Ymdibo8UcPdPbojMiHkSo3/RKrEatA8sqow3zCoKGnAStHK5kdk6l6Yh7HYF/bfRDlidBjfQ+0m91O0SExldzdNXdx47cbUNdWx3jf8bBoZqHokBimWEpTIzZu3DjweDysXLkSFhYWGDduXKkex+PxsG/fvnIFyJSMrRnIyEKNljXQc0NPXJ56Gdd/vQ6bDjao6cz6FjLyEX4vHDcX3QQA9NrciyVhTJVVphoxPp8PHo+HkJAQ1K9fH3x+6Zaq5PF4xY6qrGpUZdFvpuoiIvw99G+8+PsFjGsbY9LjSdA21lZ0WEwlk/UxC7ta7EJaZBqajGiCgUcHsiZJRqkpTY3Y+/fvAQA1a9Ys9DfDMJUDj8dDv739EBMYg+R3yTg37hyGnh7KfiQZmSEx4ezos0iLTINpPVN8s+sb9v5iqrQyJWK2trZf/ZthGNWnbaSNwScGY1/7fXh55iX8tvqhzfQ2ig6LqSQerH+A0EuhUNNSw5BTQ6BloPozozOMNErXtvgF9+7dw5w5cxAbG1vs/bGxsZgzZw4ePnwoTTEMw1QwKycr9PijBwDg2txriA6IVnBETGUQ8SAC13+7DgDotbEXLJtbKjgihlE8qRKx9evX48KFC7C0LP7DZGlpiYsXL2LDhg3lLiMqKgobN25Ejx49YGNjA01NTVhaWmLQoEHw9fUt9XG8vb3B4/G+ePtSsujv748+ffrAxMQEenp6cHZ2xrFjx8p9PgyjKpynO8PhOweI8kT4e+jfyEnNUXRIjArLTsrG38P+BokIjYc1huMktuwdwwBSziPm7++Prl27fnWfTp064dq1a+UuY8uWLVizZg3s7e3RvXt3VK9eHaGhoTh79izOnj2Lv/76C0OHDi318VxdXdG5c+ci22vVKjoLr7e3N3r27AlNTU0MHz4cRkZG+OeffzBy5Eh8+PABCxcuLPd5MYyy4/F46L+vP2IfxyL5XTIuTLiAwScGs/48TJkREc6OOYu0iDSY1jVFv9392PuIYT6RKhGLj4+XdNz/EktLS8THx5e7DGdnZ/j4+MDFxaXQ9jt37qBr1674+eef8e2330JLq3T9DDp37gwPD48S9xMKhRg/fjx4PB58fHzQsiW35Ia7uzvatWsHd3d3DBkyBPXq1SvzOTGMqtAx0cHgE4Oxv+N+vDj1Ao86P0LrKa0VHRajYh5ueIjXF15DTVMNg08OhpYh6xfGMPmkapo0NjZGeHj4V/cJCwuDvn75l9kZOHBgkSQMAFxcXODm5oakpCQEBQWV+/hfcvPmTbx9+xbff/+9JAkDAAMDAyxevBhCoRAHDhyQebkMo2xqOtdEtzXdAAD/zf4PMY9jFBwRo0oifSNx/VeuX1jPDT1Ro2UNBUfEMMpFqkSsXbt2OHPmDCIiIoq9Pzw8HGfPnkX79u2lKeaLNDQ0AADq6qWv2AsNDcXmzZuxevVq/PXXX0hMTCx2P29vbwBAjx49ityXv+327dtljJhhVFPbWW3RoH8DSX+x3LRcRYfEqIDsZK5fmFgoRqPBjeD0s5OiQ2IYpSNVIjZnzhxkZWWhQ4cOOHToEGJiuCvlmJgYHDx4EB06dEB2djZ++eUXmQT7ufDwcFy/fh2WlpZo2rRpqR937NgxzJw5EwsWLMD3338PGxsbrFu3rsh+oaGhAFBs06OJiQnMzMwk+3xJbm4u0tLSCt0YRhXxeDx8e+BbGNkYIelNEi5Ougi2OhrzNUSEc2PPITUsFSZ2Jui3l/ULY5jiSJWIubi4YPPmzYiJicHYsWNRq1YtqKuro1atWhg3bhxiY2OxadMmdOrUSVbxAgAEAgF++OEH5ObmYu3atVBTK3mBa3Nzc6xbtw4hISHIzMxEVFQUjhw5AlNTU8yfPx+7du0qtH9qaioAwMjIqNjjGRoaSvb5klWrVsHIyEhys7a2LuUZMozy0THVwaDjg8BX5yP4eDAC9wQqOiRGiflu9sWrc68k/cK0jdgKDQxTHJks+h0cHIwdO3bA398fKSkpMDY2hrOzMyZPnowmTZrIIk4JsViM0aNH48iRI5gwYQJ2794t1fGCg4Ph6OgIExMTREdHS5Zt6tGjB65du4bQ0FDUrVu3yOPs7e0RGRmJ3NwvN9Hk5uYWuj8tLQ3W1tZsiSNGpd1bdw/X519nCzUzXxTlH4X9HfZDLBCj1+ZebEJgRuUpzRJHX9KkSRNs27ZNFof6KiLChAkTcOTIEYwaNQo7d+6U+phNmjRBmzZtcOfOHbx58wb169cHUFAT9qVar/wX5Wu0tLRKPZqTYVRF+1/aI+x2GEIvheLvYX9j8rPJUNMouVaaqTouT70MsUCMhgMbwnmas6LDYRilJlXTZEUSi8X46aefsH//fowYMQJeXl6lXnS8JGZmZgCArKwsybb8vmHF9QNLTk5GYmIim7qCqZJ4fB4GHBwAXXNdJL5MxKvzrxQdEqNEovyjEO0fDTVNNfTd0Zf1C2OYEqhEIiYWizF+/HgcOHAAw4YNw+HDh0vVL6w0hEIhAgMDwePxYGNjI9nu6uoKALh69WqRx+Rvy9+HYaoa3Wq6cJzIzYzuv81fwdEwyuTR9kcAgMZDG0Ovup6Co2EY5Vempslx48aBx+Nh5cqVsLCwwLhx40r1OB6Ph3379pUrwPyaMC8vLwwZMgRHjhz5ahKWmJiIxMREmJmZSWq6AODBgwdo27ZtoaszoVCIefPmISwsDL169YKpqankvq5du8LOzg7Hjh3DjBkz0KJFCwBAeno6li1bBnV1dYwZM6Zc58QwlYHjJEfcXXUXH259QEJIAswbmis6JEbBsj5mIfh4MACg9VQ28S/DlEaZOuvz+XzweDyEhISgfv36pW4a5PF4EIlE5QrQw8MDnp6e0NfXx8yZM4udM2zAgAGSRCl/f3d390Iz6NeuXRs8Hg/t27dHzZo1kZKSAh8fH7x69Qo2Njbw8fGBra1toePeunULPXv2hJaWFkaMGAFDQ0P8888/eP/+PZYvX45FixaV6Vzk2dmPYRThxHcn8PLsS7Se1hp9tvRRdDiMgt3/4z6uzbsGy5aWmBgwkTVLMpWG0nTWf//+PQBIljXK/1uePnz4AADIyMjAihUrit2ndu3akkTsS37++WdcuXIF3t7eSExMhLq6OurWrYtFixbhl19+gYmJSZHHuLm54e7du3B3d8fJkyeRl5eHxo0bY9myZRg5cqS0p8YwKq/11NZ4efYlnh58iq4ru0LLgA1OqapITHi0g2uWbD21NUvCGKaUylQj5uPjg9q1axfqS8WUHqsRYyobIsK2htvw8dVH9NneB61/Zs1RVVXov6E41ucYtIy08Ev0L9DQ1VB0SAwjM/L8/S5TZ303Nzd4eXlJ/u7SpQsOHTok04AYhlEdPB5Psgi4/zZ/Ntt+FZY/aKPF2BYsCWOYMihTIqaurg6hUCj529vbW9J0yDBM1dR8dHNo6Gog4XkCwnzCFB0OowDJ75MRepmb6ofVijJM2ZQpEbO2tsa9e/cgFosl21g/AIap2rSNtNF0FLfea/7UBUzVErArACDArrsdqtWvpuhwGEallKmz/vDhw7Fy5UqYmJigWjXuw7ZhwwYcOHDgq4/j8Xh4+/Zt+aNkGEapOU91RuDuQIT8E4L0mHQY1DBQdEhMBRHmCBG4l1t3lE1ZwTBlV6YaMXd3dyxfvhzNmjUDj8cDj8cDEZV4+7wGjWGYyseimQVsOtpALBQjYHeAosNhKtDzU8+R/TEbhtaGqN+3vqLDYRiVU6YaMQ0NDSxcuBALFy4EwM0rNnv2bCxZskQuwTEMozpaT22N8LvhCNgVAJeFLmz9ySoiv5O+02Qn8NVVYrEWhlEqZfrU+Pj4IDw8XPK3u7s7OnfuLOuYGIZRQQ0HNoSehR4yYjLw8uxLRYfDVIDogGhE+UaBr8FHy59aKjochlFJUk1fcfv2bTZqkmEYAICaphpaTWgFgK0/WVX4b+de50aDG0HfQl/B0TCMamLTVzAMIzNOk5zAU+Mh7HYY4p/HKzocRo6yk7MRfIytK8kw0pJ6+gqGYZh8hrUM4fCtA4CC2hKmcnri9QTCHCEsmlnAur21osNhGJVVpkRs+PDhuHXrFkxMTGBnZwcA2LhxI+zs7L56s7e3l0vwDMMon/zakWeHniE3LVfB0TDyQGKSzBnH1pVkGOmw6SsYhpGp2m61YeZghryMPDw9/FTR4TBy8O76OyS9SYKWoRaaft9U0eEwjEorUyKWP33FnTt38PbtWxARZs+ejffv35d4YximauDxeHCa4gSArT9ZWeUPxmg+pjk09TUVHA3DqDapJn1h01cwDFOc5j82h4aeBhJDEvHB+4Oiw2FkKCUsBa8vvgbA1pVkGFmQOhFr3749NmzYAGdnZxgaGkJdvWCO2CdPnmDKlCl4/fq11IEyDKM6tI200eyHZgDY+pOVTcCuAJCYUKdLHZg5mCk6HIZReVIlYtnZ2XBzc8PcuXMRFhYGQ0PDQs0QderUwYEDB3Do0CGpA2UYRrU4T3UGAIScCUFaVJqCo2FkQZjL1pVkGFmTKhFbuXIl7t27h1WrViE2Nhbjx48vdL+RkRFcXV3x33//SRUkwzCqp3qT6rDtZAsSEVt/spIIOR2CrIQsGNQ0QIP+DRQdDsNUClIlYidOnEDnzp0xf/58ySjK/2dnZ1doWSSGYaqO/E77gbsDIcoTKTgaRlr5nfQdJzmydSUZRkak+iSFh4ejdeuvV08bGhoiNTVVmmIYhlFRDb9rCH1LfWTEZiDkTIiiw2GkEPskFhH3I8BX56PV+FaKDodhKg2pEjEDAwMkJCR8dZ+3b9/C3NxcmmIYhlFRappqaDWR+9FmnfZVW/5KCQ0HNoRBDQMFR8MwlYd6ybt8Wdu2bXHhwgWkpqbCyMioyP2RkZG4fPkyBgwYIE0xVZpAIIBIxJp0GMVTU1ODhoZGmR/nONERd1bcQZhPGOKC4mDR1EIO0THylJOSg6CjQQBYJ32GkTWpErF58+bBzc0N3bp1w6ZNmyQLgmdlZeHBgweYPn06BAIB5syZI5Ngq5K0tDQkJiYiN5ctEcMoDy0tLZiZmcHQ0LDUjzGsaYiG3zXEi79fwH+7P77Z8Y0cI2Tk4cnBJxBkCWDe2Bw2LjaKDodhKhWpErFOnTph27ZtmDFjBlxcXCTbDQy4ams1NTVs374djo6O0kVZxaSlpSEqKgr6+vowMzODhoYGW8uNUSgigkAgQGpqKqKiogCgTMmY0xQnvPj7BZ4dfoZuq7tB20hbXqEyMkbE1pVkGHmSKhEDgMmTJ8PV1RU7d+6Er68vkpKSYGhoiDZt2mDKlClo3LixLOKsUhITE6Gvr49atWqxLz1Gaejo6MDAwACRkZFITEwsUyJWu3NtmDcyR8KLBDw99BRtpreRY6SMLL2/+R4fX3+EpoEmmo1qpuhwGKbSkToRA4CGDRti06ZNsjhUlScQCJCbmwszMzOWhDFKh8fjwcjICFFRURAIBKXuM5a//uS/0/7Fo+2P4DzNmb2/VYRkXckfm0PLQEvB0TBM5cMmglEy+R3zy9MpmmEqQv57s6yDSJr/wC0QnfgyER9ufZBDZIyspUWm4dW5VwAAp5+dFBwNw1ROMknE7t+/j4kTJ8LZ2RkNGjRA69atMWHCBNy9e1cWh6+SWG0Bo6zK+97UMtRCsx+5pq38WhZGuT3a9QgkJti62qJ64+qKDodhKiWpE7G5c+fCxcUFe/fuxaNHj/D27VsEBARg3759cHV1ZSMmGYaRaD2Fm/rg5bmXSItk608qM1GeCIF72LqSDCNvUiVihw4dwvr169GgQQP89ddfiImJgVAoRGxsLI4fPw4HBwds2rSJLfrNMAwAoHrj6rB15daffLSLTfCqzEL+CUFmXCb0a+jDYYCDosNhmEpLqkRsx44dsLa2hq+vL4YNGwYLC26ixurVq2Po0KF48OABatWqhe3bt8skWIZhVF9+7UrgHrb+pDLLn0nfcaIj1DTUFBwNw1ReUiViwcHBGDRokGTesP9naGiIgQMH4vnz59IUwzClJhaL0bx5c/Tp00fRoSiVN2/eQF1dXSkuihwGOMDAygCZcZl4cfqFosNhihEXFIfwO+HgqfHQagJbV5Jh5EnqPmJE9NX7WadzpiJ5eXnh2bNn8PDwKPZ+f39/9OnTByYmJtDT04OzszOOHTtWpjJq164NHo9X7G3y5MlljrlZs2bg8XgIDw8v82NLq27duhg5ciQ8PDyQlqbYvllqGmz9SWUnWVfyu4YwrFn6+eIYhik7qeYRa9KkCU6fPo1ly5ZBX1+/yP3p6ek4ffo0m9SVqRAikQienp5wdXWFs7Nzkfu9vb3Rs2dPaGpqYvjw4TAyMsI///yDkSNH4sOHD1i4cGGpyzIyMsKsWbOKbHdyKtsQ/5ycHISEhMDMzAw2NvJdOmbevHk4dOgQNm/ejN9//12uZZXEcYIj7iy/g/C74Yh7FgeLZmz9SWWRm5aLZ4efAeBWRGAYRs5ICl5eXsTj8ahJkyb0999/U0JCAhERJSQk0KlTp6hJkybE5/PJy8tLmmIqjdTUVAJAqampX9wnOzubXrx4QdnZ2RUYWeVw/vx5AkB79+4tcp9AICB7e3vS0tKiwMBAyfa0tDRq3Lgxqaur0+vXr0tVjq2tLdna2sok5ocPHxIA6tGjh0yOV5LmzZuTjY0NiUSich9DVu/Rk0NOkgc86PzE81Idh5Et3y2+5AEP2tpwK4nFYkWHwzBKoTS/3+UlVdPk6NGjMXPmTDx//hxDhw6FhYUFNDQ0YGFhgWHDhuH58+eYNm0aRo8eLYuckamCMjIy4OHhgQYNGkBbWxt169bFnj17AAA3btwAj8fDpUuXAHDNkjweD4MGDSpynJs3b+Lt27f4/vvv0bJlS8l2AwMDLF68GEKhEAcOHKiYk/pMYCA3PUBFrcc6dOhQhIeH48aNGxVS3tfkd9oPOhKEnJQcBUfDAFxXk/xmydZT2LqSDFMRpF7iaMOGDRg0aBAOHDiAJ0+eIC0tDYaGhmjZsiVGjx5daDFwhimLyMhIdOvWDW/fvsXQoUPRt29fHD16FJMmTYKTkxNWrFiB1q1bo2/fviAieHt7w8HBAcbGxkWO5e3tDQDo0aNHkfvyt92+fbvUseXm5uLgwYOIioqCiYkJ2rdvj+bNm5f5HAMCAgAArVpVTIfodu3aAeAS0+7du1dImV9i28kW5o3NkfA8AU8OPkHbmW0VGg8DfPD+gMSQRGjoaaDZD2xdSYapCDJZa7Jjx47o2LGjLA7FMAC40Y+DBg3Cq1evcO7cOfTv3x8A0KdPH3Tv3h1r167FrVu3cPnyZQBASEgIkpKS0Lt372KPFxoaCgCoV69ekftMTExgZmYm2ac0YmNjMWbMmELbevXqhcOHD8PMzKzUx6noGrH8Pmz379+vkPK+hsfjofXU1rg85TIebX+ENtPbgMdnNTCKlD94otkPzaBtpK3gaBimapBJIsZUECIgK0vRUZRMVxeQsknj/Pnz8PPzw7BhwyRJGFCQSBw/fhxt2rSRJF6RkZEAIJnL7v+lpqYC4DrZF8fQ0FByjJKMGzcOrq6uaNy4MbS0tPDixQt4enri33//Rf/+/XHv3r1SNenk5eUhODgYJiYmqFOnTqnKlpaBgQG0tbVLfa7y1mxUM1z/9To+vv6I9zffw66bnaJDqrLSo9MRciYEAND6ZzaTPsNUFKkSsXv37uH06dOYP38+LC0ti9wfGxuLtWvXYujQoWjbljU7SC0rCyhmdKrSycgA9PSkOkT+lBIzZswotF1TU1Pyf09PT8n/P378CICr3ZK3JUuWFPq7TZs2uHjxIlxdXXH37l1cvnwZffv2LfE4z549g0AgkHmz5Pjx40FE2LdvX7H3m5qaIjExUaZllpeWgRaa/9gc/tv84b/NnyViChSwOwAkIth0tGGjWBmmAknVWX/9+vW4cOFCsUkYAFhaWuLixYvYsGGDNMUwVZCPjw9MTEyKJPD0ad66du3aoWfPnpLtOjo6AIDs7Oxij5dfE5ZfM/b/0tLSvlhbVhp8Ph9jx44FwF2glIa8miX9/f2/Oo1GdnY2dHV1ZVqmNPLXn3x1/hVSw4t/fRj5EglECNjN9Vdk60oyTMWSKhHz9/cvsW9Yp06d8PDhQ2mKYfLp6nK1Tcp+k/JHPjU1FXFxcWjQoAH4/MJv0fw+Yd98802h7ebm5gCApKSkYo+Z3zesuH5gycnJSExMLLb/WFnk9w3LKmXz8Zc66guFQnh4eMDe3h7a2tqoWbNmoTnOFixYgIYNG0JXVxc1a9aU1Azm5ORATU0Nz549w5QpU8Dj8TBs2LBCxxaLxUhNTZU8X8rAvJE5arvVBonZ+pOK8vLsS2TEZEDPQg8NBzZUdDgMU6VI1TQZHx+PmjVrfnUfS0tLxMfHS1MMk4/Hk7rJTxXkJzJqaoXXt8vJycFvv/0GAFBXL/zWbdy4Mfh8/hc73Lu6umLVqlW4evUqhg8fXui+q1evSvaRhq+vLwBu5v3S+FKN2LJly/Dvv//i4MGDsLa2xrt375CQkAAAEAgE0NLSgpeXFywtLeHr64sxY8bAxcUFnTt3xvXr19G1a1eEhoZCT08Pev/3fgkNDYVYLEbTpk2lOldZaz21NT7c+oDAPYFwXeIKdS3WfbUi5XfSbzWhFdQ02bqSDFORpKoRMzY2LnFZlrCwsGJn3WeYLzE3N4e2tjYCAwMRFhYm2T579my8e/cOAIok98bGxmjWrBkePXpU7LJbXbt2hZ2dHY4dO4YnT55Itqenp2PZsmVQV1cvMgry7du3ePnyJQQCgWTbixcvkJKSUuT4d+/exfr166GlpYWBAweWeI4CgQBBQUEwNDSEvb19ofuuXbuGgQMHomPHjrC1tYWbmxuGDh0KANDQ0ICHhwfatGkDW1tbDB06FC1btsSrV6/A5/MRExOD2rVrw97eHpaWlkXWgc1PFqVNOmXN4Vtu/cmshCyEnA5RdDhVSvzzeHzw/gAenwfHiRUzepdhmAJSJWLt2rXDmTNnEBERUez94eHhOHv2LNq3by9NMUwVo66ujpEjRyI7OxsuLi6YOXMmunTpgp07d2Lp0qXQ1dXFvn378PvvvyMzM1PyuAEDBiA1NRX+/v7FHnPv3r0Qi8VwcXHBxIkTMXfuXDRv3hzPnz+Hh4cH6tevX+gxXbt2RcOGDREVFSXZdvLkSVhZWaFfv36YPn065s6di169eqFTp04QCATYunVrqZYqev78OXJzc9GqVasiIyz79u2L33//Hf3798eRI0eQkZEhue/du3eYNGkSGjVqBGNjY+jr6+Phw4eoUaMGAG4AwNfmM7t27RrU1NSKNO0qGl+dD8dJXBLgv63o68fIz6MdXG1Yg28bwMi6/P0kGYYpJ2mm5ffx8SE+n0/W1tZ08OBBio6OJiKi6Oho8vLyolq1apGamhrdvn1b2hUAKgW2xFHpZWRk0NSpU8nS0pI0NDSoVq1atGnTJiIiOnDgAJmZmZGBgUGhJVgiIyNJTU2Npk+f/sXj+vr6Uq9evcjIyIh0dHTIycmJjhw5Uuy+tra2BIDev38v2ebt7U1Dhw6lunXrkoGBgSS24cOHk6+vb6nPb+/evQSA5syZU+z9z58/p+XLl5OdnR1ZW1tTSkoKxcXFUbVq1ejHH3+k69ev04sXL+jmzZsEgN6+fUtERL169aIlS5YUe8zMzEzS19enAQMGlDrO4sjrPZoWnUZL1ZeSBzwo5nGMTI/NFC8nLYdWGqwkD3jQ22tvFR0OwygteS5xJFUiRkS0detWUldXJz6fT3w+n9TU1CT/V1dXp61bt8oizkqBJWLyN2LECKpWrRplZGQoOhSZiIqKIgAUHBxMe/fupRo1ahS6f/78+WRoaChJSGvWrEmnTp0q9lj79u0jAFJfGMnzPXpq2CnygAedG39O5sdmivLb7kce8KAt9beQWMTWlWSYL1HatSYBYOrUqXj8+DEmT54MR0dH2NnZwdHRET///DMeP36MqVOnSlsEw5TaihUrkJGRgW3btik6lHJZs2YNjh49ilevXuHFixdYvHgxHBwc4ODggGrVquHjx4+4cuUKXr16BXd3d+zatQvNmjWTNG8KBAI8efIEMTExSE9PlxxXKBRi5cqV6N+/Pzp16qSo0yuRZP3Jo0HITi5+KhJGNohI0knfaYoTW9WAYRRE6kQMAJo0aYJt27bBz88Pr1+/hp+fH7Zu3YomTZrI4vAMU2p16tTBwYMHi4wWVBU5OTnw9PRE8+bN4ebmhqysLFy5cgVqamro378/fvjhBwwZMgRdu3aFmpoaunXrhhYtWkgev2LFCuzbtw9WVlbYuXOnZHtkZCRGjRqF9evXK+CsSs+mow2qN60OYbYQTw8+VXQ4lVrkw0jEB8dDQ1cDLUa3UHQ4DFNl8YiKGWLGyEX+pKGpqakwNDQsdp+cnBy8f/8ederUgbY2W+uNUT7yfo/6bvHFlRlXYOtqizHeY2R+fIZzY9EN3F15F02GN8GgvwYpOhyGUWql+f0uL5nUiDEMw8hKvT7cxLoR9yKQm56r4Ggqr3dXualg7HvZl7AnwzDyxBIxhmGUiqm9KUzsTCAWivHh1gdFh1MpZSVmITogGgBg34MlYgyjSCwRYxhG6dj35JKDt1ffKjiSyund9XcAAdWbVodBDYOSH8AwjNywRIxhGKUjScT+Y4mYPOQnuPnPM8MwisMSMYZhlE4dtzrgq/OR9CYJye+SFR1OpUJEkgSXNUsyjOKxRIxhGKWjZaiFWu1qAWDNk7KW8CIB6dHpUNdWh03HkpfjYhhGvlgixjCMUsqvrWGJmGzl14bZutpCQ0dDwdEwDKMu7QHy8vJw9uxZ+Pv7IyUlBSKRqMg+PB4P+/btk7YohmGqEPue9ri1+Bbe33gPkUAENQ01RYdUKUj6h7FmSYZRClIlYmFhYejevTvevn2Lr80LyxIxhmHKqkarGtAx1UF2UjaifKNYM5oMCLIFCLsdBoB11GcYZSFV0+Ts2bPx5s0bjBo1Crdu3UJoaCjev39f5Pbu3btylxEVFYWNGzeiR48esLGxgaamJiwtLTFo0CD4+vqW+7gCgQAtWrQAj8eDg4NDsfvUrl0bPB6v2NvkyZPLXTbDMCXjq/Fh190OAGuelJXwu+EQ5ghhUNMA5o3MFR0OwzCQskbs5s2b6Nq1Kw4ePCireIrYsmUL1qxZA3t7e3Tv3h3Vq1dHaGgozp49i7Nnz+Kvv/7C0KFDy3zcZcuW4c2bNyXuZ2RkhFmzZhXZ7uTkVOYyGYYpG/ue9nh+4jne/vcWbkvdFB2Oyvt8tGT+QvEMwyiWVImYWCxGy5YtZRVLsZydneHj4wMXF5dC2+/cuYOuXbvi559/xrfffgstLa1SHzMwMBCrVq3C+vXrMWPGjK/ua2xsDA8Pj/KEzjCMlOy7c81nUf5RyE7Kho6pjoIjUm2sfxjDKB+pmibbtWuHkJAQWcVSrIEDBxZJwgDAxcUFbm5uSEpKQlBQUKmPl5eXhzFjxqBt27aYNm2aLENlGEbGDGsZwryxOUCfZoNnyi09Oh3xQfEAD7DrZqfocBiG+USqRGz16tW4desW/v77b1nFUyYaGtzQa3X10lfseXh4IDQ0FPv27StV1Xxubi4OHjyIlStXYseOHXj69Gm542XkTywWo3nz5ujTp4+iQ5GJN2/eQF1dHdu3b1d0KArDprGQjbfXuOfPytEKuma6Co6GYZh8UjVNXrhwAW5ubhg2bBhcXV3RsmVLGBkZFdmPx+Nh8eLF0hRVRHh4OK5fvw5LS0s0bdq0VI/x9/fH2rVrsXLlStSvX79Uj4mNjcWYMWMKbevVqxcOHz4MMzOzrz42NzcXubm5kr/T0tJKVSZTfl5eXnj27Bn27NlT5L4jR47gzp07CAgIQFBQEPLy8nDgwIEir29pNWvWDEFBQQgLC4ONjXxG9NWtWxcjR46Eh4cHRo0aBUNDQ7mUo8zse9rj4YaHePsfNzqb9W0qn3dXuRpFNlqSYZQMSYHH45XqxufzpSmmiLy8POrUqRMBoEOHDpXqMTk5OdSoUSNycnIioVAo2Q6AGjRoUOxjPD09ydvbmxISEigtLY0ePnxIvXv3JgDUrl07EovFXy3T3d2dABS5paamfvEx2dnZ9OLFC8rOzi7VeTEFhEIh2djYkKura7H329raEgAyMzOT/P/AgQPlKis7O5vU1dXJzMys/AGXUlBQEAGgZcuWyb2s0qjo92heZh4t01pGHvCg+BfxFVJmZSMWiWmt2VrygAd9uP1B0eEwjMpJTU0t8fe7vKSqEbt165Y0Dy8XsViMcePGwcfHBxMmTMAPP/xQqsctXrwYoaGhCAgIgJpa6SaGXLJkSaG/27Rpg4sXL8LV1RV3797F5cuX0bdv3y8+fsGCBZgzZ47k77S0NFhbW5eqbKbsLl++jPDw8CKvW769e/eiXr16sLW1xerVq7FgwYJyl/X06VMIhUK0atWq3McorSZNmqB58+bYs2cPFi5cCD6/ai2IoaGrAdtOtnh37R3e/vcW5g3ZtAtlFfskFlmJWdDU10SttrUUHQ7DMJ+R6hvd1dW11DdZICJMmDABR44cwahRo7Bz585SPS4wMBDr16/HokWLSt2M+SV8Ph9jx44FANy7d++r+2ppacHQ0LDQjSmbjIwMeHh4oEGDBtDW1kbdunUlzY43btwAj8fDpUuXAHDNkjweD4MGDSr2WN26dYOtra1M4goMDAQAODo6yuR4JRk6dCjCw8Nx48aNCilP2Uj6if3H+omVx5v/uKl66nSpAzVNtkIBwygTmVxa379/HxMnToSzszMaNGiA1q1bY+LEiSUmKmUhFovx008/Yf/+/RgxYgS8vLxKXTPw7NkziEQieHh4FJmYFQBevXoFHo8HY2PjUh0vv29YVlZWuc6FKZ3IyEg4OTlhxYoVcHJywpQpU5Ceno5Jkybh8ePHWLFiBVq3bo2+ffuCiODt7Q0HB4dSv47SCAgIAIAKqREDuBHKADd3X1WU36/pw+0PEOYIFRyN6snvH2bXg42WZBhlI/Vak3PnzsWGDRskSxzx+XyIxWIEBARg3759mDlzJtavXy9VGWKxGOPHj8eBAwcwbNgwHD58uNTNiwBQv359/PTTT8Xet2/fPhgZGWHw4MHQ1S3dSKL8Gf1r165d6hiYshGLxRg0aBBevXqFc+fOoX///gCAPn36oHv37li7di1u3bqFy5cvAwBCQkKQlJSE3r17V0h8FV0jlj+B8P379yukPGVTvUl16NfQR0ZMBsLvhrPpF8ogLyMP4ffCAQB1e9ZVcDQMwxQhTQezgwcPEo/Ho4YNG9Lx48cpNjaWiIji4uLoxIkT1KhRI+Lz+XTw4MFylyESiWjMmDEEgIYMGUICgeCr+yckJFBISAglJCSU6vj4Qmf958+fU3JycpHtd+7cIW1tbdLS0qKwsLBSlZGvNJ39StMROiMjgzIyMgoNFsjNzaWMjAzKyckpdl+RSCTZlpeXRxkZGUXKKMu+mZmZlJGRUWjgQ0mvTVmcOXOGANCwYcMKbU9OTpYMemjTpo1k+3///UcAaM6cOaU6/qpVq8rdWT83N5c0NDTIxMSkzI+Vhra2NtnZ2VVomcVR1ICSM6PPkAc86Oq8qxVarqp7deEVecCDNtbZWOIAI4ZhiifPzvpSNU3u2LED1tbW8PX1xbBhw2BhYQEAqF69OoYOHYoHDx6gVq1aUs2BtHTpUnh5eUFfXx/169fH8uXL4eHhUej25MkTyf5bt25Fw4YNsXXrVmlODSdPnoSVlRX69euH6dOnY+7cuejVqxc6deoEgUCArVu3ym3KgpLo6+tDX18fiYmJkm3r1q2Dvr5+kUlqq1evDn19fYSHh0u2bdu2Dfr6+kVqCWvXrg19ff1Ck/TmP/fDhw8vtG+jRo2gr68vqRkCgBMnTsjk/ADg2LFjAFBk5QNNTU3J/z09PSX///jxIwDAxMREZjF8ybNnzyAQCGTeLDl+/Pgv1twCgKmpaaHXvKrJb55k/cTKJr9/mH1PtqwRwygjqZomg4ODMWHCBBgYGBR7v6GhIQYOHIi9e/eWu4wPHz4A4Dptr1ixoth9ateujRYtWpS7jOK4ubkhJCQEgYGBuH37NnJycmBhYYFhw4Zh9uzZcHZ2lml5TGE+Pj4wMTFB27ZtC22nT03g7dq1Q8+ePSXbdXS4pW+ys7PlHpu8miX9/f2/uph8dnZ2qZvPKyO7bnYAD4h7Fof0mHQY1Cj+e4cpTDJ/GFvWiGGUktSd9fN/GL9E2iswLy8vENFXb59PyOnh4QEiKvX6kESEly9fFtnu6uqKEydOIDQ0FGlpacjLy0NERAT++usvhSdhGRkZyMjIKDSh7Lx585CRkVGkJjA+Ph4ZGRmFau+mTp2KjIwM7Nu3r9C+Hz58QEZGBho2bCjZNmbMGGRkZOD48eOF9n3x4gUyMjIK1QoNGzZMJueXmpqKuLg4NGjQoMiAjPw+Yd98802h7ebm3JQGSUlJMonha4rrqC8UCuHh4QF7e3toa2ujZs2aWLhwoeT+BQsWoGHDhtDV1UXNmjUL1ebl5ORATU0Nz549w5QpU8Dj8Yo8l2KxGKmpqZLzrIr0zPVQo1UNAMC7a2y5o9JI+ZCCj68/gqfGQ50udRQdDsMwxZCqRqxJkyY4ffo0li1bBn19/SL3p6en4/T/2rvzqKaO9g/g3xAg7IuggIIIqFBRcEEQFRGrQrVS64JVXKhatVpt3WtflaCir/5arbZ1aVGhorW2WrUt7opW3HBFwSKLgiwqiIYdApnfH7xJTcMSSCAEns85nANz7515bm4gD3Pnzhw+DGdnZ0WaIf+ir68vU6atrS112662fbW0tCTLQzV03+p6Zuqz1FRtxE+j/vuBjNLSUnz++efVtuXs7AwNDQ0kJSUpJYbaVNcjtnbtWpw4cQIRERGwsbFBamoqcnJyAABCoRA8Hg/h4eGwtLTE9evXERQUBC8vLwwZMgTa2to4e/Ys3n77bSQlJUFfX1/mWiQlJUEkEik8/Yq6cxjugOxb2Ug5lQLXqa6qDqfZEy8LZd3PGjrGOiqOhhBSHYV6xObMmYOMjAx4enri8OHDkvErubm5+PXXX9G/f39kZGTg448/VkqwpHVo27YtdHR0cPv2baSlpUnKFy5ciNTUqp6QFy9eSB1jYmICFxcX3Lx5s85eWkUIhULcv38fRkZGcHD451bPmTNnMGbMGAwcOBC2trbw8fFBQEAAgKpkls/nw8PDA7a2tggICECvXr2QmJgIoOpJ4+zsbHTq1AkODg6wtLSUud0vflJXWXPyqSvJOLEzKWCixrvOLYU4EaNljQhpvhRKxKZNm4ZPP/0U8fHxCAgIgIWFBbS0tCRjqeLj4/HJJ59g2rRpyoqXtAKampoIDAxESUkJvLy88Omnn2LIkCHYuXMn1qxZAz09PezevRsrV65EUVGR5LjRo0dDIBAgNja22nrDwsIQFBSEoKAg/PLLLzJlR48erTO2+Ph4lJWVoXfv3lK33UeOHImVK1fC398fkZGRKCwslGxLTU3F7Nmz0a1bN5iYmMDAwADXrl2DlZWVZJ+4uDi4utbcw3PmzBlwuVyZW7KtjY2nDbQNtFGcU4xnd5+pOpxmTVQhQupZGh9GSHOn8BixLVu24NKlSwgKCkLPnj0lA+c//PBDXLx4EVu3blVGnKSV2bp1K+bNmwehUIgdO3YgKSkJW7duxapVq/Ddd99BU1MT27Ztk7pFOnPmTHC5XERGRlZb5+XLlxEREYGIiAjJ7cWYmBhJ2ZtP39akpolc//Of/yAuLg4eHh4IDg5Gt27dIBAI8OLFC7i7u6O0tBTffPMNrl69it9//x0ikQguLi6S4+/duyf185uKi4tx9OhRjBo1Cu3bt68zxpaMq81FJ59OAP7p7SHVy4zNRJmgDDqmOmjv1rrfN4Q0ZxzWmPdxiJT8/HwYGxtDIBDUuNxRaWkpHj9+DDs7O+jo0JiO+po0aRJOnz6NtLS0ase8NYWsrCx06NABDx48wLVr17Bq1SpkZWVJti9fvhw7d+7E69evJb1q1tbW+PrrrzFu3DiZ+vbs2YMZM2bg4sWLGDRoUJOdR01U/R698d0NnPjkBDoN7oRpF6i3vSbR/GhcDLmIbuO7Yfyh8aoOhxC1Js/nd0O1rtWDSYsXGhqKwsJCfPfdd03W5saNG7F//34kJiYiISEBq1atgpOTE5ycnGBmZoaXL1/i5MmTSExMRHBwMHbt2gUXFxepW5tCoRB3795FdnY2CgoKJOUVFRVYv349/P39m0US1hyIb7Olx6SjvLBcxdE0X5LxYXRbkpBmjRIx0qLY2dkhIiKiSXvDSktLERISAldXV/j4+KC4uBgnT54El8uFv78/pkyZgvHjx+Ptt98Gl8vF0KFDZea9Cw0Nxe7du9G+fXupxewzMjIwefJkhZcJa0nadG4DEzsTiIQiPIl+oupwmqWSVyXIvJ4JgBIxQpq7et2anD59OjgcDtavXw8LCwtMnz5dvkY4HJk5q1ojujVJWoLm8B79Y84fuLXrFtznu+OdbU2zvqg6STicgF/G/QJzJ3PMezhP1eEQovYa89ZkvSZ+Cg8PB4fDwfLly2FhYYHw8HC5jqNEjBCiTA6+Dri16xYtd1QD8etC01YQ0vzVKxF7/PgxAKBDhw5SPxNCSFOyG2IHDpeDl49e4vWT1zDpZKLqkJoNxhiNDyNEjdQrEbO1tZX6mcPhwMTEpNZuuoKCArx69aph0RFCSDV0jHVg3c8aT2OeIuV0CvrMUu66n+rs5aOXEKQJwNXmwtbbtu4DCCEqpdBgfTs7uzrnCdu+fTvs7GiNM0KIcklm2afbk1LEvWEdB3aEtr7ssmeEkOZFoURMvOh2XfsQQoiyiW+7pZ5LhahCpOJomg8aH0aIemn06SsyMjJk1s0jhBBFtXdrDx1THZQJypB5I1PV4TQLFWUVeHLhCQAaH0aIuqjXGDEAWLNmjdTP0dHR1e5XWVmJjIwMHDx4EB4eHg0KjhBCaqLB1YD9UHsk/JKAlNMpsOlvo+qQVC7jagaExULot9OHhYuFqsMhhMih3okYn8+XfM/hcBAdHV1jMgYA7du3x8aNGxsSGyGE1MrB16EqETuVgsH8waoOR+WSTyUDqOoN42hw6tibENIc1DsRu3DhAoCqsV9DhgxBUFAQpk2TXe+Ny+WiTZs2cHJygoYGTeBPCFE+8e23zBuZKHlVAl1TXRVHpFqpp1MBAPbD7VUcCSFEXvVOxLy9vSXfBwcHw8fHh9bAI4SohLGNMczfMkfuw1w8PvcY3cZ1U3VIKlP0ogjZt7MBAA7DaHwYIepCoa6q4OBgSsIIISol7hUT35ZrrVLPVvWGWbhawMDSQMXREELkpbR7hhUVFUhISMDVq1eRkJCAiooKZVVNCCE1Ek/TkHo6tVVPl0PTVhCinhROxHJycvDRRx/BxMQEPXr0wMCBA9GjRw+YmJhg1qxZyMnJUUachBBSrU7encDV5kKQLsDLxJeqDkclaFkjQtSXQolYZmYm+vbti927d0NfXx++vr6YOnUqfH19oa+vj7CwMLi7uyMzk+b4IU1DJBLB1dUVI0aMUHUoKpGcnAxNTU1s375d1aE0GS09LXT06gjgn1nlW5sX91+g8FkhNHU10XFgR1WHQwipB4USsWXLliE9PR0hISFIS0tDVFQU9u7di6ioKKSlpYHP5yMtLQ3Lly9XVryE1Co8PBxxcXFS06yIderUCRwOp9qvOXPmVFtfbGwsRowYAVNTU+jr68Pd3R0HDhxoUGwuLi7gcDhIT09v0PHy6Ny5MwIDA8Hn85Gfn99o7TQ3rX25I3EC2mlwJ2jy6v0MFiFEhRT6jT158iT8/PywatUqmW06OjpYvXo1rly5ghMnTijSDCFyqaysREhICLy9veHu7l7tPsbGxvjss89kyt3c3GTKoqOj4evrC21tbXzwwQcwNjbGkSNHEBgYiCdPnuCLL76QO7bS0lI8fPgQ5ubm6NixcXssli5dih9//BHbtm3DypUrG7Wt5sJhuAPOLjuLJ9FPUFFW0eqSERofRoj6UuivVXl5OXr37l3rPn369EFMTIwizRAil6ioKKSnp2P16tU17mNiYlJtb9m/VVRUYObMmeBwOLh06RJ69eoFoOpJYU9PTwQHB2P8+PHo0qWLXLHdu3cPFRUVdf6+KEP37t3h6uqKH374AV988UWrmMfPwsUC+hb6KHpehKcxT2E3xE7VITUZYbEQaX+lAaDxYYSoI4X+Qvfp0wd///13rfv8/fff6NOnjyLNkFassLAQfD4fjo6O0NHRQefOnfHDDz8AAM6dOwcOh4M///wTQNVtSQ6Hg7Fjxyrc7vnz55GSkoJJkyZJkjAAMDQ0xKpVq1BRUYG9e/fKXd/t27cBoMl+FwICApCeno5z5841SXuqxuFwWu00Fml/paGyrBJG1kYwdzJXdTiEkHpSKBFbu3Yt/vjjD4SHh1e7fc+ePYiKisK6desUaYa0UhkZGXBzc0NoaCjc3Nwwd+5cFBQUYPbs2bhz5w5CQ0PRt29fjBw5EowxREdHw8nJCSYmJjXWWVZWhoiICKxfvx47duzAvXv3qt1PvGzX8OHDZbaJyy5evCj3udy6dQsAmqRHDAA8PT0BVCWUrcWb01i0Jm/eluRwaFkjQtSNQrcmL1y4AB8fH8yYMQObNm3CgAED0K5dO7x48QIxMTFITEzE8OHDcf78eakPBA6HU+24MkLERCIRxo4di8TERBw7dgz+/v4AgBEjRmDYsGHYtGkTLly4gKioKADAw4cPkZeXh3feeafWep89e4agoCCpMj8/P+zbtw/m5v/0JiQlJQFAtbceTU1NYW5uLtlHHk3dIyYe83blypUmaa85EM8m/+zuMxQ+L4SBReuY1JSmrSBEvSmUiL051ubvv/+u9jblqVOncOrUKakySsQahjGGYmGxqsOok56WnsL/mR8/fhw3btzAhAkTJEkY8E+CcfDgQXh4eEgSr4yMDACAhYVFjXVOnz4d3t7ecHZ2Bo/HQ0JCAkJCQnDixAn4+/sjJiZGErdAIABQNbi/OkZGRpI261JeXo4HDx7A1NQUdnZNM3bJ0NAQOjo6csfYEui304dlL0s8u/MMqWdS4TLZRdUhNbr8jHzkxOcAHMB+KK0vSYg6UrhHjDSdYmExDDY0///yC1cUQl9bX6E6xFNELFiwQKpcW1tb8n1ISIjk+5cvqybyNDU1rbHOfw/i9/DwwB9//AFvb29cvnwZUVFRGDlypEJxVycuLg5CoVDptyVnzpwJxhh2795d7fY2bdogNzdXqW02dw6+Dnh25xlSTqW0ikQs5UxVb1iHvh2g26Z1L3hOiLpSKBF7cwFwQpTp0qVLMDU1Rb9+/aTKxUvYeHp6wtfXV1Kuq1v1IVRSUlKvdjQ0NPDhhx/i8uXLiImJkSRi4p4wcc/Yv+Xn59fYW/ZvjXVbMjY2tsb5z4Cq10JPT0+pbTZ3DsMdEPPfGKScSQETMXA0WvaYKZq2ghD1p1AixuVy8cEHH2D//v3KiofUQk9LD4UrClUdRp30tBT78BcIBHj+/Dn69esnM/WCeEzYu+++K1Xetm1bAEBeXl692xOPDSsu/ue2r3hsWFJSkkwC9erVK+Tm5qJ///5y1V/TQP2KigqsW7cO+/btQ2ZmJszMzDBt2jSsX78eALBixQocPXoUaWlpMDU1xaxZsxAcHIzS0lLo6+tDJBJh7ty5mDt3LgICAvDzzz9L6haJRBAIBHB2dq7nq6HebPrbQEtfC0XPi/A87jkse1qqOqRGI6oUIfVM1YMJND6MEPWlUCJmZGQEGxsbZcVC6sDhcBS+5acOxAkRl8uVKi8tLcXnn38OANDUlH7rOjs7Q0NDo14D6MWuX78OoGrmfTFvb29s2LABp0+fxgcffCC1/+nTpyX7yKOmHrG1a9fixIkTiIiIgI2NDVJTUyVrswqFQvB4PISHh8PS0hLXr19HUFAQvLy8MHjwYJw9exZvv/02kpKSoK+vD3196fdFUlISRCIRevToIf8L0QJo8jTRaXAnJP2ZhORTyS06Ecu+nY2SvBLwjHjo4NFB1eEQQhpIoekr3N3da3z8n5CGatu2LXR0dHD79m2kpaVJyhcuXIjU1KoegBcvXkgdY2JiAhcXF9y8eVNy+/JNCQkJeP36tUz55cuXsXnzZvB4PIwZM0ZS/vbbb8Pe3h4HDhzA3bt3JeUFBQVYu3YtNDU1ZZ6+rI5QKMT9+/dhZGQEBwfpXoszZ85gzJgxGDhwIGxtbeHj44OAgAAAgJaWFvh8Pjw8PGBra4uAgAD06tULiYmJ0NDQQHZ2Njp16gQHBwdYWlrC0NBQqm5xctkahw+0lmksxE9L2g2xA1eLW8fehJDmSqFELCQkBOfPn0dERISy4iEEmpqaCAwMRElJCby8vPDpp59iyJAh2LlzJ9asWQM9PT3s3r0bK1euRFFRkeS40aNHQyAQIDY2VqbOQ4cOoX379hg1ahTmz5+PJUuWwM/PD4MGDYJQKMS3334rtfSQpqYmwsLCIBKJ4OXlhVmzZmHJkiVwdXVFfHw8+Hw+unbtWue5xMfHo6ysDL1795Z5knTkyJFYuXIl/P39ERkZicLCf247p6amYvbs2ejWrRtMTExgYGCAa9euwcrKCkDVAwCurq41tnvmzBlwuVyZW7itgfg2XfrldJQXlas4msYjHh9mP5yeliREnSl0a/L06dMYPHgwpk+fjm+++Qbu7u6wsLCQ+cCh6SpIfW3duhU6Ojo4fPgwduzYAQsLC2zduhULFiyAjY0Nli5dim3btmHt2rWSY2bOnIm1a9ciMjJSZq1JHx8fPHz4ELdv38bFixdRWloKCwsLTJgwAQsXLqx2bUofHx9cvnwZwcHBOHToEMrLy+Hs7Iy1a9ciMDBQrvOobSLX//znP3j//ffx22+/ITg4GF988QXu37+PsrIyuLu7Y+TIkfjmm2/Qvn17PHv2DEOGDIGLS9WTgPfu3atxPc3i4mIcPXoUo0aNQvv27eWKsyUx62oGY1tjCNIESLuYhi4j5FuGSp2U5Zch42rV1CSdfTurOBpCiCI4rLr7OHKSdw07DoeDysrKhjbTYoiftBMIBDAyMqp2n9LSUjx+/Bh2dnbQ0dFp4gjV36RJk3D69GmkpaXJjJtqzrKystChQwc8ePAA165dw6pVq5CVlSXZvnz5cuzcuROvX78Gh8OBtbU1vv76a4wbN06mrj179mDGjBm4ePEiBg0apPRY1eE9+vvs33H7+9twX+COd7bWPsmvOko8noiD7x2EqYMpFiQvqPsAQohC5Pn8biiaR4y0KKGhoThy5Ai+++47LFu2TNXh1Gjjxo2wtraGm5sbKisr8dVXX8HJyQlOTk5ISkrCy5cvcfLkSdjZ2eHAgQPYtWsXXFxcJL3NQqEQd+/exYABA2BgYCAZI1ZRUYH169fD39+/UZIwdeEw3AG3v7/dYseJidfTpGkrCFF/NI8YaVHs7OwQERHR7CcyLS0tRUhICNLT02FsbIwhQ4bg5MmT4HK58Pf3x5QpUzB+/HgYGxtj1qxZGDp0qGR8GFCVcK5atQqhoaHYtGkTli5dCqBqhYHJkydjypQpqjq1ZsH+bXtwNDjI/TsXgnQBjDvKN+ebuhAnmDRtBSHqT6Fbk6R+6NYkaQnU5T26u/9uZFzNwKgfRqH3zKZZbL0pvEp9hW0O26ChqYFlL5eBZ8RTdUiEtHiNeWtSoacmxa5cuYJZs2bB3d0djo6O6Nu3L2bNmoXLly8ro3pCCKk38W078dOFLYV42gprT2tKwghpARROxJYsWQIvLy+EhYXh5s2bSElJwa1btxAWFgZvb28sWrRIGXESQki9iG/bpZ5NhahCpOJolIeWNSKkZVEoEfvxxx+xefNmODo64qeffkJ2djYqKirw7NkzHDx4EE5OTti6dSt+/PFHZcVLCCFy6dC3A3RMdFD6uhRZN7PqPkANVAor8fj8YwA0PoyQlkKhRGzHjh2wsbHB9evXMWHCBFhYWAAA2rVrh4CAAFy9ehXW1tbYvn27UoIlhBB5aWhqwH5o1WSn4qcM1V3m9UyU5ZdB10wXVr2t6j6AENLsKZSIPXjwAGPHjpVZXkXMyMgIY8aMQXx8vCLNEEJIg4hnnW8p01iIx4fZD7WHBlcpQ3wJISqm8G9yXQ9d/nuWfUIIaSri23cZ1zNQ+rpUxdEoTjI+jG5LEtJiKJSIde/eHYcPH5ZaI+9NBQUFOHz4MJydnRVphhBCGsTE1gRmjmZglUwytkpdleSVIDM2EwAlYoS0JAolYnPmzEFGRgY8PT1x+PBhySSaubm5+PXXX9G/f39kZGTg448/VkqwhBBSX+KnC9V9nFjquVSAAW27tYWRtXLnMSKEqI5CM+tPmzYNd+/exdatWxEQEACgav1JkajqUXHGGObPn49p06YpHikhhDSAw3AH3Nh2AymnUsAYU9vhEjRtBSEtk0KJGABs2bIFY8eOxd69e3H37l3k5+fDyMgIvXr1wrRp0+Dl5aWMOAkhpEE6De4EDS0NCNIEyEvKg1lXM1WHVG+MMclAfbotSUjLonAiBgADBw7EwIEDlVEVIYQolba+NjoO7IgnF54g+VSyWiZiuX/nIv9pPrg8LmwH2ao6HEKIEjXK88+MMSQlJSEjI6MxqieEkHoR385T12ksxL1htl620NLTUnE0hBBlUigRO3bsGKZPn45Xr15Jyp48eYIePXrAyckJtra2CAwMlIwZI4QQVejs2xkA8PjCY1SWV6o4mvqj8WGEtFwKJWI7d+5EbGwsTE1NJWWfffYZEhIS4OPjAxcXFxw8eBB79+5VOFBC5CESieDq6ooRI0aoOhS1kJycDE1NzRa/+oWFiwX02+lDWCTE0ytPVR1OvVSUVeBJ9BMAND6MkJZIoUQsPj4e7u7ukp8FAgGioqIwYcIEnD17Fjdu3MBbb72F3bt3KxwoIfIIDw9HXFwc+Hy+zLbIyEjMnj0bbm5u4PF44HA4CA8Pr7W+2NhYjBgxAqamptDX14e7uzsOHDig9GOq4+LiAg6Hg/T09HofK6/OnTsjMDAQfD4f+fn5jdaOqnE0OJIkRt2msUi/nI6KkgoYWBqgXY92qg6HEKJkCiViOTk5sLL6Z72zy5cvo6KiAhMnTgQAaGlpYdiwYUhOVq8/fEQ9VVZWIiQkBN7e3lL/IIitXLkS33//PdLS0qTetzWJjo7GwIED8ddff2HcuHH4+OOPkZubi8DAQKxfv15px1SntLQUDx8+hLm5OTp27Cj3cQ2xdOlS5OTkYNu2bY3ajqqp63JHbz4tqa5TbxBCaqZQImZkZISXL19Kfo6OjoaGhobUlBVaWlooKipSpBlC5BIVFYX09HRMmTKl2u1hYWF48uQJcnJyMGfOnFrrqqiowMyZM8HhcHDp0iX88MMP+PLLL3Hv3j04OzsjODgYSUlJCh9Tk3v37qGiogK9e/eW7+QV0L17d7i6uuKHH35o0eM5HYZV9Yhl385GfqZ69P4xxvDo+CMA/ySShJCWRaFEzMnJCb///jvy8vIgEAhw8OBB9O7dW2rMWFpaGiwsLBrcRmZmJr7++msMHz4cHTt2hLa2NiwtLTF27Fhcv369wfUKhUL07NkTHA4HTk5ONe6nrNtMpGEKCwvB5/Ph6OgIHR0ddO7cGT/88AMA4Ny5c+BwOPjzzz8BVN2W5HA4GDt2bLV1DR06FLa28j36f/78eaSkpGDSpEno1auXpNzQ0BCrVq1CRUWFzNjHhhxTk9u3bwMA+vTpI9f+igoICEB6ejrOnTvXJO2pgoGlATp6VfUuxn4Xq+Jo5JNyOgW5f+dCS18LXUd2VXU4hJBGoFAitmDBAmRlZaFDhw6wsbFBVlaWVE9DZWUlLl++DFdX1wa38c0332DhwoVITU3FsGHDsHjxYgwcOBDHjh1D//79cejQoQbVu3bt2jpvmSrrNhNpmIyMDLi5uSE0NBRubm6YO3cuCgoKMHv2bNy5cwehoaHo27cvRo4cCcYYoqOj4eTkBBMTE4Xbjo6OBgAMHz5cZpu47OLFiwofU5Nbt24BQJP0iAGAp6cngKpksiXzXFx1njd33ER5YbmKo6nb1S+vAgB6f9QbOiY6Ko6GENIomIK2b9/O+vTpw/r06cM2btwote3UqVPMxMSE7dy5s8H1Hz58mF26dEmm/NKlS0xLS4u1adOGlZaW1qvOW7duMU1NTbZt2zYGgDk6OsrsIxQKmYODA+PxeOz27duS8vz8fObs7Mw0NTXZo0eP6tWuQCBgAJhAIKhxn5KSEpaQkMBKSkrqVXdLU1lZydzd3RkAduzYMUn5mTNnGAD2wQcfMAAsKiqKMcZYfHw8A8ACAwPlqn/Dhg0MANu7d2+128eNG8cAsJs3b1a73dzcnLVt21bhY2rSq1cvBoClpqbKtb+i8vPzGQA2aNCgOvdV5/eoqFLEtnXZxvjgs2tbr6k6nFpl381mfPBZiEYIe/X4larDIaRVk+fzu6EUntD1448/xs2bN3Hz5k0sW7ZMatvw4cPx6tUrzJ49u8H1jxkzptplkry8vODj44O8vDzcv39f7vrKy8sRFBSEfv364ZNPPqlxP2XeZlIWxoCioub/xZji53r8+HHcuHEDEyZMgL+/v6Tczc0NAHDw4EF4eHjgnXfeAQDJ5MGK3AZ/k0AgAAAYGxtXu93IyEiyjyLHVKe8vBwPHjyAqakp7Ozs6hN2gxkaGkJHR6fFT8LM0eBIesWubbkGUUXzHRN39auq3rBu47vBpJOJaoMhhDQapSxxBFQNVH706BEEAgGMjY3RtWtXaGoqrfpqaWlVzTBdn3b4fD6SkpJw7969Wp9AUuZtJmUpLgYMDJq0yQYpLAT09RWrQzwOb8GCBVLl2traku9DQkIk34sfGnlzfKK6iouLg1AobJTbkjNnzgRjrNopZdq0aYPc3Fylt9ncuE51xYWVF/D6yWs8PPIQzgHOqg5JRn5GPh789AAA0H9JfxVHQwhpTAr3iOXk5OCjjz6CiYkJevTogYEDB6JHjx4wMTHBrFmzkJOTo4w4ZaSnp+Ps2bOwtLREjx495DomNjYWmzZtQkhICLp2rX3gq/jpti5dushsMzU1hbm5udxPwJH6u3TpEkxNTdGvXz+pcva/7jZPT0/4+vpKynV1dQEAJSUlSmlf3KtVUw9Wfn6+TM9XQ46pTmMO1I+NjZX0Kv5bSUkJ9PT0lN5mc6Olq4W+n/QFAFz58orkPdWcXN92HaIKEWy9bdHerb2qwyGENCKFErHMzEz07dsXu3fvhr6+Pnx9fTF16lT4+vpCX18fYWFhcHd3R2ZmprLiBVD1xOOUKVNQVlaGTZs2gcvl1nlMWVkZgoKC0KtXLyxevLjO/ZVxm6msrAz5+flSX4rQ06vqbWruX4p+lgsEAjx//hyOjo7Q0JB+i0ZFRQEA3n33Xanytm3bAgDy8vIUa/x/xAl4dcn2q1evkJubK5OkN+SY6tQ0UL+iogJ8Ph8ODg7Q0dFBhw4d8MUXX0i2r1ixAm+99Rb09PTQoUMHqR7D0tJScLlcxMXFYe7cueBwOJgwYYJku0gkgkAgkLyOLV3fuX2hqaOJrNgspP/VeBPmNkRZfhlu7ap6D1BvGCEtn0KJ2LJly5Ceno6QkBCkpaUhKioKe/fuRVRUFNLS0sDn85GWlobly5crK16IRCJMnz4dly5dwkcffVTjnFH/tmrVKiQlJWHPnj1yJW7KsGHDBhgbG0u+bGxsFKqPw6m65dfcvxSdc7K4uBgAZK5TaWkpPv/8cwCyt6OdnZ2hoaGhtF5Kb29vAMDp06dltonLxPsockx1auoRW7t2LaKiohAREYHExERERkaiZ8+eAKr+OeHxeAgPD8fDhw+xZcsWbNy4UfIUpLa2Ns6ePQsOh4Pk5GRkZ2cjLCxMUndSUhJEIpHcvcvqTr+tPlyDqp7mvvJ/V1QcjbTbu2+jLL8MZo5m6DKi7sSdEKLmFBnp36ZNG/bOO+/Uuo+vry9r06aNIs1IiEQiNn36dAaATZ48mVVWVsp13K1btxiXy2V8Pl9mG2p4alIZT8CVlpYygUAg+Xr69Ck9NSkHoVDIdHR0mK6uLnvy5ImkfM6cOQwAA8AWL14sc1zPnj2ZsbExE4lEdbZR11OTQqGQ2dvbMx6Px+7cuSMpf/Op2cTERIWP+bfy8nLG4/GYkZGRzHl4enqyDRs21HluYv3792fbt2+X/Lx//35mZ2dX7b4REREMANu1a1ed9baU92huYi7jc/iMDz57kfBC1eEwxhirKK9gWzpuYXzw2c3vq//bQwhpes32qcny8vI6BxT36dMH5eWKz9cjEokwY8YM7NmzBxMnTkR4eLjMbauaxMXFobKyEnw+HxwOR+oLABITE8HhcKTmn1LGbSYejwcjIyOpL1I3TU1NBAYGoqSkBF5eXvj0008xZMgQ7Ny5E2vWrIGenh52796NlStXSq3aMHr0aAgEAsTGVj9ZZ1hYGIKCghAUFIRffvlFpuzo0aNSMYSFhUEkEsHLywuzZs3CkiVL4Orqivj4ePD5fJlxhg055t/i4+NRVlaG3r17yzxMMnLkSKxcuRL+/v6IjIxEYWGhZFtqaipmz56Nbt26wcTEBAYGBrh27ZrUUk5xcXE1zul35swZcLlcmVu+LZlZVzM4vVc1mfPVzVdVHE2VhF8TIEgXQL+dPlynNHz+RUKIGlEki/P29mZjx46tdZ8xY8Ywb29vRZphlZWV7MMPP2QA2IQJE1hFRUW9jo+JiWEzZsyo9gsAMzY2ZjNmzGDz58+XHHPy5EkGgH344Ycy9R08eJABYCtWrKhXHDSPmPwKCwvZvHnzmKWlJdPS0mLW1tZs69atjDHG9u7dy8zNzZmhoaFUr1FGRgbjcrlS1/FN06ZNk/SoVfcVHBwsc8z169eZn58fMzY2Zrq6uszNzY1FRkbWGntDjhELCwtjANiiRYuq3R4fH8/WrVvH7O3tmY2NDXv9+jV7/vw5MzMzY1OnTmVnz55lCQkJ7Pz58wwAS0lJkRzr5+fHVq9eLVNnUVERMzAwYKNHj5Yrxpb0Hk27nMb44LO12mtZQXaBSmMRiURsV+9djA8+i14TrdJYCCHSGrNHTKFE7NKlS4zH49V4e2f37t1MR0eH/fXXXw1uo7KykgUFBTEAbPz48UwoFNa6f05ODnv48CHLycmRq37UMqGroreZ/o0SscY3ceJEZmZmxgoLC1UdSqPKzMxkANiDBw9YWFgYs7Kyktq+bNkymdubHTp0YL/88otMXbt372YA2MWLF+VquyW9R0UiEQvrF8b44LNzK8+pNJbHFx4zPvhsne46VpRTpNJYCCHSGjMRq9dEX2vWrJEp8/HxwYwZM7Bp0yYMGDAA7dq1w4sXLxATE4PExEQMHz4cFy5cwMCBAxvSYYc1a9YgPDwcBgYG6Nq1K9atWyezz+jRoyWDlr/99luEhIQgODgYfD6/QW0C/9xm8vX1hZeXFyZOnAgjIyMcOXIEjx8/xrp16+q8zUSaXmhoKI4cOYLvvvtOZoJhdbZx40ZYW1vDzc0NlZWV+Oqrr+Dk5AQnJyckJSXh5cuXOHnyJOzs7HDgwAHs2rULLi4uUrc3hUIh7t69iwEDBsDAwACGhoaoqKjA+vXr4e/vj0GDBqnwDFWDw+Gg/9L+ODT2EG5uv4mBnw+Etr523Qc2gitfVj000DOoJ/TMW/40IoSQKvVKxGpLbP7++2/8/fffMuWnTp3C6dOnsWrVqnoHBwBPnjwBULX4c2hoaLX7dOrUSZKIKZOPjw8uX76M4OBgHDp0COXl5XB2dsbatWsRGBio9PaI4uzs7BAREdHiJiYtLS1FSEgI0tPTYWxsjCFDhuDkyZPgcrnw9/fHlClTMH78eBgbG2PWrFkYOnSo1PgwoCpJXbVqFUJDQ7Fp0yYsXboUGRkZmDx5stxPH7dEju85wtTBFK9SXuFu+F24z3Nv8hhyEnKQ9GcSwAH6LexX9wGEkBaDw5j8sxkqMpO8PI/tt3TiCT0FAkGNA/dLS0vx+PFj2NnZQUeHFvklzU9LfI/Gbo9F1LwomNqb4pNHn0CDq/Bc1/VyfOZx3Nl9B07vO2HCkQl1H0AIaVLyfH43VL16xBqaTFVUVDToOEIIaQo9g3riwuoLeJX6Cn8f/RvdxnZrsrYLnxUibl8cAJrAlZDWqFH/7UtISMDixYthbW3dmM0QQohCtPS00Hfu/5Y9+r+mXfboxnc3UFleCet+1rDpr9ikz4QQ9aP0RKywsBBhYWHw9PREjx49sGXLFrx+/VrZzRBCiFL1ndcXXB4Xmdcz8fTK0yZps7yoHDe33wQAeC7xbJI2CSHNi9ISscuXL2P69OmwsrLC7Nmzcf36dfTs2RPbtm1DVlaWspohhJBGYWBhANepVZOoXv2yaSZ4vRt+FyV5JTC1N4XTaKcmaZMQ0rzUa4zYvz1//hwRERHYs2cPkpKSwBiDpaUlioqKMHXqVISHhyspTEIIaXyeizxx+4fb+PvY33j56CXMupo1WluiShGubb4GAOi3qF+TPyBACGke6v2bLxKJ8Pvvv2P06NGwsbHB559/jvT0dAQEBODPP//E06dVXfra2qqZi4cQQhrK3MkcXUd1BRhwdUvj9oolHkvEq9RX0G2ji55BPRu1LUJI81XvHjFra2s8f/4cADBgwABMnToVAQEBtI4iIaRF6L+kPx79/gj3wu/BZ40P9NvqN0o74glc3T52U9kksoQQ1at3j9izZ8/A4XCwZMkSHD9+HDNnzqQkjBDSYnT06oj2fdujorQCsd9Vv4C8op5eeYqMqxnganPh/knTTyBLCGk+6p2ITZ48GTo6Ovjyyy9hZWWF8ePH4/jx4zRXGCGkReBwOJL5vGK/i4WwWKj0NsS9YS5TXGBgaaD0+gkh6qPeidiPP/6I7OxsbN++HT169MDhw4fx/vvvw9LSEp988gmuXbvWGHESQkiTeWvMWzDpZILi3GLc+/GeUuvOS87D30erloPzXERTVhDS2jXoMR1DQ0PMnj0bN27cQFxcHObPnw8Oh4Pt27djwIAB4HA4SExMRHp6urLjJYSQRqehqSFZ8/Hq5qsQVYqUVvfVLVcBBnQZ2QVtu7VVWr2EEPWk8PPS3bt3x9dff42srCwcPHgQw4YNA4fDwV9//QV7e3sMGzYMP/30kzJiJYSQJtNrei/omOggLykPj35/pJQ6i3OLcXfvXQCA52LqDSOEKHFCVy0tLQQEBODkyZN48uQJ+Hw+OnbsiHPnzmHy5MnKaoYQQpqEtoE23D52A/DPmC5Fxe6IRUVJBax6W6HT4E5KqZMQot4aZQZBa2trrF69GqmpqTh9+jQmTJjQGM0QIkMkEsHV1RUjRoxQdShqITk5GZqamti+fbuqQ2mW3Oe7Q0NLA09jnuLpVcWWPaoorUDst1VPYXou8QSHw1FGiIQQNdfoUzkPHToUBw4caOxmCAEAhIeHIy4uDnw+X2Zbp06dwOFwqv2aM2dOtfXFxsZixIgRMDU1hb6+Ptzd3et8PzfkmOq4uLiAw+E06ljLzp07IzAwEHw+H/n5+Y3WjroytDKEy2QXAMDVrxSb4DUuMg5FL4pg3NEY3cZ1U0Z4hJAWQKEljghpTiorKxESEgJvb2+4u1c/N5OxsTE+++wzmXI3NzeZsujoaPj6+kJbWxsffPABjI2NceTIEQQGBuLJkyf44osvlHJMdUpLS/Hw4UOYm5ujY8eOch3TUEuXLsWPP/6Ibdu2YeXKlY3aljryXOyJu3vv4uGRh8hLyUMbhzb1roOJmCSR8/jMA1wtrrLDJISoK0aajEAgYACYQCCocZ+SkhKWkJDASkpKmjCyluH48eMMAAsLC6t2u62tLbO1tZWrLqFQyBwcHBiPx2O3b9+WlOfn5zNnZ2emqanJHj16pPAxNbl27RoDwIYPHy7X/opydXVlHTt2ZJWVlXXu2xrfo/tH7Gd88Nmf8/5s0PGJvycyPvhsg9EGViooVXJ0hJDGJs/nd0PRKrOkWSssLASfz4ejoyN0dHTQuXNn/PDDDwCAc+fOgcPh4M8//wRQdVuSw+Fg7NixCrd7/vx5pKSkYNKkSejVq5ek3NDQEKtWrUJFRQX27t2r8DE1uX37NgCgT58+Cp+LPAICApCeno5z5841SXvqxnNJ1ROOd/bcQfHL4nofLx7s32d2H/CMeEqNjRCi3igRI81WRkYG3NzcEBoaCjc3N8ydOxcFBQWYPXs27ty5g9DQUPTt2xcjR44EYwzR0dFwcnKCiYlJjXWWlZUhIiIC69evx44dO3DvXvWTdUZHRwMAhg8fLrNNXHbx4kWFj6nJrVu3AAC9e/eWa39FeXpWJRrnz59vkvbUTafBnWDV2woVJRW4ueNmvY7NupmFtItp0NDUgMcCj0aKkBCirmiMGGmWRCIRxo4di8TERBw7dgz+/v4AgBEjRmDYsGHYtGkTLly4gKioKADAw4cPkZeXh3feeafWep89e4agoCCpMj8/P+zbtw/m5uaSsqSkJABAly5dZOowNTWFubm5ZB9FjqlJU/eIicfIXbminGkaWhoOhwPPJZ44MukIbnxzA/2X9Iemjnx/PsVjw7pP7A4ja1qXlxAijRIxNcIYa5R175RNS09L4Ufzjx8/jhs3bmDChAmSJAz4J2E4ePAgPDw8JIlXRkYGAMDCwqLGOqdPnw5vb284OzuDx+MhISEBISEhOHHiBPz9/RETEyOJWyAQAKga3F8dIyMjSZtiDTmmOuXl5Xjw4AFMTU1hZ2dX5/7KYGhoCB0dHbnia626jeuGc5+fgyBdgLjIOPSeWXdv5esnrxH/SzwAmsCVEFI9SsTUiLBYiA0GG1QdRp1WFK6Atr62QnWIp3tYsGCBVLm29j/1hoSESL5/+fIlgKqep5qsXr1a6mcPDw/88ccf8Pb2xuXLlxEVFYWRI0cqFLcyxMXFQSgUNsptyZkzZ4Ixht27d8tsa9OmDXJzc5XeZkvB1eLC4zMPnF50Gle/uope03uBo1H7PxzXtl4Dq2SwH2oPS1fLJoqUEKJOaIwYaZYuXboEU1NT9OvXT6qcMQagakyTr6+vpFxXVxcAUFJSUq92NDQ08OGHHwIAYmJiJOXiXi1xL9e/5efny/R8NeSY6jTmbcnY2Nhqp+oAql47PT09pbfZkvSe2Rs8Yx5y/87Foz9rX/ao9HUp7oTdAfDPYH9CCPk36hFTI1p6WlhRuELVYdRJS09LoeMFAgGeP3+Ofv36QUND+n8F8Ziwd999V6q8bduqxZPz8vLq3Z54bFhx8T9Pw4nHeSUlJckkRK9evUJubi769+8vVd6QY6pT00D9iooKrFu3Dvv27UNmZibMzMwwbdo0rF+/HgCwYsUKHD16FGlpaTA1NcWsWbMQHBwMoGpeMn19fYhEIsydOxdz585FQEAAfv75ZwBVY/IEAgGcnZ3rjK814xny0Gd2H1zZdAVXv7wKx1GONe576/tbKC8sR7vu7eAw3KEJoySEqBPqEVMjHA4H2vrazf5L0fFh4oSIy5We9LK0tBSff/45AEBTU/p/CGdnZ2hoaMg9GP5N169fB1A1876Yt7c3AOD06dMy+4vLxPsockx1auoRW7t2LaKiohAREYHExERERkaiZ8+eAAChUAgej4fw8HA8fPgQW7ZswcaNGyVPQWpra+Ps2bPgcDhITk5GdnY2wsLCJHUnJSVBJBKhR48edcbX2nks8ICGpgbSLqUh80ZmtftUllfi+taq9xUtZ0QIqZXSZyYjNaIJXeUjFAqZjo4O09XVZU+ePJGUz5kzhwFgANjixYtljuvZsyczNjZmIpFIZlt8fDx79eqVTPlff/3FdHR0GI/HY2lpaVIx2NvbMx6Px+7cuSMpf3Ny1sTERJm463vMv5WXlzMej8eMjIxkzsPT05Nt2LCh1uPf1L9/f7Z9+3bJz/v372d2dnbV7hsREcEAsF27dtVZL71HGftt6m+MDz77JeCXarffjbjL+OCzL62+ZBVlFU0cHSFE2WhCV9KqaGpqIjAwECUlJfDy8sKnn36KIUOGYOfOnVizZg309PSwe/durFy5EkVFRZLjRo8eDYFAgNjYWJk6Dx06hPbt22PUqFGYP38+lixZAj8/PwwaNAhCoRDffvut1FJCmpqaCAsLg0gkgpeXF2bNmoUlS5bA1dUV8fHx4PP56Nq1q0zc9T3m3+Lj41FWVobevXvL9KKMHDkSK1euhL+/PyIjI1FYWCjZlpqaitmzZ6Nbt24wMTGBgYEBrl27BisrK8k+cXFxcHV1rbbdM2fOgMvlytzyJdUTPwGZ8GsCXj1+JbWNMSaZwNVjgQe42rScESGkFkpP7UiNqEdMfoWFhWzevHnM0tKSaWlpMWtra7Z161bGGGN79+5l5ubmzNDQUKrXKCMjg3G5XDZ//nyZ+qKjo1lAQADr3LkzMzQ0lNT5wQcfsOvXr9cYx/Xr15mfnx8zNjZmurq6zM3NjUVGRtYae0OOEQsLC2MA2KJFi6rdHh8fz9atW8fs7e2ZjY0Ne/36NXv+/DkzMzNjU6dOZWfPnmUJCQns/PnzDABLSUmRHOvn58dWr14tU2dRUREzMDBgo0ePlitGeo9W2Td8H+ODz6IWREmVJ59OZnzwWah+KCvOK1ZRdIQQZWrMHjEOY/97DI00OvFTcwKBAEZG1U/sWFpaisePH8POzg46OjpNHKH6mzRpEk6fPo20tDTo6+urOpxGk5WVhQ4dOuDBgwe4du0aVq1ahaysLMn25cuXY+fOnXj9+rWkZ83a2hpff/01xo0bJ1XXnj17MGPGDFy8eBGDBg2qs216j1ZJOZOCyOGR0NLXwsL0hdBtU/XkbqRvJFJOp8DjUw/4fe2n4igJIcogz+d3Q9GtSdKihIaGorCwEN99952qQ1GqjRs3Yv/+/UhMTERCQgJWrVoFJycnODk5wczMDC9fvsTJkyeRmJiI4OBg7Nq1Cy4uLlK3N4VCIe7evYvs7GwUFBQAqHoSc/369fD395crCSP/sB9qDwsXCwiLhLi5q2rZo+dxz5FyOgUcDQ76fdavjhoIIYQSMdLC2NnZISIiosX1hpWWliIkJASurq7w8fFBcXExTp48CS6XC39/f0yZMgXjx4/H22+/DS6Xi6FDh0qeqBQLDQ3F7t270b59e+zcuRNA1YoEkydPxubNm1VwVupNvOwRANzYdgMVZRWS5Yy6jesGk04mKoyOEKIu6NZkE6Jbk6QloPfoPyrLK7HVfisKMgvgs9YHF9dchEgowswbM9GhbwdVh0cIURK6NUkIIc0QV5sLj089AAAXVl2ASCiC7SBbSsIIIXKjRIwQQhTQZ1YfaBv+swYqLWdECKkPSsQIIUQBOsY66DOrahUEM0czdB1Z+1xxhBDyJlprkhBCFDRo1SCIKkXoMakHOBq0nBEhRH6UiBFCiIJ0jHXgt4XmDCOE1B/dmiSEEEIIURFKxJopmlWENFf03iSEEOWhRKyZ4XKrFggWCoUqjoSQ6onfm+L3KiGEkIajRKyZ0dLSAo/Hg0AgoJ4H0uwwxiAQCMDj8aClpaXqcAghRO3RYP1myNzcHJmZmcjIyICxsTG0tLSk1gwkpKkxxiAUCiEQCFBYWIgOHWjCUkIIUQZKxJoh8fIJubm5yMzMVHE0hPyDx+OhQ4cOSl/igxBCWitKxJopIyMjGBkZQSgUorKyUtXhEAIul0u3IwkhRMkoEWvmtLS06MOPEEIIaaFosD4hhBBCiIpQIkYIIYQQoiKUiBFCCCGEqAglYoQQQgghKkKJGCGEEEKIitBTk01IPFN+fn6+iiMhhBBCiLzEn9uNseINJWJNqKCgAABgY2Oj4kgIIYQQUl8FBQUwNjZWap0cRgsaNhmRSISsrCwYGhoqbcmi/Px82NjY4OnTpzTbuZqia9gy0HVUf3QN1V9jXUPGGAoKCtC+fXtoaCh3VBf1iDUhDQ0NWFtbN0rd4pn4ifqia9gy0HVUf3QN1V9jXENl94SJ0WB9QgghhBAVoUSMEEIIIURFKBFTczweD8HBweDxeKoOhTQQXcOWga6j+qNrqP7U8RrSYH1CCCGEEBWhHjFCCCGEEBWhRIwQQgghREUoESOEEEIIURFKxAghhBBCVIQSsSbw7NkzzJw5E1ZWVtDR0UHXrl2xZs0alJeX17uuU6dOYfDgwTAyMoKhoSEGDx6MU6dOKaXtFy9eYMOGDRg3bhzs7OzA4XDkXgHgt99+w7Bhw2BmZgZdXV3Y2dlh4sSJePr0ab3PsTlSl2sIVK3g8O2338LFxQW6urpo27YtAgICkJSUVO3+jDEcOXIEPj4+sLKygp6eHhwdHTF79mykpqbW+/xULTY2FiNGjICpqSn09fXh7u6OAwcO1KuO+r6GDWk3Pz8fixYtgq2tLXg8HmxtbbFo0aJa16I9cOAA3N3doa+vD1NTU4wYMQI3b96s17mpg5Z6DUtKSrB582b07t0bpqamMDExgaurK0JDQyEQCOp1fs2dOlzDu3fv4osvvoCvry/atm0LDoeDwYMH1xlXeXk5Nm/eDDc3NxgaGsLQ0BDdu3fHvHnz6nV+Eow0quzsbNaxY0fG4XDY+++/z5YvX84GDBjAADA/Pz9WWVkpd12RkZEMADM3N2effPIJmz9/PrOwsGAAWGRkpMJtX7hwgQFgHA6Hde3alenp6bG63iIikYjNmjWLAWAODg5s7ty5bPny5WzKlCmsY8eO7K+//pL7/JordbqGjDH20UcfMQCsW7dubOnSpWzq1KmMx+MxY2NjFh8fL7P/okWLGABmZWXF5syZw5YtW8Z8fX0Zh8NhhoaG7P79+/V7wVTowoULTFtbmxkYGLCZM2eyxYsXMzs7OwaAhYaGyl1PfV/D+rZbWFjIevbsyQCwYcOGseXLlzM/Pz8GgPXs2ZMVFhbKHBMaGsoAsI4dO7JFixaxWbNmMSMjI6atrc0uXLhQr9epOWup17C8vJx5eHhItn/66afss88+Y66urgwAc3Z2ZkVFRfV/wZohdbmGwcHBDADT1tZm3bt3ZwCYt7d3rTHl5eUxd3d3BoD179+fLV68mC1evJiNGTOGmZmZyX1ub6JErJFNnTqVAWDbt2+XlIlEIjZt2jQGgO3Zs0euevLy8piJiQkzNzdn6enpkvKsrCxmaWnJTExMWF5enkJtP3v2jF28eJHl5+czxhhzdHSsMxHbunUrA8DmzZvHKioqZLYLhUK5zq85U6dreP78eQaAeXl5sdLSUkn52bNnGYfDYYMGDZLaPzs7m2loaLBOnToxgUAgtW3Lli0MAPvwww/lOj9VEwqFzMHBgfF4PHb79m1JeX5+PnN2dmaamprs0aNHddZT39ewIe2uXr2aAWDLli2rtnz16tVS5Y8ePWKampqsa9eu7PXr15LyBw8eMD09Pebg4NAiftda8jX8+eefGQA2ZswYmXhHjx7NALCIiIg6z625U6dr+ODBA3br1i1WXl7OsrOz5UrE3n//fcbhcNj+/furPfeGoESsEeXn5zMej8fs7e2ZSCSS2paVlcU0NDSYp6enXHXt2rWLAWAhISEy2/773/8yAGzXrl1KbbuuRKy4uJi1adOG2dvbt4gPgeqo2zWcOHEiA8AuXrwo04b4v/XExERJ2dWrVxkAFhgYKLP/o0ePGAA2cuRIuc5P1U6dOlVj4njw4EEGgK1YsaLOeur7Gta3XZFIxNq3b88MDAxkek1KSkqYqakp69Chg9Q1X7FiRY0f1HPmzGEA2KlTp+o8t+auJV/DDRs2MADshx9+kGnj+++/ZwDY//3f/9V5bs2dulzDf5MnEbt27RoDwKZMmVJn/PVBY8Qa0dWrV1FWVoZhw4bJjLWysrJCjx49cP36dZSWltZZV3R0NABg+PDhMtt8fX0BABcvXmyUtmty5swZ5OXlYfTo0aisrMSRI0fw3//+Fzt37kRycnKD621O1O0aRkdHQ19fHwMGDJCrjS5dukBbWxsxMTEoKCiQ2j8qKgoAMGTIkDrPrTmo7fUVl7157rXVU5/XsL7tJiUlISsrCwMGDIC+vr7U/jo6Ohg0aBAyMzOlfofq+95RVy35Gjo7OwMATp48KdPGiRMn5B6f1NypyzVsiJ9//hkAMH78eOTm5mLPnj3YsGEDIiMj8fLlywbXq6lQVKRW4gGFXbp0qXZ7ly5dcO/ePaSmpqJbt24Nrktc9uYARmW2XRPxIGFNTU24uroiMTFRsk1DQwMLFy7El19+2aC6mwt1uoZFRUXIzs5G9+7dweVy5WrDzMwMoaGhWLp0Kd566y34+/vD0NAQ9+/fx9mzZzFr1izMnz+/1vNqLmp7vUxNTWFubl7rIF8ADXoN69uuPNdVvN+b3xsYGMDS0lKumNRVS76G7777LkaNGoXDhw+jT58+8Pb2BlCVQCQnJ2P79u1wc3Or9dzUgbpcw4YQf+YlJydjypQpUg9YGBgYICwsDBMmTKh3vdQj1ojEF8nY2Lja7UZGRlL7NbQufX19cLlcqXqU2XZNXrx4AQD46quvYGRkhBs3bqCgoACXLl1C165d8dVXX2HHjh0Nrr85UKdr2NBYlyxZgv3790MgEGDHjh3YtGkTTpw4gb59+2Ly5MnQ0tKq89yaA3nOv67r1JDXsL7tNrSNxvxdbi5a8jXkcDj47bffsGTJEty5cwdbtmzBli1bcOfOHYwePRp+fn61npe6UJdr2BDiz7ylS5fivffeQ0pKCl69eoXIyEhoaGhgypQpiIuLq3e9lIjJwdzcXDKVgzxf4i7Slk4kEgEAtLW1cfToUfTt2xcGBgbw8vLCr7/+Cg0NDXz11VcqjrIKXcOarVu3DkFBQVixYgWePn2KwsJCXL58GRUVFfDx8cGRI0dUHSIhaq+kpARjxozBvn37cODAAeTm5uLly5c4dOgQzpw5g759+yIlJUXVYZJaiD/zXFxcEB4eDnt7e5iYmCAwMBAbN26EUCjEtm3b6l0v3ZqUw8SJE2XGz9RGfPtAnJnXlIWL55qpKYN/05t1mZmZSW0rKipCZWWlVD3KbLuumNzc3NC+fXupbc7OzrC3t0dycjJev34NExOTBrejDK3hGjYk1vPnz2PVqlVYuHAhvvjiC0n5gAED8Mcff8De3h4LFy7EmDFj6jw/VZPn/Ou6Tg15DevbbkPbaMzf5eaiJV/DDRs24Pjx4zh27Bj8/f0l5ePHj4ehoSHeeecdrFmzBhEREbWeX3OnLtewIcTHv/vuuzLjdkeNGoWPP/64QfP6USImh2+++aZBx9U1diMpKQkaGhqwt7eXq66bN28iKSlJ5kO8unvjymy7Jo6OjgBQY5IlLi8pKVF5ItYarqG+vj6srKzw+PFjVFZWyoytqK6NP//8EwDg4+MjU3/btm3Ro0cPXL16Fbm5uTA3N6/zHFXpzderT58+UttevXqF3Nxc9O/fv9Y6GvIa1rddea5rdW1cvXoVz549kxknVtd4JXXSkq9hbb9rPj4+4HA4uHXrVq3npg7U5Ro2hKOjI27evFnt59mbn3f1RbcmG1G/fv3A4/Fw5swZMMaktmVnZ+P+/fvw8PCAjo5OnXWJB3aePn1aZpt4VnbxPspuuybiPygPHz6U2SYUCpGcnAx9fX20bdu2wW2omrpdQ29vbxQVFSEmJkauNsSz8+fk5FQbs7icx+PVeX6qVtvrKy5789xrq6c+r2F92+3SpQvat2+PmJgYFBUVSe1fWlqKS5cuoX379ujcubNcbVQXk7pqydewtt+13NxcMMbU4vesLupyDRtC/AR5QkKCzDZxWadOnepfsVInwyAy6jshZ1FREXv48CFLS0uTKs/Ly2PGxsaNOhnov8kzoevw4cOrnRtnzZo1DACbPHlyrcerA3W6hm9OglhWViYpr2kSxJ9++onhf7N6vzlRKGOMhYeHMwCsT58+8rxMKicUCpm9vT3j8Xjszp07kvI3J3R8c96hnJwc9vDhQ5aTkyNVT31fw/q2y1j9JwNNTExsNRO6ttRrOHv2bAaATZ06VWry68rKSjZ9+nQGgC1evFi+F6oZU6dr+CZ55hETCATM3Nyc6ejosLi4OEl5WVkZe+eddxgAFhYWVuPxNaFErJFlZWUxGxsbxuFw2JgxY9jnn38uWaLG19e3xmWGqnsz7Nu3j+GN5XEWLFggWR5n3759CrfNGGPTpk2TfBkZGTEAUmX//mVJTk5m7dq1Y/jfxJ+LFy9mQ4YMYQCYra0ty87OVuwFbAbU7RrOnDmTQc5lQSoqKtjgwYMZANa2bVs2Y8YMtmTJEjZs2DAGgPF4PLVapur8+fNMS0uLGRgYsI8++khqiZN169ZJ7Ste3iQ4OFimnvq8hvVtlzHZ5XE+//xzyR/ympY4WrduHcMbSxzNnj2bGRkZMS0tLXb+/PmGv2jNTEu9hunp6czKykryj8/8+fPZggULWI8ePRgA1qlTJ/bixQvFXrxmQl2u4cOHDyWfbQEBAQwAs7CwkJRVlxj/9ttvjMvlMj09PTZ16lT26aefMmdnZwaAjRgxotoVZupCiVgTyMrKYtOnT2cWFhZMW1ubde7cmYWEhEgt2yBW24c4Y4ydOHGCDRo0iBkYGDADAwM2aNAgdvLkSaW0zRhjAGr9evz4scwx6enpLCgoiFlaWjItLS1mY2PD5s2bx54/fy7X66MO1OkaVlZWsm3btjFnZ2fG4/GYmZkZGzduXI3/DZaWlrKNGzey3r17Mz09Paapqck6dOjAJk2apFbrTIpdv36d+fn5MWNjY6arq8vc3NyqXceztg+A+r6G9WlX7PXr12zhwoXMxsZG8nuzcOFCmZ7JN0VGRjI3Nzemq6vLjI2NmZ+fH7tx40btL4gaaqnXMDs7m82fP5917tyZaWtrMx6Px7p27coWLVrEcnNz635h1Ig6XEPx3+qavmxtbas97vLly8zPz4+ZmJgwbW1t5uzszDZu3NjgXmkOY/8afEIIIYQQQpoEDdYnhBBCCFERSsQIIYQQQlSEEjFCCCGEEBWhRIwQQgghREUoESOEEEIIURFKxAghhBBCVIQSMUIIIYQQFaFEjBBCCCFERSgRI4QQQghREUrECCGEEEJUhBIxQkir98svv+Cdd95Bhw4dwOPxoKenBxcXF3z77bcy+167dg0cDge7d+9WQaSEkJZGU9UBEEKIKn355ZdYunQp3nrrLYwdOxbGxsYoKipCfHw8YmJi8Mknn0jtf+zYMWhoaODdd99VUcSEkJaEFv0mhLRqlpaWEIlEyMrKgqam9P+mZWVl4PF4UmXdunWDqakpYmJimjJMQkgLRbcmCSGtmrm5OXJycjBp0iT8+uuvyMnJkWz7dxKWnJyMhw8f4r333pOpJzg4GBwOB99//32NbcXExIDD4SAgIEB5J0AIUWuUiBFCWrWvv/4aLi4uGDx4MC5evIg+ffpg8ODBuH//vsy+R48eBQCZRCw5ORn//e9/0atXL8ycObPGtnr16gUAuHLlivJOgBCi1ujWJCGk1YqIiMB3332H8+fPw8DAAACQl5eHIUOGIDExEXfv3oWjo6Nkfy8vL7x48QKJiYlS9UyaNAk//fQTjh8/jlGjRtXaZrt27ZCTk4OCggJJm4SQ1ot6xAghrdKFCxcwffp0fPvtt1IJUZs2bbBs2TKUlpbixx9/lJTn5ubiypUrGD16tFQ9T58+xaFDh2Bvby/XAH49PT0AwOvXr5VyHoQQ9UaJGCGk1WGM4ZNPPoGjoyPc3d1ltpuZmQEAXrx4ISn7/fffIRKJZG5L/vzzz6isrERAQAA4HI6k/ODBg5g8eTIKCgqk9i8pKQEAaGtrS8q2bNkCGxsb6OrqYsiQIXj06JHiJ0kIUQuUiBFCWp07d+4gISEBQ4YMqXZ7amoqAMDKykpSduzYMbRr1w79+vWT2vfSpUsAqm5bvik8PBy///47DA0NJWX5+fl48eIFjIyM0LZtWwDAgQMH8MUXX2Djxo2IjY2Fqakp/Pz8UFZWpviJEkKaPUrECCGtTlxcHICqqSuq8/vvvwMABgwYAKCqF+vMmTMYNWoUNDSk/2wmJSUBgNRYstLSUly+fBk2NjZS+/71118AgIEDB0p6z7Zs2YK5c+di0qRJ6N69O8LDw5GdnY1jx44pepqEEDVAiRghpNUR3x6sbpzW3bt3cebMGVhbW0t6zM6cOYPi4uJqp60QCAQAAH19fUnZ0aNHUVRUJDP9xcGDBwEA48ePBwCUl5fjzp07Uj1zhoaG8PDwwLVr1xQ4Q0KIuqBEjBDS6rz11lsAgN9++w3FxcWS8ufPnyMwMBAVFRXYvHkztLS0AFTdltTT08PQoUNl6hKPJ3vw4AEAoLi4GCEhIXB1dcXjx48hFAoBANevX8dPP/0EGxsbTJw4EUDVAwCVlZVo166dVJ3t2rXD8+fPlXzWhJDmiJY4IoS0Ot7e3hgwYABiYmLQu3dvvPfee3j9+jV++eUXCAQCbN68WdJrJRKJ8Mcff2D48OHQ1dWVqeu9997DgwcP8OGHH2LSpEk4deoUysvLsWnTJvj7+8Pf3x8ODg6IiIiApqYmfvzxR5meMkJI60U9YoSQVofD4eCPP/7AvHnzUFBQgM2bN+PYsWMYOnQorl69ioULF0r2vXr1Kl68eFHtbUkA+M9//oNZs2ahqKgIO3bsgI2NjWQ82YIFC3Djxg1ERERgwIABuHz5MgYPHiw51tzcHFwuV+rpTKDqaU0LC4tGOXdCSPNCE7oSQkgtli1bhq+++grPnz+Hubm50uvv27cvvL298eWXXwIACgsL0bZtW0RERNBSSIS0ApSIEUJILRwdHdGuXTvJE4/Ktn//fnz00UfYs2cPunfvjpCQENy6dQsJCQnQ0dFplDYJIc0HjREjhJBa/Hs5I2ULDAzEixcvsGTJEuTm5sLT0xMnTpygJIyQVoJ6xAghhBBCVIQG6xNCCCGEqAglYoQQQgghKkKJGCGEEEKIilAiRgghhBCiIpSIEUIIIYSoCCVihBBCCCEqQokYIYQQQoiKUCJGCCGEEKIilIgRQgghhKgIJWKEEEIIISpCiRghhBBCiIpQIkYIIYQQoiL/D/NykyjN0/67AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "frequency_window = np.loadtxt('frequency_window_21pt.txt')\n",
    "\n",
    "alpha_gas_0p5Isat = np.loadtxt('alpha_gas_0p5Isat.txt')\n",
    "alpha_gas_1Isat = np.loadtxt('alpha_gas_1Isat.txt')\n",
    "alpha_gas_50Isat = np.loadtxt('alpha_gas_50Isat.txt')\n",
    "alpha_gas_100Isat = np.loadtxt('alpha_gas_100Isat.txt')\n",
    "alpha_gas_500Isat = np.loadtxt('alpha_gas_500Isat.txt')\n",
    "\n",
    "coeff_ad_hoc = 158\n",
    "\n",
    "##Plot\n",
    "\n",
    "font_size = 14\n",
    "font_size_lgd = 14\n",
    "tick_font_size = 14\n",
    "#PLot population\n",
    "plt.plot((frequency_window[0:len(frequency_window)]-omega_0)/omega_0, coeff_ad_hoc*alpha_gas_0p5Isat[0:len(alpha_gas_0p5Isat)], \"-\", c=\"r\", label=\"$\\\\alpha (0.5 ~I_{sat})$\")\n",
    "#plt.scatter((frequency_window[0:len(frequency_window)]-omega_0)/omega_0, coeff_ad_hoc*alpha_gas_0p5Isat[0:len(alpha_gas_0p5Isat)], marker=\"x\",c=\"lightcoral\", label=\"$\\\\alpha (0,5 ~I_{sat})$\")\n",
    "plt.plot((frequency_window[0:len(frequency_window)]-omega_0)/omega_0, coeff_ad_hoc*alpha_gas_1Isat[0:len(alpha_gas_1Isat)], \":\", c=\"k\", label=\"$\\\\alpha (1 ~I_{sat})$\")\n",
    "#plt.scatter((frequency_window[0:len(frequency_window)]-omega_0)/omega_0, coeff_ad_hoc*alpha_gas_1Isat[0:len(alpha_gas_1Isat)], marker=\"x\",c=\"gold\", label=\"$\\\\alpha (1 ~I_{sat})$\")\n",
    "plt.plot((frequency_window[0:len(frequency_window)]-omega_0)/omega_0, coeff_ad_hoc*alpha_gas_50Isat[0:len(alpha_gas_50Isat)]+0.025, \"-\", c=\"g\", label=\"$\\\\alpha (50 ~I_{sat})$\")\n",
    "#plt.scatter((frequency_window[0:len(frequency_window)]-omega_0)/omega_0, coeff_ad_hoc*alpha_gas_0p5Isat[0:len(alpha_gas_50Isat)], marker=\"x\",c=\"limegreen\", label=\"$\\\\alpha (50 ~I_{sat})$\")\n",
    "plt.plot((frequency_window[0:len(frequency_window)]-omega_0)/omega_0, coeff_ad_hoc*alpha_gas_100Isat[0:len(alpha_gas_100Isat)]+0.05, \"-\", c=\"b\", label=\"$\\\\alpha (100 ~I_{sat})$\")\n",
    "#plt.scatter((frequency_window[0:len(frequency_window)]-omega_0)/omega_0, coeff_ad_hoc*alpha_gas_0p5Isat[0:len(alpha_gas_100Isat)], marker=\"x\",c=\"cyan\", label=\"$\\\\alpha (100 ~I_{sat})$\")\n",
    "plt.plot((frequency_window[0:len(frequency_window)]-omega_0)/omega_0, coeff_ad_hoc*alpha_gas_500Isat[0:len(alpha_gas_500Isat)]+0.24, \"-\", c=\"purple\", label=\"$\\\\alpha (500 ~I_{sat})$\")\n",
    "#plt.scatter((frequency_window[0:len(frequency_window)]-omega_0)/omega_0, coeff_ad_hoc*alpha_gas_0p5Isat[0:len(alpha_gas_500Isat)], marker=\"x\",c=\"magenta\", label=\"$\\\\alpha (500 ~I_{sat})$\")\n",
    "plt.xlabel(\"$\\delta / \\omega_0$\", fontsize=font_size), plt.ylabel(\"Absorption coefficient $\\\\alpha ~[1/a_{0}]$\", fontsize=font_size)\n",
    "plt.xticks(fontsize=tick_font_size)  # X-axis graduation font size\n",
    "plt.yticks(fontsize=tick_font_size)  # Y-axis graduation font size\n",
    "# Set the number of ticks on the x-axis to a maximum of 3\n",
    "ax = plt.gca()\n",
    "ax.xaxis.set_major_locator(MaxNLocator(nbins=5))  # Set the maximum number of ticks\n",
    "plt.legend(fontsize=font_size_lgd)\n",
    "plt.savefig('CHAPTER_4_LambDip.png', dpi=300, bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7bd0093",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}