SESS6067 Wind Farm Design Ammended.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Design of wind farm - Ammended by Thomas Le Moignan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Python created on Fri Feb 03 20:53:32 2017\n",
    "\n",
    "@author: J.Wellicome, modified by S.Turnock and Artur K. Lidtke for SESS6058 coursework\n",
    "\n",
    "Computes Power & thrust data for Turbines of specified geometry.\n",
    "Made into Juypter notebook in 2020 by SRT\n",
    "Modified for MTEC TCE1 Module June/July 2020\n",
    "Modified for sess6067 Feb 2021\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.interpolate import interp1d\n",
    "import statistics\n",
    "\n",
    "def readInputData(filename):\n",
    "    \"\"\" Imports data from a file. See the example file for details of the accepted format\n",
    "    \n",
    "    Arguments:\n",
    "        filename - path to the file\n",
    "    \n",
    "    Returns:\n",
    "        data - dictionary of all variables and arrays defined in the data file\n",
    "        variables - labels of the variables (useful for arrays only)\n",
    "    \"\"\"\n",
    "    # get a string representation of the input file\n",
    "    with open(filename) as f:\n",
    "        s = f.readlines()\n",
    "    \n",
    "    # keep each line separately, ignore any comments and empty lines\n",
    "    s = [l.replace('\\n','') for l in s if (l and l[0] != '#' and l!='\\n')]\n",
    "    \n",
    "    data,variables = {},{}\n",
    "    currentData = ''\n",
    "    readingData = False\n",
    "    \n",
    "    # go over each line and extract information\n",
    "    for l in s:\n",
    "        # a new data array\n",
    "        if l.startswith('data') or l.startswith('DATA'):\n",
    "            data[l.split('=')[1].replace('\"','').strip()] = []\n",
    "            currentData = l.split('=')[1].replace('\"','').strip()\n",
    "            readingData = True\n",
    "        \n",
    "        # variables for the current data array\n",
    "        elif l.startswith('variables') or l.startswith('VARIABLES'):\n",
    "            variables[currentData] = [v for v in l.split('=')[1].replace('\"','').split(',')]\n",
    "        \n",
    "        # line of actual data for a data array\n",
    "        elif readingData:\n",
    "            # keep all data as strings, let the other parts of the code worry about conversion\n",
    "            data[currentData].append([v for v in l.split()])\n",
    "            \n",
    "        else:\n",
    "            # individual variable - keep the value as string and let the other parts\n",
    "            # of the code worry about conversion\n",
    "            data[l.split('=')[0].strip()] = l.split('=')[1].strip()\n",
    "            variables[l.split('=')[0].strip()] = l.split('=')[0].strip() # this label is meaningless but keep it for consistency\n",
    "    \n",
    "    # return all the data\n",
    "    return data,variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def GOLDK(N,X,TPHI):\n",
    "    \"\"\" Computes the Goldstein correction factor\n",
    "    \n",
    "    Arguments:\n",
    "        N - number of blades\n",
    "        X - non-dimensional radius\n",
    "        TPHI - tan of the PHI angle\n",
    "    \n",
    "    Returns:\n",
    "        K - Goldstein correction factor\n",
    "    \"\"\"\n",
    "    K = 1.\n",
    "    if np.abs(TPHI) > 0.001:\n",
    "        F = N/2.0/X/TPHI-0.50\n",
    "        if np.abs(F) > 85.:\n",
    "            K = 1.\n",
    "        else:\n",
    "            K = 2.0*np.arccos(np.cosh(X*F)/np.cosh(F))/np.pi\n",
    "    \n",
    "    return K"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Rayleigh(U,Umean):\n",
    "    \"\"\"\n",
    "    calculates the probablity of rayleigh function\n",
    "    \"\"\"\n",
    "    A=np.pi*U/(2*Umean*Umean)\n",
    "    B=-np.pi*np.square(U/(2*Umean))\n",
    "    prob=A*np.exp(B)\n",
    "    return prob"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def interp(ALPHA,A,C):\n",
    "    \"\"\" Performe Lagrange polynomial interpolation\n",
    "    \n",
    "    Arguments\n",
    "        ALPHA - value at which to interpolate\n",
    "        A - reference x data\n",
    "        C - reference y data\n",
    "        \n",
    "    Returns\n",
    "        COEFT - interpolated value of C at ALPHA\n",
    "    \"\"\"\n",
    "    \n",
    "    N = A.shape[0]\n",
    "    M = 4\n",
    "    \n",
    "    for i in range(1,N):\n",
    "        if ALPHA < A[0]:\n",
    "            COEFT = C[0]\n",
    "            break\n",
    "        \n",
    "        elif ALPHA >= A[-1]:\n",
    "            COEFT = C[-1]\n",
    "            break\n",
    "        \n",
    "        elif A[i-1] <= ALPHA and A[i] > ALPHA:\n",
    "            \n",
    "            I = i-M/2\n",
    "            if I < 0:\n",
    "                I = 0\n",
    "            if I > (N-M):\n",
    "                I = N-M\n",
    "\n",
    "            I = int(I)   #added to eliminate integer error\n",
    "            \n",
    "            S,V = np.zeros(M),np.zeros(M)\n",
    "            \n",
    "            for J in range(1,M+1):\n",
    "                V[J-1] = C[J+I-1]\n",
    "                S[J-1] = A[J+I-1]\n",
    "            \n",
    "            for J in range(1,M):\n",
    "                for I in range(J+1,M+1):\n",
    "                    V[I-1] = (V[I-1]-V[J-1])/(S[I-1]-S[J-1])\n",
    "    \n",
    "            U = V[M-1]\n",
    "            for I in range(M-1,0,-1):\n",
    "                U = (ALPHA-S[I-1])*U+V[I-1]\n",
    "            \n",
    "            COEFT = U\n",
    "            break\n",
    "    \n",
    "    return COEFT\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def TurbinePerformance(dB,dA,Bpitch,SF,NB,Dm,Vw,BRPM):\n",
    "    \"\"\" Performe Analysis of Turbine at specifc condition\n",
    "\n",
    "     Arguments\n",
    "        dB Data definig blade shape\n",
    "        dA Data defining section performance\n",
    "        Bpitch = set blade pitch\n",
    "        SF = set scale factor\n",
    "        NB= number of blades\n",
    "        Dm- diameter in m\n",
    "        Vw=speed in m/s\n",
    "        BRPM = rat eof rotation per minute\n",
    "               \n",
    "    Returns\n",
    "        TSR = tip speed ratio\n",
    "        CPOWER = power coeffiencit\n",
    "        CTHRUST thrust coefficient\n",
    "    \"\"\"\n",
    "    \n",
    "    # Set initial relaxation factor (MU) for updating ALPHA in iteration loop\n",
    "    MU = 0.9\n",
    "    # IMCRIT=Count of secs MACH>MCRITA\n",
    "    IMCRIT = 0\n",
    "    A1 = 0.\n",
    "    # Mach number corrections\n",
    "    # print(Bpitch,SF,NB,Dm,Vw,BRPM)\n",
    "    DIFCP = np.zeros(dB.shape[0])\n",
    "    DIFCT = np.zeros(dB.shape[0])\n",
    "    ANGLE = np.zeros(dB.shape[0])\n",
    "    LIFT = np.zeros(dB.shape[0])\n",
    "    DRAG = np.zeros(dB.shape[0])\n",
    "    PHIS = np.zeros(dB.shape[0])\n",
    "    AAXIAL = np.zeros(dB.shape[0])\n",
    "    ATANG = np.zeros(dB.shape[0])\n",
    "    GOLD = np.zeros(dB.shape[0])\n",
    "    MACHNO = np.zeros(dB.shape[0])\n",
    "    MCRITS = np.zeros(dB.shape[0])\n",
    "    DIFFS = np.zeros(dB.shape[0])\n",
    "    \n",
    "    #   set up blade parameters  \n",
    "    for IXCOUNT in range(dB.shape[0]):\n",
    "        \n",
    "        X = dB[IXCOUNT,0]\n",
    "        THETA = dB[IXCOUNT,2]/180.*np.pi + Bpitch\n",
    "        SIGMA = SF*dB[IXCOUNT,1]*Nblades # solidity\n",
    "        TSR = Dm*np.pi*BRPM/(Vw*60.0)\n",
    "        MACH=Vw*np.sqrt(1.0+(TSR*X)**2)/VSOUND   \n",
    "        if MACH < 0.97:\n",
    "            COMPFAC=1.0/(np.sqrt(1.0-MACH*MACH))\n",
    "        else:\n",
    "            COPMFAC=4.0  \n",
    "        # loop over angle of attack\n",
    "        # initial estimate\n",
    "        ALPHA=np.arctan(0.8/(TSR*X))-THETA\n",
    "        for ICOUNT in range(nAlphaItersMax):\n",
    "            PHI=ALPHA+THETA\n",
    "            TPHI=np.tan(PHI)\n",
    "            \n",
    "            # interpolate the force coefficients\n",
    "            CLIFT = COMPFAC*interp(ALPHA/np.pi*180,dA[:,0],dA[:,1])\n",
    "            CDRAG = COMPFAC*interp(ALPHA/np.pi*180,dA[:,0],dA[:,2])\n",
    "            GAMMA = np.arctan(CDRAG/CLIFT)\n",
    "            \n",
    "            # Goldstein correction due to tip loss\n",
    "            K = GOLDK(NB,X,np.tan(TPHI))\n",
    "            \n",
    "            # axial inflow factor\n",
    "            AA=4.0*np.pi*X*K*np.sin(PHI)*TPHI/(1.0+TPHI*CDRAG/CLIFT)\n",
    "            AB=SIGMA*CLIFT\n",
    "            AA=AB/(AB+AA)\n",
    "            \n",
    "            # tangential inflow factor\n",
    "            AT=AA*np.tan(PHI-GAMMA)/TSR/X\n",
    "            \n",
    "            # derivative of power coefficient\n",
    "            DCPDX=2.0*X*TSR*SIGMA*CLIFT*(1.0-np.tan(GAMMA)/TPHI)/np.sin(PHI)/np.pi\n",
    "            AB=1.0-AA\n",
    "            DCPDX=DCPDX*AB*AB\n",
    "            \n",
    "            # calculate PHI from the velocity vector and use it to calculate ALPHA\n",
    "            PHI=np.arctan((1.0-AA)/(1.0+AT)/TSR/X)\n",
    "            A0=PHI-THETA\n",
    "            \n",
    "            DIFF=np.abs(A0-ALPHA)\n",
    "            A2=A1\n",
    "            A1=ALPHA\n",
    "            ALPHA=(1.0-MU)*ALPHA+MU*A0\n",
    "            \n",
    "            if ((ICOUNT == 5)  or (ICOUNT == 10)):\n",
    "                # Apply Aitken Extrapolator to esimate a likely final ALPHA value as an aid\n",
    "                # to selecting a suitable relaxation factor. See notes. \n",
    "                # MU re-estimated only at 5th and 10th iterations\n",
    "                A3=(ALPHA*A2-A1*A1)/(ALPHA-2.0*A1+A2)\n",
    "                # Reset relaxation factor, limiting the allowable range of values.\n",
    "                AMU=(A3-A1)/(A0-A1)\n",
    "                if (AMU > 2.0):\n",
    "                    AMU = 2.0   \n",
    "                elif (AMU < 0.4):\n",
    "                    AMU=0.4\n",
    "                MU=AMU\n",
    "                # Reset next ALPHA value\n",
    "                ALPHA=A3\n",
    "            ### end of if  \n",
    "             \n",
    "            if DIFF < 0.00001:\n",
    "                break\n",
    "    \n",
    "        # DCTDX=Derivative of Thrust coefficient from momentum    \n",
    "        DCTDX=8.0*X*K*AA*(1.0-AA)\n",
    "        \n",
    "        MCRITA = interp(ALPHA/np.pi*180,dataAero[:,0],dataAero[:,3])\n",
    "        if (MACH > MCRITA):\n",
    "            IMCRIT=IMCRIT+1\n",
    "        \n",
    "        # store the data for this section at this condition\n",
    "        DIFCP[IXCOUNT]=DCPDX\n",
    "        DIFCT[IXCOUNT]=DCTDX\n",
    "        ANGLE[IXCOUNT]=ALPHA\n",
    "        LIFT[IXCOUNT]=CLIFT\n",
    "        DRAG[IXCOUNT]=CDRAG\n",
    "        \n",
    "        PHIS[IXCOUNT] = PHI\n",
    "        AAXIAL[IXCOUNT] = AA\n",
    "        ATANG[IXCOUNT] = AT\n",
    "        GOLD[IXCOUNT] = K\n",
    "        MACHNO[IXCOUNT] = MACH\n",
    "        MCRITS[IXCOUNT] = MCRITA\n",
    "        \n",
    "        DIFFS[IXCOUNT] = DIFF\n",
    "    \"\"\"\n",
    "    # print a summary of the flow equilibrium found for this condition\n",
    "    outputString += '     X  |  PHI  |  A-AXIAL | A-TANG  |  GOLD  |  MACH   |  MCRIT  |  DIFF\\n'\n",
    "    for IXCOUNT in range(dataBlade.shape[0]):\n",
    "        outputString += '{:7.2f} {:7.2f} {:9.4f} {:9.4f} {:9.4f} {:9.4f} {:9.4f} {:9.6f}\\n'.format(dataBlade[IXCOUNT,0],PHIS[IXCOUNT]/np.pi*180.,\n",
    "            AAXIAL[IXCOUNT],ATANG[IXCOUNT],GOLD[IXCOUNT],MACHNO[IXCOUNT],MCRITS[IXCOUNT],DIFFS[IXCOUNT])\n",
    "        outputString += '\\n'\n",
    "    \n",
    "        outputString += '     X | AoA [deg] |   CL   |   CD   |  DCPDX  |  DCTDX\\n'\n",
    "    for IXCOUNT in range(dataBlade.shape[0]):\n",
    "        outputString += '{:7.2f} {:7.2f} {:10.4f} {:9.4f} {:9.4f} {:9.4f}\\n'.format(dataBlade[IXCOUNT,0],ANGLE[IXCOUNT]/np.pi*180.,\n",
    "            LIFT[IXCOUNT],DRAG[IXCOUNT],DIFCP[IXCOUNT],DIFCT[IXCOUNT])\n",
    "        outputString += '\\n'\n",
    "    \"\"\"\n",
    "    # integrate the coefficients to get turbine performance\n",
    "    CTHRUST = np.trapz(DIFCT,dB[:,0])\n",
    "    CPOWER = np.trapz(DIFCP,dB[:,0])\n",
    "    \n",
    "    return TSR,CTHRUST,CPOWER\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# find pitch angle for given TSR blade etc\n",
    "\n",
    "def find_pitch_angle(dB,dA,Bpitch,Maxpitch,cp_target,SF,NB,Dm,Vw,BRPM):\n",
    "    \"\"\" Pserachto find pitch angle that meets cp_target\n",
    "    simple algorithm that steps in 0.1 deg of pitch\n",
    "         Arguments\n",
    "        dB Data definig blade shape\n",
    "        dA Data defining section performance\n",
    "        Bpitch = initial blade pitch\n",
    "        SF = set scale factor\n",
    "        NB= number of blades\n",
    "        Dm- diameter in m\n",
    "        Vw=speed in m/s\n",
    "        BRPM = rate of rotation per minute\n",
    "               \n",
    "    Returns\n",
    "        found = boolean (1 = found , 0 if not)\n",
    "        new_pitch\n",
    "        \n",
    "    \"\"\"\n",
    "    search_inc = 0.01 # rad\n",
    "    FOUND = False\n",
    "    no_search = (int)((Maxpitch-Bpitch)/search_inc)\n",
    "    \n",
    "    for ic in range(no_search):\n",
    "        apitch = Bpitch + search_inc*ic\n",
    "        tsr,cp,ct = TurbinePerformance(dB,dA,apitch,SF,NB,Dm,Vw,BRPM)\n",
    "        if cp < cp_target:\n",
    "            FOUND = True\n",
    "            new_pitch = apitch*180.0/np.pi\n",
    "            break\n",
    "    return FOUND,new_pitch\n",
    "    \n",
    "    \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\tlemo\\AppData\\Local\\Temp/ipykernel_28384/1256137609.py:9: DeprecationWarning: `np.float` is a deprecated alias for the builtin `float`. To silence this warning, use `float` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.float64` here.\n",
      "Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations\n",
      "  dataAero = np.array(dataIn['aero']).astype(np.float)\n",
      "C:\\Users\\tlemo\\AppData\\Local\\Temp/ipykernel_28384/1256137609.py:10: DeprecationWarning: `np.float` is a deprecated alias for the builtin `float`. To silence this warning, use `float` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.float64` here.\n",
      "Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations\n",
      "  dataBlade = np.array(dataIn['blade']).astype(np.float)\n",
      "C:\\Users\\tlemo\\AppData\\Local\\Temp/ipykernel_28384/1256137609.py:11: DeprecationWarning: `np.float` is a deprecated alias for the builtin `float`. To silence this warning, use `float` by itself. Doing this will not modify any behavior and is safe. If you specifically wanted the numpy scalar type, use `np.float64` here.\n",
      "Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations\n",
      "  dataConditions = np.array(dataIn['conditions']).astype(np.float)\n"
     ]
    }
   ],
   "source": [
    "\n",
    "inputFilename = 'input1.txt'\n",
    "outputFilename = 'output2.txt'\n",
    "outputString = ''\n",
    "\n",
    "# read the input file\n",
    "dataIn,varsIn = readInputData(inputFilename)\n",
    "\n",
    "# convert data types\n",
    "dataAero = np.array(dataIn['aero']).astype(np.float)\n",
    "dataBlade = np.array(dataIn['blade']).astype(np.float)\n",
    "dataConditions = np.array(dataIn['conditions']).astype(np.float)\n",
    "Nblades = int(dataIn['Nblades'])\n",
    "D = float(dataIn['diameter'])\n",
    "SF = float(dataIn['scalefactor'])\n",
    "Uten = float(dataIn['Uten'])\n",
    "HH = float(dataIn['HubHeight'])\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "* item lines above read in default values from input1.txt to alter design you need to uncomment lines below by removing hashtag \n",
    "\n",
    "* D = 120.0  # chosen design diameter of turbine (m)\n",
    "* HH = hub height of wind turbine above sea level, need to add in depth of wtaer to get total tower height above seabed\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7.061761459072345  m/s\n"
     ]
    }
   ],
   "source": [
    "D=362.5 # m\n",
    "HH=206.8 # m\n",
    "clearance = HH - 0.5*D\n",
    "# clearance is minimum height blade is above sea level - obvs tsi need sto be high enough to allow safe pasage of vessels\n",
    "\n",
    "Uten=4.621 # m/s\n",
    "A=0.14* np.log(HH/10.0)\n",
    "vrated = Uten*np.exp(A)\n",
    "print(vrated,' m/s')\n",
    "\n",
    "# this determines the rated design speed for wind turbine\n",
    "# when using for tidal turbine need to remove hash tag and select apropriate max tidal speed for rated\n",
    "\n",
    "# vrated = 2.5 # speed in m/s for tidal assignment\n",
    "no_hours=365*24  # number of hours in a year\n",
    "CF=0.35 # a typical capacity factor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "current tip chord is: 14.1375  m \n",
      " Scale Factor 0.7073386383731212  Gear Box ratio 241.93548387096774\n"
     ]
    }
   ],
   "source": [
    "select_clearance_limit = 25.0\n",
    "if clearance<select_clearance_limit :\n",
    "    outputString+=' Clearance {:6.1f} m not sufficient'.format(clearance)\n",
    "outputString += ' At rated Speed {:6.1f} ms-1 \\n'.format(vrated)\n",
    "\n",
    "VSOUND = 342.0 # m/s only if mach number check switched on to increase drag as mach 1 approached\n",
    "nAlphaItersMax = 20 # max no iterations on alpha loop\n",
    "\n",
    "Density_air=1.223 # kg/m^3\n",
    "Density_salt_water=1025.0 # kg/m^3\n",
    "# choose fluid density\n",
    "Density=Density_air\n",
    "# Density=Density_salt_water\n",
    "\n",
    "\n",
    "\n",
    "# choose operational conditions\n",
    "\n",
    "# D = 182.0 ### change this to a specific value of duiameter in m\n",
    "pitch= -1 #degrees dataConditions[0,1] # set pitch value to first one or enter own value\n",
    "RPM = 3.1  # dataConditions[0,2] # likewise here\n",
    "GBRatio=750/RPM\n",
    "Nblades = 3 # change number of blade\n",
    "ctipchord = D*dataBlade[dataBlade.shape[0]-1,1] \n",
    "print('current tip chord is:',ctipchord,' m ')\n",
    "actual_chord = 10 #  enter student value\n",
    "ScaleFactor = actual_chord/ctipchord\n",
    "print(' Scale Factor', ScaleFactor,' Gear Box ratio',GBRatio)\n",
    "\n",
    "## adjust geometry\n",
    "    \n",
    "\n",
    "   \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x214e1c740d0>]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAzgElEQVR4nO3dd3iUZdbH8e+PKk16h9ARpEsU7Li6q2LBLlZ0UXTXtb/WLfZd26511UURsaFYUGwoVkRECL1LC733Gkhy3j+eJzjGlAEyJcn5XFeuzNxPO5Nkcua5q8wM55xzDqBMogNwzjmXPDwpOOec28uTgnPOub08KTjnnNvLk4Jzzrm9yiU6gANRp04da968eaLDcM65YmXixInrzKxuXtuKdVJo3rw5aWlpiQ7DOeeKFUmL89vm1UfOOef2illSkNRU0jeSZkuaKenGsPxeScslTQm/ekccc5ek+ZLmSjo5VrE555zLWyyrjzKBW81skqRqwERJo8JtT5jZ45E7SzoU6At0ABoBX0pqa2ZZMYzROedchJjdKZjZSjObFD7eCswGGhdwSB/gLTPLMLNFwHzgiFjF55xz7rfi0qYgqTnQDfgpLPqLpGmSXpZUMyxrDCyNOGwZeSQRSQMkpUlKW7t2bSzDds65UifmSUFSVeA94CYz2wI8D7QCugIrgX/n7JrH4b+Zrc/MBppZqpml1q2bZ48q55xz+ymmSUFSeYKE8IaZvQ9gZqvNLMvMsoEX+aWKaBnQNOLwJsCKWMbnnHPu12LZ+0jAIGC2mf0norxhxG5nAzPCxyOAvpIqSmoBtAHGxyo+55wrrl77MZ2x89fF5Nyx7H10NHAZMF3SlLDsbuAiSV0JqobSgWsAzGympGHALIKeS9d5zyPnnPu1kTNW8Y8RMzmjcyOOal2nyM8fs6RgZmPIu53g0wKOeQh4KFYxOedccTZ5yUZuensyXZvW4NHzOsfkGj6i2TnnioGlG3Zw1ZA06lU7iJcuT+Wg8mVjch1PCs45l+Q279jDFYPHk5ltDL7ycGpXrRiza3lScM65JJaRmcXVr6WxdMNOBl7WnVZ1q8b0esV6llTnnCvJsrON/3tnGuMXbeDpi7rRo2XtmF/T7xSccy5JPfr5XD6auoI7T23HmV0axeWanhSccy4JvTZuMS98t4BLe6ZwzXEt43ZdTwrOOZdkPp+5ins+nMFJ7etx7xkdCMYCx4cnBeecSyITF2/ghqGT6dykBk9f1I1yZeP7b9qTgnPOJYn5a7bRf0gajWpUYlC/VCpXiH9fIE8KzjmXBNZs2UW/l8dTrowYcuURMR2LUBDvkuqccwm2Zdce+g2ewMYdu3lrQE9SaldOWCx+p+Cccwm0a08WA15NY97qrbxwaXc6N6mR0Hj8TsE55xIkK9u4ZdgUxi3cwJMXduW4tolfOMzvFJxzLgHMjPs+msmn01fx197tOatbQUvYx48nBeecS4Bnvp7Pqz8u5upjW3B1HAenFcaTgnPOxdlr4xbzn1E/c85hjbnr1PaJDudXYrkcZ1NJ30iaLWmmpBvD8sckzZE0TdJwSTXC8uaSdkqaEn69EKvYnHMuUT6etoJ/fDiDE9vV45FzO1OmTPxGK0cjlncKmcCtZtYe6AlcJ+lQYBTQ0cw6Az8Dd0Ucs8DMuoZf18YwNueci7sx89Zx89tTSG1Wk/9echjl4zxaORoxi8jMVprZpPDxVmA20NjMvjCzzHC3cUCTWMXgnHPJYvKSjQx4LY1WdavyUr/DY7Zy2oGKS5qS1BzoBvyUa9Mfgc8inreQNFnSd5KOzedcAySlSUpbu3ZtbAJ2zrkiNHfVVq4YPIG61Sryav8jqF6pfKJDylfMk4KkqsB7wE1mtiWi/K8EVUxvhEUrgRQz6wbcArwp6eDc5zOzgWaWamapdesmvk+vc84VZMn6HVw26CcOKl+G1/v3oF61gxIdUoFimhQklSdICG+Y2fsR5f2A04FLzMwAzCzDzNaHjycCC4C2sYzPOediafWWXVwyaBy7s7J5rX8PmtZK3PQV0Yo6KUiqIinqSjAFE4APAmab2X8iyk8B7gDONLMdEeV1c84vqSXQBlgY7fWccy6ZbNy+m8sG/cT6bbt55cojaFu/WqJDikq+01xIKgP0BS4BDgcygIqS1gKfAgPNbF4B5z4auAyYLmlKWHY38DRQERgVLhwxLuxpdBxwv6RMIAu41sw2HMBrc865hNi6aw/9Bo8nff0OXrnycLo2rZHokKJW0NxH3wBfEnQZnWFm2QCSagEnAA9LGm5mr+d1sJmNAfLqgPtpPvu/R1DV5JxzxdbO3Vn0H5LGrBVb+N9l3TmqVZ1Eh7RPCkoKJ5nZntyF4af394D3wjYD55xzwO7MbK59fSIT0jfwdN9unNi+fqJD2mf5JoXIhCCpJtA0cn8zm5RX0nDOudIoMyubG4ZO5ruf1/LwOZ04o0ujRIe0XwqdOlvSA8AVBL2BLCw24HexC8s554qPrGzj1nemMnLmKv5++qH0PSIl0SHtt2jWU7gAaGVmu2MdjHPOFTfZ2cZfh0/nwykruO3kQ+h/TItEh3RAoumSOgOoEeM4nHOu2DEz7v94Fm9NWMr1v2vNdSe0TnRIByyaO4V/AZMlzSDolgqAmZ0Zs6iccy7JmRkPfzaHV8amc9UxLbjl9yVjrG00SWEI8AgwHciObTjOOZf8zIzHv5jL/0Yv5LKezfjrae0Jx10Ve9EkhXVm9nTMI3HOuWLi6a/m899vFnDREU2578wOJSYhQHRJYaKkfwEj+HX10aSYReWcc0nqv9/M54kvf+a87k146KxOSbdIzoGKJil0C7/3jCjzLqnOuVJn4OgFPPb5XPp0bZSUq6YVhUKTgpmdEI9AnHMumb04eiH//HQOp3duyL/P70LZEpgQoIAuqZIuDSfFy297K0nHxCYs55xLHi99v5CHPp3NaZ0a8uSFXSmXhMtoFpWC7hRqE3RFnQhMBNYCBwGtgeOBdcCdMY/QOecS6OUxi3jwk9n07tSAJ/uW7IQABc999JSkZwnaDo4GOgM7CdZavszMlsQnROecS4yXxyzi/o9ncUqHBjzVtxvlS3hCgELaFMwsCxgVfjnnXKnx0vcLefCT2ZzSoQHPXFw6EgJE1/vIOedKlRdHB20IvTuVnjuEHDF7pZKaSvpG0mxJMyXdGJbXkjRK0rzwe82IY+6SNF/SXEknxyo255zLz/++WxA0KnduWOoSAsQwKQCZwK1m1p5gjMN1kg4laJz+yszaAF+Fzwm39QU6AKcAz+3LmtDOOXeg/vvNfP712RzO6NKIpy7sWuoSAhS8RvMtBR1oZv8pZPtKYGX4eKuk2UBjoA/QK9xtCPAtcEdY/paZZQCLJM0HjgB+jOaFOOfc/jIznvxyHk99NY+zuzXmsfM6l/heRvkpqE2hWlFdRFJzgpHRPwH1w4SBma2UVC/crTEwLuKwZWGZc87FjJnx6Odzef7bBZzfvQkPn9u5xA5Mi0ZBXVLvK4oLSKpKsKbzTWa2pYCJo/LaYL/ZSRoADABISSm+qxs55xLPzHjwk9kMGrOIS3qk8ECfjiVy6op9Uej9kaS2kr4K11NAUmdJf4vm5JLKEySEN8zs/bB4taSG4faGwJqwfBnBOtA5mgArcp/TzAaaWaqZpdatWzeaMJxz7jeys41/fDiTQWMWccVRzXnwLE8IEF1D84vAXcAeADObRtAgXCAFtwSDgNm52h9GAP3Cx/2ADyPK+0qqKKkF0AYYH82LcM65fZGVbdz+3jReG7eYa45ryT1nHFqipr8+ENGMU6hsZuNz/cAyozjuaOAyYLqkKWHZ3cDDwDBJ/YElwPkAZjZT0jBgVnj+68LBc845V2T2ZGVz89tT+HjaSm4+qS03nNjaE0KEqBbZkdSKsH5f0nmEvYoKYmZjyLudAODEfI55CHgoipicc26f7dqTxV/enMyXs1dzd+92DDiuVaJDSjrRJIXrgIFAO0nLgUXAJTGNyjnnitiO3Zlc89pEvp+3jgf6dOCyI5snOqSkFE1SWGxmJ0mqApQxs62xDso554rS5p17uHLweKYs3cTj53fhvO5NEh1S0oqmoXmRpIEEo5K3xTge55wrUuu2ZdB34DimL9/Mc5cc5gmhENEkhUOALwmqkRZJetYX13HOFQfLN+3kghd+ZNG6bQzqdzindGyY6JCSXqFJwcx2mtkwMzuHYFTywcB3MY/MOecOwIK127jghR9ZuzWD1/v34Li2Pq4pGlFN7iHpeEnPAZMIVl+7IKZROefcAZi+bDPnv/AjGZlZDB3Qk9TmtRIdUrFRaEOzpEXAFGAYcJuZbY91UM45t79+XLCeq19No3ql8rx+VQ9a1KmS6JCKlWh6H3Uxsy0xj8Q55w7Q5zNXcf3QyTSrVZnX+vegQfWDEh1SsRNN9VGD/Z37yDnn4mVY2lL+9PpEDm14MO9ce6QnhP0Us7mPnHMuHsyM/34zn9vfncbRrevwxlU9qFG5QqLDKrZiOfeRc87FVHa2cf/Hs3hlbDpndW3Eo+d1oUK50rk4TlGJ2dxHzjkXS7szs7n1nal8NHUF/Y9pwV97t/epr4vA/s59dGlMo3LOuQJs3bWHP70+iTHz13Hnqe245riWPtNpESk0KZjZQsDnPnLOJYXVW3ZxxeAJzFu91ecxioF8k4KkW/IpByDXwjnOORdz89dso9/L49m4YzeDrjic432UcpEr6E6hWtyicM65QqSlb+CqV9MoV6YMbw84kk5Nqic6pBIp36RgZvfFMxDnnMvPyBkrufGtKTSqUYkhVx5BSu3KiQ6pxIpZ3y1JL0takzPoLSx7W9KU8Cs9Z5lOSc0l7YzY9kKs4nLOFS+DxiziT29M4tBGB/Pen47yhBBj0fQ+2l+vAM8Cr+YUmNmFOY8l/RvYHLH/AjPrGsN4nHPFSFa28UA4BuGUDg14sm9XDipfNtFhlXgxSwpmNlpS87y2KWitvgD4Xayu75wrvnbuzuLGtybzxazV9D+mBXf3bk9ZH4MQF1ElBUmnAR0Ips0GwMzuP4DrHgusNrN5EWUtJE0GtgB/M7Pv84llADAAICUl5QBCcM4lo7VbM7jq1TSmLdvEPWccypVHt0h0SKVKNFNnvwBUBk4AXgLOA8Yf4HUvAoZGPF8JpJjZekndgQ8kdchrdlYzG0gwmI7U1FQ7wDicc0nk59VbuXLwBNZvz+CFS7tzcocGiQ6p1ImmofkoM7sc2Bj2SDoSaLq/F5RUDjgHeDunzMwyzGx9+HgisABou7/XcM4VP9/PW8u5z41ld1Y2w6450hNCgkSTFHaG33dIakQwW+qB3M+dBMwxs2U5BZLqSiobPm4JtAEWHsA1nHPFyJs/LeGKwRNoXLMSH1x3NJ2b1Eh0SKVWNG0KH0uqATxGsBynEVQjFUjSUKAXUEfSMuAeMxtEMO320Fy7HwfcLykTyAKuNbMN0b4I51zxlJVtPDJyDgNHL+T4tnV59uJuVDuofKLDKtVkFn21vKSKwEFmtrnQneMgNTXV0tLSEh2Gc24/bMvI5Ka3JvPl7DVcfmQz/nH6oZQr69Nex4OkiWaWmte2aBqaKwO3EjQEXy0pRdKxZvZxUQfqnCsdlm3cwVVD0pi3Zhv39+nA5Uc2T3RILhRN9dFgYCJBAzPAMuAdwJOCc26fTVy8kWteSyMjM5vBVxzOcT6pXVKJ5l6tlZk9yi/Lce4EfBSJc26fDZ+8jIteHEeViuUY/uejPCEkoWjuFHZLqsQvK6+1AjJiGpVzrkTJyjYeHTmH/41eSI8WtXjh0u7UrOLrKCejaJLCPcBIoKmkN4CjgStiGZRzruTYsmsPNw6dzDdz13JpzxTuOaMD5b1BOWlFs/LaKEmTgJ4E1UY3mtm6mEfmnCv20tdt56pX00hft50HzurIZT2bJTokV4hC07Wko4FdZvYJUAO4W5L/Zp1zBfp+3lr6/PcH1m/L4LX+PTwhFBPR3MM9TzCauQtwG7CYiOmwnXMukpkxcPQC+r08ngYHH8SH1x3Dka1qJzosF6Vo2hQyzcwk9QGeNrNBkvrFOjDnXPGzc3cWd7w3jRFTV9C7UwMeO68LVSrGctkWV9Si+W1tlXQXcClwXDhHkY9Dd879yrKNOxjw6kRmr9rCbScfwp97tSJYOsUVJ9EkhQuBi4H+ZrZKUgrBPEjOOQfA2Pnr+MvQyezJyublfodzQrt6iQ7J7adoeh+tAv4T8XwJ3qbgnCOn/WAhj4ycQ8u6VRl4WXda1q2a6LDcAfDKPufcftmWkcnt707l0+mr6N2pAY+e14Wq3n5Q7Plv0Dm3zxas3cY1r01k4dpt3N27HVcf29LbD0qIaNdorkQwS+rcGMfjnEtyI2es4v/emUqFcmV4vX8PjmpdJ9EhuSIUzeC1M4ApBFNdIKmrpBExjss5l2T2ZGXz0CezuPb1ibSqV5WPrz/GE0IJFM3gtXuBI4BNAGY2BWhe2EGSXpa0RtKMiLJ7JS2XNCX86h2x7S5J8yXNlXTyvr0M51wsrd6yi4tfHMeL3y/i8iObMeyanjSqUSnRYbkYiHbw2ub9qC98BXiW3/ZUesLMHo8skHQowTKdHYBGwJeS2ppZ1r5e1DlXtMYuWMcNQyezPSOLp/p2pU/XxokOycVQNElhhqSLgbKS2gA3AGMLO8jMRktqHmUcfYC3zCwDWCRpPsHdyY9RHu+cK2LZ2cbz3y3g31/MpUWdKgy9uidt6ldLdFguxqKpPrqe4BN8BvAmsBm46QCu+RdJ08LqpZphWWNgacQ+y8Ky35A0QFKapLS1a9ceQBjOufxs2L6bK1+ZwGOfz6V3p4Z8+JdjPCGUEtEkhUPM7K9mdnj49Tcz27Wf13seaAV0BVYC/w7L86qbsrxOYGYDzSzVzFLr1vVVm5wramnpG+j91Pf8uHA9D57VkWcu6ubjD0qRaH7T/5HUkGBd5rfMbOb+XszMVuc8lvQiv6zzvAxoGrFrE2DF/l7HObfvsrONgd8v5LHP59KkZiXe/9NRdGxcPdFhuTgr9E7BzE4AegFrgYGSpkv62/5cLEwuOc4GcnomjQD6SqooqQXQBhi/P9dwzu27Ddt303/IBB7+bA6ndGjAx9cf4wmhlIrqnjCc/+hpSd8AtwP/AB4s6BhJQwmSSR1JywiW9ewlqStB1VA6cE14/pmShgGzgEzgOu955Fx8/LhgPTe9PZmN2/dwf58OXNazmY9OLsVklmfV/S87SO0JZko9H1gHvAW8Z2ZrYh9ewVJTUy0tLS3RYThXLGVmZfP01/N55ut5tKhdhWcu7kaHRn53UBpImmhmqXlti+ZOYTAwFPi9mXk9v3MlwMrNO7nxrSmMX7SBcw9rwv19OvhiOA6IbursnpIqAG0l1QLmmtme2IfmnIuFUbNWc9u7U9mdmc0TF3bh7G5NEh2SSyKFJgVJxxOMSk4n6DraVFI/Mxsd49icc0Vo5+4sHvxkFm/8tIQOjQ7mmYu6+doH7jei6pIK/CFnhlRJbQmqk7rHMjDnXNGZtWILN7w1mflrtjHguJbc+oe2VCxXNtFhuSQUTVIoHzlltpn9LMnXaHauGDAzBv+QzsOfzaFG5fK81v8Ijm3jgz5d/qJJCmmSBgGvhc8vASbGLiTnXFFYs3UXt787jW/nruWk9vV45NzO1K5aMdFhuSQXTVL4E3AdwUR4AkYDz8UyKOfcgRk5YxV3vT+NHbuzeOCsjlzaI8XHHrioRNP7KEPSs8BXQDZB76PdMY/MObfPtmVkct+ImbwzcRkdGx/Mkxd2pXU9n8jORS+a3kenAS8ACwjuFFpIusbMPot1cM656KWlb+DmYVNYvnEn153QihtPbEuFctHMeencL6KpPvo3cIKZzQeQ1Ar4BPCk4FwS2J2ZzVNf/czz3y6gcc1KDLvmSFKb10p0WK6YiiYprMlJCKGFQMKnuHDOweyVW7hl2FRmr9zC+d2bcM+ZHXyaa3dA8v3rkXRO+HCmpE+BYQQT2Z0PTIhDbM65fGRmZfO/0Qt58sufqV6pAi9ensrvD62f6LBcCVDQR4ozIh6vBo4PH68Fav52d+dcPCxcu41b35nK5CWb6N2pAQ+e1YlaVSokOixXQuSbFMzsyngG4pwrWHa2MeTHdB4ZOYeK5cryVN+unNmlkXc1dUUqmt5HBwH9CdZpPiin3Mz+GMO4nHMR0tdt5/Z3pzE+fQO9DqnLI+d2pv7BBxV+oHP7KJoWqdeAOcDJwP0EI5pnxzIo51wgO9sYPDadxz6fQ/myZXjsvM6c172J3x24mIkmKbQ2s/Ml9TGzIZLeBD4v7CBJLwOnE/Re6hiWPUbQVrGbYNzDlWa2SVJzgkSTM8fSODO7dt9fjnMlx8K127j93WmkLd7I79rV459nd6JBdb87cLEVzciWnLUTNknqCFQHmkdx3CvAKbnKRgEdzawz8DNwV8S2BWbWNfzyhOBKraxs46XvF3LqU9/z8+qt/Pv8Lgzql+oJwcVFNHcKAyXVBP4GjACqAn8v7CAzGx3eAUSWfRHxdBxwXvShOlfyzVm1hTvem87UpZs4sV09/nlOJ287cHEVzdxHL4UPRwMti/DafwTejnjeQtJkYAvwNzP7Pq+DJA0ABgCkpKQUYTjOJU5GZhb//WYBz30zn4MrlfeeRS5hEjL0UdJfgUzgjbBoJZBiZusldQc+kNTBzLbkPtbMBgIDAVJTUy1eMTsXKxMXb+SO96Yxf802zu7WmL+ffqiPO3AJE/ekIKkfQQP0iWZmEMzECmSEjydKWgC0BdLiHZ9z8bItI5PHP5/LkB/TaXjwQQy+8nBOOKReosNypVw04xQqhv+0CyyLhqRTgDuA481sR0R5XWCDmWVJagm0IZhjybkSadSs1fzjwxms2rKLy3o24/ZT2vmcRS4pRPNX+CNwWBRlvyJpKNALqCNpGXAPQW+jisCosK40p+vpccD9kjKBLOBaM9uwD6/DuWJh1eZd3DtiJiNnruKQ+tV49uLD6N7MZ41xyaOgCfEaAI2BSpK6EaylAHAwULmwE5vZRXkUD8pn3/eA9wqN1rliKivbeH3cYh77fC57srK5/ZRDuPrYlpQv6+sduORS0J3CycAVQBPgPxHlW4G7YxiTcyXKrBVbuGt40M302DZ1ePCsjjSrXSXRYTmXp4ImxBsCDJF0bvhJ3jm3D7bu2sMTo+bxythF1KxcgScv7Eqfrt7N1CW3aNoUPpZ0McEo5r37m9n9sQrKueLMzPh42koe+HgWa7dlcPERKdx28iHUqOzdTF3yiyYpfAhsBiYSdht1zuVt4dpt/OPDmYyZv46OjQ9m4OWpdG1aI9FhORe1aJJCEzPLPYeRcy7Czt1ZPP/tfF74biEVy5fh/j4duKRHM8qW8aoiV7xEkxTGSupkZtNjHo1zxYyZ8fnMVTzw8WyWb9rJ2d0ac1fvdtSr5vMVueKpoC6p0wnWZC4HXClpIUH1kQALZzp1rtSav2Yb9300k+/nraNdg2q8NaAnPVvWTnRYzh2Qgu4UTo9bFM4VI1t37eGZr+fz8phFVK5QlnvPOJRLezajnI85cCVAQV1SFwNIqpXH5q0xi8i5JJWdbXw4dTn//HQO67ZlcGFqU247+RBqV62Y6NCcKzLRtClMApoCGwmqjmoAKyWtAa42s4mxC8+55DB5yUbu+2gWU5ZuokvTGrx0eSpdvFeRK4GiSQojgeFm9jmApD8QrKg2DHgO6BG78JxLrFWbd/HIyDkMn7ycetUq8vj5XTinW2PKeK8iV0JFkxRSI5fHNLMvJP3TzG6R5PfNrkTauTuLgaMX8sJ3C8gy47oTWvHnXq2p4jOZuhIumr/wDZLuAN4Kn18IbJRUFsiOWWTOJYCZ8dG0lTz86WxWbN7FaZ0acuep7Whaq9A5IJ0rEaJJChcTTHv9AUGbwpiwrCxwQcwicy7OJqRv4KFPZjNl6SY6NDqYJy7sSg/vYupKmWjWaF4HXJ/P5vlFG45z8bdo3XYe+WwOI2euov7BFXnsvM6cc1gTH43sSqWCBq89aWY3SfqIYBDbr5jZmTGNzLkY27B9N09/NY/Xxy2mQrky3PL7tlx1bAsqV/B2A1d6FfTX/1r4/fH9ObGklwkGwK0xs45hWS3gbYIZV9OBC8xsY7jtLqA/wcprN+T0dnKuqO3ak8WQsek8+818tmdkcuHhKdz8+zY+NYVzFDx4bWL4/TtJlYAUM5u7D+d+BXgWeDWi7E7gKzN7WNKd4fM7JB0K9AU6AI2ALyW1NbOsfXo1zhUgK9t4f9Iynhj1Mys276LXIXW569T2HNKgWqJDcy5pFHqfLOkMgruFCkALSV2B+wurPjKz0ZKa5yruQ7BuM8AQ4FvgjrD8LTPLABZJmg8cQbAWtHMHxMz4es4aHhk5h59Xb6NT4+o8dn4Xjm5dJ9GhOZd0oqk8vZfgH/S3AGY2JY9/9tGqb2Yrw/OslFQvLG8MjIvYb1lY9huSBgADAFJSUvYzDFdaTFy8gYc/m8OE9I00r12ZZy/uRu+ODX3wmXP5iCYpZJrZ5hgvIZjXyX/TuA1gZgOBgQCpqal57uPcvNVbefTzuYyatZo6VSvy4FkdufDwppT3SeucK1A0SWFGuBxnWUltgBuAsft5vdWSGoZ3CQ2BNWH5MoL5lXI0AVbs5zVcKbZ4/Xae/HIeH0xZTpUK5fi/P7Tlj8d4jyLnohXNO+V64K8EaykMBT4HHtjP640A+gEPh98/jCh/U9J/CBqa2wDj9/MarhRasWknz3w9n3fSllK2jLj62JZce3wralXxdZGd2xfRDF7bQZAU/rovJ5Y0lKBRuY6kZQSjoh8GhknqDywBzg+vMVPSMGAWkAlc5z2PXDTWbs3guW/n88ZPSzAzLu6Rwl9OaE29g717qXP7I5reR22B/yMYW7B3fzP7XUHHmdlF+Ww6MZ/9HwIeKiwe5wA27djNwNELGfxDOhmZWZzXvQnX/66Nz1Hk3AGKpvroHeAF4CWCgWXOJcymHbsZNGYRg39IZ1tGJmd0acTNJ7WhZd2qiQ7NuRIh2t5Hz8c8EucKkDsZ9O7UgBtObEO7BgcnOjTnSpSC5j7KWYbzI0l/BoYTNDYDYGYbYhybc54MnIuzgu4UJhKMFcgZQ3BbxDYDWsYqKOc8GTiXGAXNfdQinoE4B7Bm6y4Gfb+I18ctZvvuLE8GzsWZj+hxSWH5pp0M/G4Bb01Yyp6sbE7r3IjrTmjlycC5OPOk4BJq0brtPP/tfN6ftByAcw5rzJ96taZFnSoJjsy50smTgkuIOau28N9vFvDJtBWUL1uGS3qkMOD4VjSuUSnRoTlXqkWVFCSdAxxD0MA8xsyGxzQqVyKZGeMXbeB/oxfy9Zw1VKlQlquPa0n/Y1r4AjfOJYloRjQ/B7QmmPcI4BpJJ5nZdTGNzJUYWdnGqFmreOG7hUxZuolaVSpw80lt6XdUM2pU9rmJnEsm0dwpHA90NDMDkDQEmB7TqFyJsGtPFsMnL+fF0QtZuG47KbUq80CfDpzXvSmVKpRNdHjOuTxEkxTmAinA4vB5U2BazCJyxd7mnXt446fFDP4hnbVbM+jY+GCevbgbp3RoQDlfz8C5pFbQiOaPCNoQqgOzJY0Pn/dg/9dTcCXY4vXbGfxDOu+kLWX77iyObVOHJy/sylGtahPjRZqcc0WkoDuFx+MWhSu2zIwJ6RsZNGYhX8xaTVmJM7o0ov8xLejYuHqiw3PO7aOCRjR/F89AXPGyJyubT6ev5KXvFzF9+WZqVC7Pn3u14vIjm1Pf1zJwrtiKpvfRVn5ZL7kCUB7YbmY+1LQU2rxjD2+OX8KQsems2rKLlnWr8OBZHTn3sCbeeOxcCRDNymvVIp9LOgs4IlYBueQ0Z9UWhoxNZ/jk5ezak83RrWvzz3M60qttPcqU8fYC50qKfR7RbGYfSLpzfy8o6RDg7YiilsA/gBrA1cDasPxuM/t0f6/jDlxmVjZfzFrNkLHp/LRoAxXLleGsro3pd1RzDm3kN4rOlUTRVB+dE/G0DJDKL9VJ+8zM5gJdw3OXBZYTrNVwJfCEmXkDd4Kt25bBW+OX8MZPS1i5eRdNalbirlPbceHhTX2wmXMlXDR3CmdEPM4E0oE+RXT9E4EFZrbYuywm3tSlmxgyNp2Pp61kd1Y2x7Suw/19OvK7dvUo61VEzpUK0bQpXBnD6/fll+kzAP4i6XIgDbjVzDbmPkDSAGAAQEpKSgxDKx127M7ko6kreOOnJUxbtpkqFcrS94imXH5kM1rXq1b4CZxzJYrC2St+u0F6hgKqiczshgO6sFQBWAF0MLPVkuoD68JrPgA0NLM/FnSO1NRUS0tLO5AwSq05q7bw5k9LGD5pOVszMmlbvyqX9GjGOYc1ptpB5RMdnnMuhiRNNLPUvLYVdKcQ+d/2PuCeIo0KTgUmmdlqgJzvAJJeBD4u4uuVerv2ZPHJtJW8OX4JExdvpEK5MpzWqSGX9Eihe7OaPurYOVfg4LUhOY8l3RT5vIhcRETVkaSGZrYyfHo2MKOIr1dqzV+zjTd/WsJ7k5axeeceWtapwt9Oa8+5hzWhZhVvOHbO/SLaLqn73dsoL5IqA78HrokoflRS1/Ba6bm2uX20Y3cmn0xbyTtpyxifvoHyZcXJHRpwcY8UjmzpcxE55/KWkJXXzGwHUDtX2WWJiKUkMTMmL93EsAlL+WjqCrbvzqJlnSrccUo7zk9tQp2qFRMdonMuyRU0S2rk9BaVJW3J2QSYT3ORPNZty2D4pOUMS1vKvDXbqFS+LKd1bsiFhzcl1dsKnHP7oKA2Be+PmMQys7IZPW8twyYs48vZq8nMNg5LqcHD53Ti9C6NqFrRl992zu07/89RzMxZtYXhk5bzwZTlrN6SQe0qFbjy6OZckNqUNvU9jzvnDownhWJgzdZdjJiygvcnLWfWyi2UKyN6HVKP+85szO/a1adCOV/NzDlXNDwpJKmdu7P4YtYq3p+0nO/nrSXboEuT6tx3ZgdO79yQ2t5o7JyLAU8KSSQ72xi3aD3vT1rOyBmr2JaRSeMalfhTr1ac3a0JretVTXSIzrkSzpNCgpkZM1dsYcTUFXw8dQUrNu+iasVy9O7UgLO7NaFHi1q+XoFzLm48KSTI/DVbGTF1JR9PXcHCddspV0Yc17Yud5zajj8c2sBXMXPOJYQnhThaumEHH09byYipK5i9cgsSHNmyNlcf15JTOzbwtQqccwnnSSHG1mzdxSfTVvLR1BVMWrIJgG4pNbjnjEM5rVND6vki9865JOJJIQZWb9nF5zNX8dn0Vfy0aD3ZBu0aVOP2Uw7hjM6NaFqrcqJDdM65PHlSKCLLN+3ks+krGTljFROXbMQMWteryl9OaM0ZXRr5wDLnXLHgSeEALF6/nc9mrOKz6SuZumwzAO0bHszNJ7Xl1I4NPBE454odTwr7aP6arXw2fRWfzljF7JXBHIGdm1TnjlPacWrHBjSvUyXBETrn3P7zpFCI7GxjyrJNfDlrNV/MWs38NdsA6N6sJn87rT2ndGxAk5reRuCcKxk8KeRh5+4sxsxfx5ezVvPVnDWs25ZB2TKiR4taXNazGSd3aECD6t5ryDlX8iQkKUhKB7YCWUCmmaVKqgW8DTQnWHntAjPbGK+Y1m7N4Os5qxk1aw1j5q9l155sqlUsR6929TipfT16ta1H9cq+oL1zrmRL5J3CCWa2LuL5ncBXZvawpDvD53fE6uJmxvw12xg1ezWjZq1mytJNmEHjGpXoe3gKJ7WvzxEtavkMpM65UiWZqo/6AL3Cx0OAb4lRUpi2bBPXD53M4vU7gKCh+OaT2nJS+/q0b1jNVypzzpVaiUoKBnwhyYD/mdlAoL6ZrQQws5WS6uV1oKQBwACAlJSU/bp4k5qVaVGnClcf25IT29ejYfVK+3Ue55wraWRmhe9V1BeVGpnZivAf/yjgemCEmdWI2GejmdUs6DypqamWlpYW22Cdc66EkTTRzFLz2paQCnMzWxF+XwMMB44AVktqCBB+X5OI2JxzrjSLe1KQVEVStZzHwB+AGcAIoF+4Wz/gw3jH5pxzpV0i2hTqA8PDxtxywJtmNlLSBGCYpP7AEuD8BMTmnHOlWtyTgpktBLrkUb4eODHe8TjnnPuFd8J3zjm3lycF55xze3lScM45t5cnBeecc3slZPBaUZG0FlhcBKeqA6wrdK/E8ziLlsdZtDzOohXLOJuZWd28NhTrpFBUJKXlN7ovmXicRcvjLFoeZ9FKVJxefeScc24vTwrOOef28qQQGJjoAKLkcRYtj7NoeZxFKyFxepuCc865vfxOwTnn3F6eFJxzzu1VqpKCpHRJ0yVNkfSb1XkkVZf0kaSpkmZKujJBcdaQ9K6kOZJmSzoy13ZJelrSfEnTJB2WpHFeEsY3TdJYSb+ZCDEZ4ozY73BJWZLOi3eM4fULjVNSr/Dvd6ak75IxzmR4H0k6JPw55XxtkXRTrn0S/j6KMs74vo/MrNR8AelAnQK23w08Ej6uC2wAKiQgziHAVeHjCkCNXNt7A58BAnoCPyXo51lYnEcBNcPHpyZrnGF5WeBr4FPgvGSME6gBzAJSwuf1kjTOpHgf5frdriIYsBVZnhTvoyjijOv7KFFrNCcrA6opWOyhKsEfc2Y8A5B0MHAccAWAme0GdufarQ/wqgV/JePCT24NLVzjOlniNLOxEU/HAU3iFV+OKH+eECwJ+x5weNyCixBlnBcD75vZknCfuK9OGGWcCX8f5XIisMDMcs9+kPD3US55xhnv91Gpqj4i+GP9QtJESQPy2P4s0B5YAUwHbjSz7HgGCLQE1gKDJU2W9FK4Ql2kxsDSiOfLwrJ4iibOSP0JPpXFW6FxSmoMnA28kID4ckTz82wL1JT0bfg3fHn8w4wqzmR4H0XqCwzNozwZ3keR8oszUszfR6UtKRxtZocR3IJdJ+m4XNtPBqYAjYCuwLPhJ6N4KgccBjxvZt2A7cCdufZRHsfFu29xNHECIOkEgj/mO+IX3l7RxPkkcIeZZcU5tkjRxFkO6A6cRvC3+ndJbeMaZXRxJsP7CABJFYAzgXfy2pxHWUL66BcSZ84+cXkflaqkYGYrwu9rgOHAEbl2uZLg9tzMbD6wCGgX3yhZBiwzs5/C5+8SvAlz79M04nkTgk9l8RRNnEjqDLwE9LFgdb14iybOVOAtSenAecBzks6KW4SBaH/vI81su5mtA0aTxyqGMRZNnMnwPspxKjDJzFbnsS0Z3kc5Coozru+jUpMUJFWRVC3nMfAHYEau3ZYQLgkqqT5wCLAwnnGa2SpgqaRDwqITCRoXI40ALg97T/QENse7HjSaOCWlAO8Dl5nZz/GML0c0cZpZCzNrbmbNCf7J/dnMPki2OIEPgWMllZNUGegBzI5jmNHGmfD3UYSLyL9KJuHvowj5xhn391EiW9vj+UVQFzo1/JoJ/DUsvxa4NnzcCPiCoB50BnBpgmLtCqQB04APgJq54hTwX2BBGGtqksb5ErCRoCphCpCWjHHm2vcVEtf7qNA4gdsI/gnPAG5KxjiT6H1UGVgPVI8oS8b3UWFxxvV95NNcOOec26vUVB8555wrnCcF55xze3lScM45t5cnBeecc3t5UnDOObeXJwWXJ0lPRM7WKOlzSS9FPP+3pFsknSkpz5HMBZz7lbxmIpXUU9JP4WyRsyXdeyCvIYo4rpD07D7s31VS73y29ZL0ca6yPF9nrn3ulfR/UVx7W67nhcYeze8mr7gjtt0UjofI79h3JbUs6Pz5HDcxHMGb17bHJf1uX8/pio4nBZefsQSzMyKpDFAH6BCx/SjgBzMbYWYPF9E1hwADzKwr0BEYVkTnLSpdCWbWLBaK4HdzE0Ef+t+Q1AEoa2b7NChNUnNguQUT6eXlGfKZLsXFhycFl58fCJMCQTKYAWyVVFNSRYIJzyZHfmINPxk/Hc75vjDnU3I4YvRZSbMkfQLUy+ea9YCVAGaWZWazwuPvlfSapK8lzZN0dc4Bkm6TNEHBXPP3RZRfKml8eNfxP0llw/IrJf2sYC2Co/MKQtIR4WuYHH4/JPxkez9wYXjOC/flh6lgLY864eNUSd9GbO6S12vbh3PXlfRe+HOYIOnosDzyd9NK0rhw+/257jyq6pf1Ed4If183EAxC+0bSN3lc9hKCEdY5MWyT9Eh4F/Bl+DP8Nvw7ODPiuFOBkZLKhn8vMxSscXIzgAUzhNaW1GBffw6uaPjU2S5PZrZCUqaCIfZHAT8SzCB5JLAZmGZmu6XfzCnWEDiGYK6bEQTTRpxNMNVBJ6A+wYjcl/O47BPA3PAf5khgiJntCrd1JpjzvgpBMvqE4G6iDcEcVgJGKJjkcC1wIcEEiHskPQdcImkUcB/BpHKbgW+AyXnEMQc4zswyJZ0E/NPMzpX0D4JRr3/J58d2rKQpEc9TgDyrZnL5zWuzcJ6uCJVynbsWwc8X4CngCTMbE/6+PidI2pGeAp4ys6GSrs21rRtB4l9B8GHgaDN7WtItwAkWzLOU29H8elqGKsC3ZnaHpOHAg8DvgUMJ7gBzYj0FuJngrquxmXWEYOGeiHNNCs//Xh7XdTHmScEVJOdu4SjgPwRJ4SiCf6hj8znmAwumSZ6lYN4bCObfH2rBLKQrJH2d14Fmdr+kNwjmpbqYYD6YXuHmD81sJ7Az/OR6BEHy+QO//GOvSpAkOhP8458QJq1KwBqCuYK+NbO1AJLeJpiOOrfqwBBJbQhmzSyf3w8ol+/N7PScJ5JeifK4vF7bB7n22RlWq+Wc+wqCifwATgIOjUjQByuc5yvCkcBZ4eM3gccjto03s2XheacAzYExhcTckCD55thNkMghmDIiI0zI08Pz5cwE2sTMFkqqCbSU9AzwCcG0GDnWENyluATwpOAKktOu0Img+mgpcCuwhbw/6QNkRDyOvI2Iaj4VM1sAPC/pRWCtpNr5HG/h+f9lZv+L3CDpeoK7jLtylZ8VZRwPAN+Y2dlhHfi30cReiEx+qa49KNe2vF7bvigDHBkmlr3yuIvLT+TvLIvo/i/s5NevY4/9MmdOds45zSxbUs75jiVMNma2UcGykicD1wEXAH8M9zsoPL9LAG9TcAX5ATgd2BDW8W8gWBLySILqpGiNBvqG9cgNgRPy2knSafrlP1kbgn9Qm8LnfSQdFCaJXsAEgmqSP0qqGh7fWFI94CvgvPAxkmpJagb8BPSSVFtSeeD8fOKtDiwPH18RUb4VyP0JPFrpBHcvAOfm2pbXa9sXXwB7q7Qkdc1jn3ER1+0b5XkLer2zgdZRnifHKYQLxITtK2XM7D3g7/x6+u22/HYGYxcnnhRcQaYT9Doal6tscz71zPkZDswLj30eyG/B+csI2hSmAK8Bl9gvC9+MJ6hmGAc8YGYrzOwLgqqQH8NqineBamED9d8IVtmbBowCcpZZvJcgoX1JUHedl0eBf0n6gWDd3BzfEFTT7HNDM0FbxlOSvidIdpF+89r28dw3AKkKGttnEcywmdtNwC2SxhNU/WyO4rwDgc/yaWj+hF+q9qLVi19+942Bb8Pf9SvAXQBhsm5NMAurSwCfJdUlPQXjFbaZ2eOF7evypmC8wU4zM0l9gYvMrM8BnK8SQZI82qJYsU5SE+BFMzu1kP3OBg4zs7/vb2zuwHibgnOlQ3eCZTFFUCX3x4J3L5iZ7ZR0D8En/iVR7L+MoDtqYcoB/z6Q2NyB8TsF55xze3mbgnPOub08KTjnnNvLk4Jzzrm9PCk455zby5OCc865vf4fRqhUWGTIPWQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEGCAYAAACZ0MnKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA2BElEQVR4nO3dd3xUVfr48c+TRggEkkASQg0gvUNooohiwYquin2xIi676u6qq+t+Ld/9fV0X265rxYarYBdFV0RAigWEhN6lBOkJvQeSPL8/7skyxiRMQiYzkzzv12tec+fMPXOfKfc+c8+99xxRVYwxxhh/RAQ7AGOMMeHDkoYxxhi/WdIwxhjjN0saxhhj/GZJwxhjjN+igh3AyWjYsKGmp6cHOwxjjAkrWVlZO1Q1uSJ1wzpppKenk5mZGewwjDEmrIjIhorWteYpY4wxfgtY0hCRZiIyXURWiMgyEbnLlT8iIptFZKG7XeBT5wERWSMiq0TkvEDFZowxpmIC2TyVD/xRVeeLSDyQJSJT3HPPqOqTvjOLSEfgaqAT0BiYKiJtVbUggDEaY4wph4DtaajqVlWd76b3AyuAJmVUGQq8q6p5qroeWAP0CVR8xhhjyq9KjmmISDrQA/jBFf1WRBaLyOsikujKmgAbfaptooQkIyIjRCRTRDJzc3MDGbYxxphiAp40RKQu8BFwt6ruA14EWgPdga3AU0WzllD9F70pquoYVc1Q1Yzk5AqdMWaMMaaCApo0RCQaL2GMU9WPAVR1u6oWqGoh8ArHm6A2Ac18qjcFtgQyPmOMMeUTyLOnBHgNWKGqT/uUp/nMdhmw1E1PBK4WkVoi0hJoA8wNRGx7Dh3l0c+WsffwsUC8vDHGVFuBPHtqAHADsEREFrqyPwPXiEh3vKanbOB2AFVdJiLvA8vxzrwaFagzp37adYg3v8/m8NECHr+8ayAWYYwx1VLAkoaqfkvJxym+KKPO/wH/F6iYinRtmsBtA1vx8sx1XNS1Mae1aRjoRRpjTLVQY68I//3ZbWnVsA73f7yYg3n5wQ7HGGPCQo1NGrHRkYy+oiub9xzmicmrgh2OMcaEhRqbNAAy0pMY3j+dsd9n88O6ncEOxxhjQl6NThoA9w1pR/OkOO79cDGHjlozlTHGlKXGJ424mChGX9GVn3YdYvSX1kxljDFlqfFJA6BfqwbceKrXTDXHmqmMMaZUljSc+4a0o0WDOO79cBEH7GwqY4wpkSUNJy4miqeu7Mam3Yd57IsVwQ7HGGNCkiUNHxnpSYw4vRXjf/iJGatygh2OMcaEHEsaxfz+nLa0Ta3Lnz5azN5D1jeVMcb4sqRRTGx0JE8P687OA0f5n0+XnriCMcbUIJY0StC5SX3uGtyGiYu28OnCzcEOxxhjQoYljVLcMag1PZon8D+fLGXr3sPBDscYY0KCJY1SREVG8Myw7uQXKvd8sIjCwl8MImiMMTWOJY0ypDesw/9c1JHv1uzk9e/WBzscY4wJOksaJ3B172ac0zGV0V+uYvmWfcEOxxhjgiqQw702E5HpIrJCRJaJyF2u/AkRWSkii0VkgogkuPJ0ETksIgvd7aVAxVYeIsLfL+9K/bho7np3AUeOBWQwQWOMCQuB3NPIB/6oqh2AfsAoEekITAE6q2pXYDXwgE+dtara3d1GBjC2ckmqE8NTV3bjx5wDdrW4MaZGC1jSUNWtqjrfTe8HVgBNVPUrVS3q3GkO0DRQMVSmgW2TueW0lvx79gamLt8e7HCMMSYoquSYhoikAz2AH4o9dTMwyedxSxFZICIzReT0Ul5rhIhkikhmbm5uYAIuxX1D2tExrR73friI7fuOVOmyjTEmFAQ8aYhIXeAj4G5V3edT/iBeE9Y4V7QVaK6qPYA/AONFpF7x11PVMaqaoaoZycnJgQ7/Z2pFRfLsNT04cqyQu99dSIGdhmuMqWECmjREJBovYYxT1Y99yocDFwHXqaoCqGqequ5001nAWqBtIOOriFNS6vLoJZ2YvW4nL81cG+xwjDGmSgXy7CkBXgNWqOrTPuVDgD8Bl6jqIZ/yZBGJdNOtgDbAukDFdzKuzGjKRV3TeHrKarI27Ap2OMYYU2UCuacxALgBOMvnNNoLgOeAeGBKsVNrBwKLRWQR8CEwUlVDcossIjz2qy40Tojld+MXsOfQ0WCHZIwxVUJc61BYysjI0MzMzKAtf9HGPVzx0vec0TaFV37dC2/nyhhjQpuIZKlqRkXq2hXhJ6FbswTuP78DU1ds543vsoMdjjHGBJwljZN084B0zu6Qyt8mrWDhxj3BDscYYwLKksZJEhGevLIrKfGxjBo3345vGGOqNUsalSAhLobnr+tJzv4j/PF960bdGFN9WdKoJN2bJfDgBR2YtjKHV74JyTOFjTHmpFnSqETDT03nwi5pjJ68ijnrdgY7HGOMqXSWNCqRiPD45V1o0SCO345fYP1TGWOqHUsalSw+NpqXru/Fwbx8Ro2bz7GCwmCHZIwxlcaSRgC0TY3n8cu7kLlhN49PWhnscIwxptJY0giQod2bcOOp6bz27XomLtoS7HCMMaZSWNIIoD9f0IGMFon86cPFrNhq44sbY8KfJY0AiomK4IXrexIfG8Xtb2XZhX/GmLBnSSPAUuJjefH6Xmzde5i7bOAmY0yYs6RRBXq1SOThizsxc3UuT361KtjhGGNMhUUFO4Ca4rq+zVm2ZS8vzlhLx7R6XNytcbBDMsaYcrM9jSoiIjx6SWcyWiRy74eLWLp5b7BDMsaYcgvkcK/NRGS6iKwQkWUicpcrTxKRKSLyo7tP9KnzgIisEZFVInJeoGILlpioCF68vheJcTHc/lYWOw7kBTskY4wpl0DuaeQDf1TVDkA/YJSIdATuB6apahtgmnuMe+5qoBMwBHihaMzw6iQ5vhZjbshgx4E87ng7i7z8gmCHZIwxfgtY0lDVrao6303vB1YATYChwJtutjeBS930UOBdVc1T1fXAGqBPoOILpi5N6/PEld2Yl72bv0xYSjgPuWuMqVmq5JiGiKQDPYAfgFRV3QpeYgFS3GxNgI0+1Ta5smrpkm6NufOsU/ggaxOvfrM+2OEYY4xfAp40RKQu8BFwt6qWdVm0lFD2i7/gIjJCRDJFJDM3N7eywgyKu89uy/mdG/HYpBVMW7E92OEYY8wJBTRpiEg0XsIYp6ofu+LtIpLmnk8Dclz5JqCZT/WmwC86bVLVMaqaoaoZycnJgQu+CkRECE8N60bHtHrc+c4C62rEGBPyAnn2lACvAStU9WmfpyYCw930cOBTn/KrRaSWiLQE2gBzAxVfqIiLieLV4RnUjY3ilrHzyLExOIwxISyQexoDgBuAs0RkobtdADwOnCMiPwLnuMeo6jLgfWA58CUwSlVrxKlFafVr89rw3uw+dIxb/53J4aM14m0bY8KQhPOZOxkZGZqZmRnsMCrNV8u2cfvbWZzXsREvXNeTiIiSDvMYY8zJEZEsVc2oSF27IjyEnNupEQ9e0IEvl23jb5NWBDscY4z5Bet7KsTcclpLNuw8xCvfrKdpYhzDT00PdkjGGPNfljRCjIjw8MUd2br3MI9+toy0+rGc26lRsMMyxhjAmqdCUlRkBM9e04POTepz57sLWLhxT7BDMsYYwJJGyIqLieK14b1pWLcWt4ydR/aOg8EOyRhjLGmEsuT4Wrx5cx8KVRn+xlxy91uvuMaY4LKkEeJaJ9fltRt7s33fEW4eO4+DefnBDskYU4NZ0ggDPZsn8vy1PVm+dR8j387iaH5hsEMyxtRQljTCxOAOqfztsi588+MO7vtwEYWF4XtRpjEmfNkpt2FkWO9m5B7I44nJq0iIi+HhizvidfFljDFVw++kISJ1gCM1pT+oUPWbQa3ZdfAor327ngZ1Yvjd4DbBDskYU4OUmjREJAJv+NXrgN5AHlBLRHKBL4AxqvpjlURp/ktEePCCDuw+dJSnpqwmoU4MN/RrEeywjDE1RFl7GtOBqcADwFJVLQQQkSTgTOBxEZmgqm8HPkzjKyJC+PvlXdl3+BgPfbqUerFRDO1ebQc5NMaEkLKSxtmqeqx4oaruwhtY6SM3yJIJgujICJ67tifDX5/LH99fRJ2YKM7umBrssIwx1VypZ0/5JgwRSRSRriLSs+hWfB5T9WKjI3l1eAadGtfjN+Pn892aHcEOyRhTzZ3wlFsR+SuwGHgWeMrdngxwXMZP8bHRjL2pDy0b1OG2f2eStWF3sEMyxlRj/lynMQxoraqDVPVMdzsr0IEZ/yXWieGtW/qQEl+Lm96Yy9LNe4MdkjGmmvInaSwFEsr7wiLyuojkiMhSn7L3fIZ+zRaRha48XUQO+zz3UnmXV9Ol1Ivl7Vv7Eh8bzfWv/cDKbfuCHZIxphryJ2n8DVggIpNFZGLRzY96Y4EhvgWqepWqdlfV7ngH0z/2eXpt0XOqOtLP+I2PpolxjL+tL7FRkVz3yg+sydkf7JCMMdWMP0njTeDvwOMcP6bx1IkqqeosYFdJz4l3GfMw4B2/IzV+adGgDuNu64uIcO0rP7DeulQ3xlQif5LGDlV9VlWnq+rMottJLvd0YHuxiwNbisgCEZkpIqeXVlFERohIpohk5ubmnmQY1VPr5LqMv60v+YXKta/MYcNOSxzGmMrhT9LIEpG/iUj/4qfcnoRr+Plexlaguar2AP4AjBeReiVVVNUxqpqhqhnJycknGUb11TY1nnG39uXIsQKuGWOJwxhTOfxJGj2AfsBjVMIptyISBfwKeK+oTFXzVHWnm84C1gJtK7oM4+mQVo9xt/bjkCUOY0wlOWHS8DnN9sxKOuX2bGClqm4qKhCRZBGJdNOtgDbAupNYhnE6Nq7HeEscxphKUmrSEJHrXaeFpT3fWkROK+P5d4DZQDsR2SQit7inruaXB8AHAotFZBHwITDSdVdiKoFv4rh6zBwbb9wYU2GiWvJgPiJyF3AzkOVuuUAscApwBrADuD+YPd1mZGRoZmZmsBYfdpZv2cf1r/1AdKQw7tZ+nJJSN9ghGWOCQESyVDWjInXL6nvqn0BPvL2CZGCwe7wZuEFVL7eu0cNLx8b1eOe2fhQUwtVj5rBqm13HYYwpn1L3NMKB7WlUzJqcA1z36hyO5hfy9q196dS4frBDMsZUoYDsaZjq65SUurw3oj+1oyO59pUfWLhxT7BDMsaECUsaNVR6wzq8d3t/6teO5rpX5jB77c5gh2SMCQP+dI3e0p8yE36aJcXxwcj+NE6ozY1vzOXrlduDHZIxJsT5s6fxUQllH1Z2ICY4UuvF8t7t/WmbGs+If2fx2aItwQ7JGBPCSh3uVUTaA52A+iLyK5+n6uGdemuqiaQ6MYy/rS+3jM3kzncXcCAvn2v6NA92WMaYEFTWGOHtgIvwxtK42Kd8P3BbAGMyQRAfG82bN/fhjnFZPPDxEnYdPMpvBrXG65DYGGM8pSYNVf0U+FRE+qvq7CqMyQRJ7ZhIXvl1Bvd8sIgnJq9i54Gj/OXCDkREWOIwxnjK2tMoskZE/gyk+86vqjcHKigTPNGRETwzrDuJcTG8/t16dh86yugruhIdaSfaGWP8SxqfAt8AU4GCwIZjQkFEhPDwxR1pWDeGJ79azZ5DR3n+up7ExfjzczHGVGf+bAXiVPVPAY/EhBQR4bdntSGpTi3+8skSrhkzh9du7E3DurWCHZoxJoj8aXP4XEQuCHgkJiRd27c5L9+Qwart+7n8xe+th1xjariyuj7fLyL7gLvwEsdhEdnnU25qiHM6pjLu1n7sO3yMy1/8nkXW7YgxNVZZvdzGq2o9dx+hqrV9Hpc4FKupvnq1SOSjO04lrlYkV4+Zw7QVdvW4MTWRP92I9Czh1toN22pqkFbJdfnojlNpnVKH2/6dyZvfZwc7JGNMFfPnmMYLwBzgFXebA7wLrBaRc0urJCKvi0iOiCz1KXtERDaLyEJ3u8DnuQdEZI2IrBKR8yr8jkxApcTH8v7t/TmrfSoPT1zG/362nILC8O1e3xhTPv4kjWygh6r2UtVeQHdgKd5Y36PLqDcWGFJC+TOq2t3dvgAQkY54w8B2cnVeKBoz3ISeuJgoXr6hFzcNSOf179Yz8u0sDh3ND3ZYxpgq4E/SaK+qy4oeqOpyvCSyrqxKqjoL8Hec76HAu6qap6rrgTVAHz/rmiCIjBAevrgTj1zckWkrtnPVy3PYvu9IsMMyxgSYP0ljlYi8KCJnuNsLeE1TtYBjFVjmb0VksWu+SnRlTYCNPvNscmW/ICIjRCRTRDJzc3MrsHhTmW4c0JJXfp3B2twDXPLctyzZtDfYIRljAsifpHEj3j//u4HfA+tc2THgzHIu70WgNV4T11bgKVdeUudGJTaUq+oYVc1Q1Yzk5ORyLt4EwuAOqXx0x6lERURw5cvf85/FW4MdkjEmQE6YNFT1sKo+paqXqeqlqvqkqh5S1UJVPVCehanqdlUtUNVCvIPqRU1Qm4BmPrM2BWxghzDSIa0en4waQKfG9Rk1fj7/nPoj4Tz+vDGmZGVd3Pe+u1/impN+dqvIwkQkzefhZXgH1AEmAleLSC03KmAbYG5FlmGCJzm+FuNu7cuvejThmamr+d07Czh81LorM6Y6Ketai7vc/UUVeWEReQcYBDQUkU3Aw8AgEemO1/SUDdwOoKrLXJJaDuQDo1TVtjZhKDY6kqeGdaNNajyjJ69kXe5Bxvy6F00T44IdmjGmEog/TQgi0gJoo6pTRaQ2EKWq+wMe3QlkZGRoZmZmsMMwpZi+Moc7311AdGQEz1/bk/6tGwQ7JGMMICJZqppRkbr+XBF+G96Y4C+7oqbAJxVZmKlZzmyfwqejBpAYF831r/3Am99n23EOY8KcP2dPjQIGAPsAVPVHICWQQZnqo1VyXT4ZNYAz2yXz8MRl3PfhYo4cs5ZHY8KVP0kjT1WPFj1wfU7Z30Xjt/jYaMbckMGdg9vwQdYmrnjpezbuOhTssIwxFeBP0pjphnutLSLnAB8AnwU2LFPdREQIfzinLa/+OoMNOw9x8XPfMmNVTrDDMsaUkz9J434gF1iCd7bTF8BfAhmUqb7O7pjKZ789jUb1Yrlp7DyenfYjhdbhoTFhw6+zp0KVnT0Vvg4fLeDPE5YwYcFmzmqfwtPDupEQFxPssIypEU7m7KlSr9MQkSWUcexCVbtWZIHGANSOieTpYd3o2TyB//18ORc++y3/urYHPZsnnriyMSZoymqeugi4GLgEiHTTvjdjToqIcEP/dD6641QiImDYS7N59Zt1dlquMSGsrOFeN7hbNt4ZVBt8b1UXoqnuujZN4PPfnc7gDin8v/+sYMRbWew9VJEOlI0xgebPgXBjAq5+7Wheur4XD13UkRmrcrjwX9+w4KfdwQ7LGFNMWR0W/ndMcLzTbXsUKzOmUokIN5/Wkg9GnooqXPnSbF6csdbOrjImhJR69pSITC+jnqrqWYEJyX929lT1tffwMf788RL+s2Qrp53SkKeHdSOlXmywwzKmWjiZs6fslFsTslSV9zM38sjE5dSOieTJK7tyVvvUYIdlTNgLaIeFxgSLiHBV7+Z89rsBpNaL5eaxmTwycZn1XWVMEFnSMCHvlJR4JvzmVG48NZ2x32dzyXPfsnzLvmCHZUyNZEnDhIXY6EgeuaQTY2/qze5Dxxj6/Le8NHMtBXaQ3Jgq5VfSEJFficjTIvKUiFzmZ53XRSRHRJb6lD0hIivdkLETRCTBlaeLyGERWehuL1Xo3Zhqb1C7FCbfPZDB7VN5fNJKrn1lDpt2W4+5xlQVfwZhegEYiddh4VLgdhF53o/XHgsMKVY2BejsuiBZDTzg89xaVe3ubiP9Cd7UTEl1Ynjx+p48cUVXlm7ey/n/+IaP52+yK8mNqQL+7GmcAZynqm+o6hvABXhjf5dJVWcBu4qVfaWq+e7hHLxRAI0pNxHhyoxmfHn3QNo1iucP7y/i9reyyN2fF+zQjKnW/Ekaq4DmPo+bAYsrYdk3A5N8HrcUkQUiMlNETi+tkoiMEJFMEcnMzc2thDBMOGuWFMd7t/fnzxe0Z8bqXM59ZiafL94S7LCMqbbKuiL8MxGZCDQAVojIDHfB3wog+WQWKiIPAvnAOFe0FWiuqj2APwDjRaReSXVVdYyqZqhqRnLySYVhqonICGHEwNZ8cedpNE+K47fjFzBq3Hx2HrC9DmMqW6ldowNPBmKBIjIcrwfdweoaoVU1D8hz01kishZoC9iVe8Zvp6TE89Edp/LyrHX8Y+pq5qzbyf+7tDPnd0kLdmjGVBulJg1VnVnZCxORIcCfgDNU9ZBPeTKwS1ULRKQV0AZYV9nLN9VfVGQEo848hcEdUrjng0XcMW4+53duxKNDO5ESb92QGHOy/Dl7ar+I7HO3IyJSICInvLJKRN4BZgPtRGSTiNwCPAfEA1OKnVo7EFgsIouAD4GRqrqrxBc2xg/tG9Vjwm8GcO957Zi2Modznp7FB5kb7QwrY05SufueEpFLgT6q+ueARFQO1veU8cfa3APc/9Fi5mXv5vQ2DXnssi40S4oLdljGBE2V9j2lqp8AQe/h1hh/tU6uy3sj+vPXSzszf8Nuzn1mFq99u96uJjemAso6EA54V4P7PIwAMihj7HBjQlFEhHBDvxYMbp/CXz5Zyl8/X84nCzbz2GVd6NK0frDDMyZs+LOn4Tsu+HnAfmBoIIMyJlAaJ9TmteEZPHdtD7btO8LQ57/lkYnL2H/Ehpc1xh82noapsfYdOcaTk1fx1pwNpMTX4pGLOzGkcyNEJNihGRNQARmESUT+RRnNUKp6Z0UWWJksaZjKsHDjHv788RKWb93HWe1TePSSTnag3FRrgToQnglkudslPtNFN2Oqhe7NEpj42wH85cIOzFm3k7Ofnsk/p/5ogz0ZUwK/mqdEZIHr4iOk2J6GqWxb9x7m//6zgs8Xb6VZUm0euqgTZ3dIsSYrU61UxSm34Xvgw5hySKtfm+eu7cn42/oSGxXJbf/O5Kax81i/42CwQzMmJNjIfcaU4NTWDfnirtP5y4UdyMzezXnPzOKJySs5dDT/xJWNqcbKOhC+n+N7GHFAUV9RAqiqltgLbVWy5ilTFXL2H+HxSSv5eP5mUuvV4r7z2nNZjyZERFiTlQlPATl7KhxY0jBVKWvDbv738+Us2riHrk3r89BFHclITwp2WMaUW5V2I2JMTdWrRSIT7jiVf1zVnZx9eVzx0mxGjZ/Pxl02RrmpOSxpGFMOERHCpT2a8PU9Z3D32W2YtmI7g5+eyegvV3Igz453mOrPkoYxFRAXE8XdZ7dl+j2DuLBLGi/MWMsZo6fz5vfZHM0vDHZ4xgSMJQ1jTkJa/do8c1V3Ph01gDapdXl44jLOeWYm/1m81cbuMNWSJQ1jKkG3Zgm8c1s/3rixN7WiIhg1fj6XvvA9c9btDHZoxlSqgCUNEXldRHJEZKlPWZKITBGRH919os9zD4jIGhFZJSLnBSouYwJFRDizfQqT7hrI6Cu6sn3vEa4eM4ebx85j1bb9wQ7PmEoRyD2NscCQYmX3A9NUtQ0wzT1GRDoCVwOdXJ0XRCQygLEZEzCREcKwjGbMuHcQ9w1px7z1uxjyz1nc/e4Csu3KchPmApY0VHUWUHyc76HAm276TeBSn/J3VTVPVdcDa4A+gYrNmKoQGx3Jbwadwqz7zmTEwFZ8uWwbg5+eyQMfL2bLnsPBDs+YCqnqYxqpqroVwN2nuPImwEaf+Ta5sl8QkREikikimbm5uQEN1pjKkFgnhgfO78Cse8/k+r7N+TBrE4OemMEjE5eRuz8v2OEZUy6hciC8pP4YSjz1RFXHqGqGqmYkJycHOCxjKk9KvVgeHdqZ6fcM4tIejXlrzgYGjp7O379cyZ5DR4MdnjF+qeqksV1E0gDcfY4r3wQ085mvKbClimMzpko0TYxj9BXdmPL7gZzTMZWXZq5lwONfM/rLlew6aMnDhLaqThoTgeFuejjwqU/51SJSS0RaAm2AuVUcmzFVqlVyXZ69pgeT7jqdQe1TeHHmWk77+9f8bdIKdh6wZisTmgLWYaGIvAMMAhoC24GHgU+A94HmwE/Alaq6y83/IHAzkA/craqTTrQM67DQVCert+/nua/X8NniLcRGRXJ9v+aMGNia5PhawQ7NVDPWy60x1cianAM8P30Nny7cTExUBNf2acHIM1qRUi822KGZasKShjHV0PodB3l++homLNhMpAiX92rK7QNbkd6wTrBDM2HOkoYx1dhPOw/x8qy1fJC1ifyCQi7oksbIM1rTuUn9YIdmwpQlDWNqgJz9R3jju2zenr2B/Xn5DGybzB1ntKZfqyREbBRB4z9LGsbUIPuOHOPtORt4/dtsdhzIo3uzBO4Y1JpzOqTaELTGL5Y0jKmBjhwr4MOsTYyZtY6fdh2iVcM63HRaSy7v2YS4mKhgh2dCmCUNY2qw/IJCJi3dxqvfrGPRpr0kxEVzXd/m/Lp/Oql2xpUpgSUNYwyqStaG3bz6zXq+Wr6NyAjh4q6Nufm0lnbQ3PzMySQN24c1ppoQETLSk8hIT+KnnYd44/v1vD9vIx8v2Ey/VkncelorzmqfYsc9zEmxPQ1jqrG9h4/x3ryfGPtdNlv2HqFFgziu79uCKzOakhAXE+zwTJBY85QxpkzHCgr5cuk23pq9gbnZu4iNjmBotybc0L+FNV3VQJY0jDF+W75lH2/N2cAnCzZz+FgBvVok8uv+LTi/cxoxUaEyWoIJJEsaxphy23v4GB9mbeKt2dlk7zxEw7q1uLZPM67t24JG9e2sq+rMkoYxpsIKC5VZP+by1uwNfL0qhwgRzmyXwrV9m3FG2xQi7cB5tWNnTxljKiwiQhjULoVB7VLYuOsQ4374iQ+zNjF1xXbS6scyLKMZw3o3o0lC7WCHakKA7WkYY37hWEEh01ZsZ/zcjXzzYy4Ag9omc3Wf5pzVPoXoSDv2Ec6secoYEzAbdx3i/cyNvJ+5ke378kiOr8WwjKZcldGc5g3igh2eqQBLGsaYgMsvKGTGqlzemfsT01flUKjQr1USl/dsygVd0qhTy1q7w0VYJQ0RaQe851PUCngISABuA3Jd+Z9V9YuyXsuShjHBsXXvYT7M3MRH8zeRvfMQcTGRDOnciCt6NqVfqwZ21XmIC6uk8bOFi0QCm4G+wE3AAVV90t/6ljSMCa6i/q4+mr+JzxdtZX9ePk0SavOrnk24vGdTG2UwRIXz2VODgbWqusEGkTEm/Pj2d/XwxZ2YvGwbH83fzPPT1/Cvr9eQ0SKRy3s15cKuadSLjQ52uKYSBHtP43Vgvqo+JyKPADcC+4BM4I+quruEOiOAEQDNmzfvtWHDhqoL2Bjjl217jzBhwWY+mr+JNTkHiImKYHD7FC7p1pgz26cQGx0Z7BBrtLBsnhKRGGAL0ElVt4tIKrADUOCvQJqq3lzWa1jzlDGhTVVZvGkvExZs5vPFW9lxII+6taI4r1MjLunemAGtGxBlp+9WuXBNGkOBUap6bgnPpQOfq2rnsl7DkoYx4SO/oJA563bx6cLNfLlsG/uP5NOgTgwXdEljaPfG9GyeaAfQq0i4Jo13gcmq+oZ7nKaqW93074G+qnp1Wa9hScOY8JSXX8CMVblMXLSFaSu2c+RYIU0SanNRtzQu6daYjmn1sOOcgRN2SUNE4oCNQCtV3evK3gK64zVPZQO3FyWR0ljSMCb8HcjLZ8rybUxcuIVvftxBfqHSqmEdhnRuxAVd0ujU2BJIZQu7pFFZLGkYU73sOniUSUu3MmnJNmav20lBodI0sTbnd27EkM5p9GiWYE1YlcCShjGm2tl18ChTl29n0tKtfLtmB8cKlEb1YhnSuRFDOjeid3qS9cBbQZY0jDHV2t7Dx/h65XYmLdnGzNW55OUX0rBuDOd2asT5nRvRr1UD60SxHCxpGGNqjIN5+UxflcOkpduYvjKHQ0cLqF87mjPbJTO4QypntEu2CwlPwJKGMaZGOnKsgJmrc5m8zEsguw8dIypC6NsqicHtUzm7Q6r1xFsCSxrGmBqvoFBZ8NNupqzYzrQVOazJOQBA29S6nN0hlcEdUuneLMGOg2BJI9hhGGNCUPaOg0x1CWRu9i4KCpUGdWI4q30KgzukcnqbhjW2O3dLGsYYU4a9h44xY3UOU1fkMGNVDvuP5BMTFUHflkmc0TaZM9omc0pK3RpzPYglDWOM8dOxgkLmZe9i2oocZq7O/W8zVlr9WM5om8zAtskMOKUh9WtX34PpljSMMaaCNu85zKzVucxclct3a3awPy+fyAihR7MEBrq9kC5N6leriwotaRhjTCU4VlDIwo17vCSyOpclm/eiColx0Zzexksgp7dtSEp8bLBDPSmWNIwxJgB2Hsjj2zU7mLkql1k/5rLjwFEA2jeK59TWDTm1dQP6tEoKu+tCLGkYY0yAFRYqy7fuY+bqXL5fu4PM7N3k5RcSIdClaQL9WzXg1NYNyEhPJC4mtM/KsqRhjDFV7MixAhb8tIfZa3cwe91OFvy0h/xCJTpS6NEskf6tG9C/dQN6NE+gVlRojVRoScMYY4LsYF4+mRt28/3aHcxeu5Olm/dSqBAbHUFGiyT6t/b2RLo0qR/00QpPJmmE9j6UMcaEiTq1ov57zQd4nSzOXb/rv0nkicmrAKhbK4peLRLpnZ5IRnoS3ZslhNWY6banYYwxVWDngTzmrPOSyLzsXaze7l0fEh0pdGlSn97pSWSkJ5HRIpHEOjEBjSXsmqdEJBvYDxQA+aqaISJJwHtAOt7IfcNUdXdZr2NJwxgTrvYcOkrWht3My95NZvYuFm/ay9GCQgDapNQlIz2J3umJ9E5Pomli7Uq9Wj1ck0aGqu7wKRsN7FLVx0XkfiBRVf9U1utY0jDGVBdHjhWweNNe5mXvIjN7F5kbdrP/SD4AjerFkuESSO/0JNo1ij+pjheryzGNocAgN/0mMAMoM2kYY0x1ERsdSZ+WSfRpmQR4p/iuztnPvPW7mJe9m3nZu/h88VYA4mtFcXWfZjx4YccqjzNYSUOBr0REgZdVdQyQqqpbAVR1q4iklFRRREYAIwCaN29eVfEaY0yViogQ2jeqR/tG9bihfzrgdXmSmb2Luet3kVa/dlDiClbzVGNV3eISwxTgd8BEVU3wmWe3qiaW9TrWPGWMMeV3Ms1TQTlZWFW3uPscYALQB9guImkA7j4nGLEZY4wpXZUnDRGpIyLxRdPAucBSYCIw3M02HPi0qmMzxhhTtmAc00gFJrjTx6KA8ar6pYjMA94XkVuAn4ArgxCbMcaYMlR50lDVdUC3Esp3AoOrOh5jjDH+C24HKMYYY8KKJQ1jjDF+s6RhjDHGb5Y0jDHG+C2se7kVkVxgQxUtriGw44RzhYZwihXCK16LNTDCKVYIr3hLirWFqiZX5MXCOmlUJRHJrOgVlFUtnGKF8IrXYg2McIoVwiveyo7VmqeMMcb4zZKGMcYYv1nS8N+YYAdQDuEUK4RXvBZrYIRTrBBe8VZqrHZMwxhjjN9sT8MYY4zfLGkYY4zxmyUNQERiRWSuiCwSkWUi8qjPc78TkVWufLRP+QMissY9d14oxCsi3UVkjogsFJFMEekTCvG65UeKyAIR+dw9ThKRKSLyo7tP9Jk31GJ9QkRWishiEZkgIgmhGqtP+T0ioiLSMJRjDcX1q6RYQ3zdyhaRJUWxubLArV+qWuNvgAB13XQ08APQDzgTmArUcs+luPuOwCKgFtASWAtEhkC8XwHnu/ILgBmhEK+L4Q/AeOBz93g0cL+bvh/4ewjHei4Q5ab/HsqxurJmwGS8C18bhmqsobp+lRJrKK9b2UXfs09ZwNYv29MA1HPAPYx2NwXuAB5X1Tw3X9FogkOBd1U1T1XXA2vwRh8MdrwK1HPl9YEtoRCviDQFLgRe9SkeCrzppt8ELvUpD6lYVfUrVc13D+cATUM1VucZ4D6830ORUIw1JNevUmINyXWrDAFbvyxpOG53dCHeMLNTVPUHoC1wuoj8ICIzRaS3m70JsNGn+iZXFux47waeEJGNwJPAA272YMf7D7yNWKFPWaqqbgVw9ymuPBRj9XUzMMlNh1ysInIJsFlVFxWbN+RiJXTXr3/wy1jvJjTXLfAS2lcikiUiI1xZwNYvSxqOqhaoane8f5F9RKQz3iBViXhNP/fijSwoeM1Dv3iJqooVSo33DuD3qtoM+D3wmps9aPGKyEVAjqpm+VulhLKQiFVEHgTygXFFRSXMFrRYRSQOeBB4qKQqJZQF+3MNufWrjFhDbt3yMUBVewLnA6NEZGAZ8550vMEY7jWkqeoeEZkBDMHLwh+r1xg4V0QK8Tr/2oTXblykKcd3V6tUsXiHA3e5pz7g+O51MOMdAFwiIhcAsUA9EXkb2C4iaaq6VUTS8PaYQjJWVb1eRIYDFwGD3e8h5GIF3sJrp17kbXtpCsx3B21DKlb3GwjF9au0WC8m9NYtAFR1i7vPEZEJeM1NgVu/qvKATajegGQgwU3XBr7B20CMBP7XlbfF260ToBM/P5i0jqo9EF5avCuAQa58MJDlpoMar0/cgzh+YPEJfn6gbnQIxzoEWA4kF5sn5GItVp7N8QPhIRdrqK5fpcQakusWUAeI95n+3v1eA7Z+2Z6GJw14U0Qi8Zrs3lfVz0UkBnhdRJYCR4Hh6n3yy0TkfbwNST4wSlULQiDePcA/RSQKOAKMAFDVYMdbksfxmiNuAX4CroSQjfU5vJVsivsHP0dVR4ZorCUK0VhfJzTXr5LcRmiuW6nABPe7jALGq+qXIjKPAK1f1o2IMcYYv9mBcGOMMX6zpGGMMcZvljSMMcb4zZKGMcYYv1nSMMYY47eQTxoiUuB6b1wmXq+ufxCRCPdchog8G+DlXyoiHQO5jJMhIiNF5NcnUT/dnfJY0nNPuM/9iYpH6FcMCSLymwrUe0RE7glETIEkIjeKSONSnpshIhnlfL0DfswzVkSuKOfrfu/HPNni05OuT/kgETm1jHqXikhJV66Xtaypvr21ljHf6e53u1BEapdnGeGmrO9VvJ6OV4rIUrftrPB2wlfIJw3gsKp2V9VOwDl4PUw+DKCqmap6Z4CXfylez5B+c+dyl/u5ilDVl1T13wFazu1AT1W915+ZT+J9JwDlThqhRDz+rk83AiUmjVCiqqVu9P0wCCir/n3AC+V8zbfw73dyHfCk224cPtHM7nqnakVERuJtL/uoamdgICV3IVJ+VX2VZQWueDxQ7HErYKf7AAZx/IrNPnhXQy5w9+1c+Y3AJ8BnwHrgt3jdHi/A67E0yc3XGvgSyMK7wro93o9+l6u30M3zi/lc/bHA08B04KliMd+I1/XAZ8DXeFduvg7Mc3EMdfPVBt4FFgPv4XV5nlH8cwCuAMa66UeAe9z0DOAxYCbwR6CXm87C6yo7zc3XC++q0Nl4V44uLeFznwgUuPd9FdACmOZimwY0r8D7ruvqzgeW+Lzvd4HDbllPuLJ73eezGHjU5zUfBFbhdan9TtF7L7bcscBL7vtZDVzkymOBN9yyFwBnuvIvgK5uegHwkJv+K3BrafEA6XhXCr/g6rUoFsdDrs5SvHGaxX13B9x7WAjULlZnBl7363Nd7Kf7fJbP+cz3OcevUD4APOU+12kUu3Ld5zN5Fm/dWAdc4fNcaZ/1AXcf4d7jMrfcL4rq4111/qjPd9refS7bgM3uPZ5eLJa2wPRisb2I9xtaB5yBt36swP3O3XyJlPBbLfbat3J8nR3nPvMn3HewBLjKzTfILW88sLyE13kRyHTv2fczeRzvwrjFeIkJvAvnluKtU7P8/L7+jrdeTsXbds1w7/0SN0+ki7voe7ndlQvexabLgf/4fhfF4v8JaF3KZ5TN8d/YXOCUcm2TK2PDHsgbxZKGK9uNdyXkII4njXocH/PgbOAjny9vDRCP1/3GXmCke+4Z4G43PQ1o46b7Al/7/KB9V7Cy5vucEi7JdzFs4niCegy43k0n4G0c6uAls9ddeVe8KzbLmzRecNPReBuIZPf4Kp/XXgyc4aZLTBolLPMzvCt2wevp9ZMKvO8ooJ6bbui+F8HbyCz1qXcuxzeyEe71B+IluyVAnPu+11B60vjS1W3jYojFS6RvuHna461YsXjdLIxyrzkPmOzmmQ60KyOedLyeUPuV8vkl+Uy/BVzs8z1llFJnBi754u1VT/X5LEvbCClwnZt+yHe+Yp/JBy7+jsCasj5r3+8f7/f2hXu+Ed7655s0fuemfwO8Wvx3WUIsN+HzB8PF9q6LYSiwD+jilpcFdPeZ90egwQm2GWN94rscmIK3EU5133ka3rbjINCyrO/O1ZuBtz4m4SX7oouiE9z9EqBJsbITfV9FY3NMwBurIxroBix05SOAv7jpWngJrCXwK5/30xjYQ7Gkgbet213G55MNPOimf00JXdCUdQvXbkRK2s2qj9e1Rhu8LyXa57npqrof2C8ie/E2gOB92V1FpC7eXsUH7nJ88L6ony/0xPN9oKVfkj9FVXe56XPxOkUrao+PBZrjbYieBVDVxSKyuJTXKst77r4d0JnjXV9EAltFpD7eD3umm+8tvN4xT6Q/3g+2qM5on+f8fd8CPOZ64SzE65I5tYQ657rbAve4Lt7GPx6YoKqHAERkYhnxvq+qhcCPIrIOL0mcBvwLQFVXisgGvH+93wB34v07/Q9wjng9xqar6ioRua2UeH4CNqjqnFJiOFNE7sNLckl4/1o/K2VeXx+7+yy8xHQihRz/3t/2qV/cJ+4zWS4iRZ97aZ/1LJ96p+F9x4XANhGZXka8v+LE0oDcYmWfqaqKyBJgu6ouARCRZXifwUI3Xw7exnKnH8spiv0d9/vcLiIzgd54iWmuemNKlGSYeN2MR7l4O+L9uz8CvCoi/8FLBADfAWNd9xylffa+juL9qQFvG5Snqsfce0935efibZuKjlfUx/teBvq8ny0i8nUJry+cuOfad3zun/Ej5v8Ku6QhIq3wmk1ygA4+T/0VLzlcJiLpeP8OiuT5TBf6PC7E+wwigD3qdTVelhPNd7CMur7PCXC5qq7yncFt3Ev7sn3LY/1YjgDLVLV/sWUklLGM8vB9DX/f93V4e3u93EqSTcnvRYC/qerLPysUuRv/Yy8+n1J6m+48IAOveWAK3l7QbXgbwbLiSaeU9y4isXhNOhmqulFEHqHs781X0e+zgOPraD4/PwZZ1muV9hn5rgfic/+L91bMidrCS4q3LIfxNoIlvYbv+ln02Pc1Y119f5UVe2nfXUvgHqC3qu4WkbFArKrmi9dj8GDgarym7rNUdaSI9MUbuGmhiHSn7O/rmLq/+fi8X1Ut9Dn2J3h7cJOLxXYBJ1gHVHWfiBwUkVaquq602UqZPqFwOBD+XyKSjNdW/ZzPh16kPl4bKni7hn5T1X3AehG50i1HRKSbe3o/3j/cE81XHpOB34nLEiLSw5XPwtuwIt74GF196mwXkQ7uYOtlfixjFZAsIv3d60WLSCdV3QPsFZHT3HzX+Rnz93grSlGdb/2s56s+3lgFx0TkTLzjJODzGTuTgZvdnh0i0kREUvA+n8tEpLaIxON1V12aK0UkQkRa4x0HW8XPP9+2eHt3q1T1KF4Pq8PwjnN9g7fR+OYE8ZSlaCOxw9XzPcOl+Pv1RzbQ3b2nZvx8tLUIn9e/lvJ9N/68t2+By92yi5qFT6Ss97gCOKUcMeJiE7zmsWz3eJqInGgAoVnAVeINWpaM90997gnq1MNLKHvd+z3fLa8uUF9Vv8AblKm7K2+tqj+o6kPADryux7Mp/fvyx2TgDhGJdstoKyJ13Pu52r2fNLwhc0vyN+B5Eann6teT4wM0gddcXXQ/uzyBhcOeRm3xRqiLxsveb+EdeC1uNF7z1B/wDrqW13XAiyLyF7esd/EObL0LvCIid+KtmKXNVx5/xRsdbLFbEbLxujZ/EXjDNUst5Oc/7vvxdoc34h10q1vWAlT1qNu1fdY1SUW5ZS7Da1N+XUQO4f04/XGnq3MvXtPCTX7W8zUO+ExEMvHe30oX604R+U68U38nqeq9ItIBmO3y6gG8Y0DzReQ9V3cDxzfqJVmFdxJAKt4xrCMi8gLwkmsGyAduVDfUqHutwap6SES+wRtn4BsX31clxYP3z7pE6o1z8gpe80M23t5MkbEujsNAf/XjDB+8JpD17vWW4h14LnIQ6CQiWXjH7K76ZfVS4yztveX4zPYR3r/rpXjH335wyynLZ8CHIjIU7x+z73c1C3hKRKSEP39l6YXXw3C++/N0Ct5B77JMwGtaXYT3j/o+Vd0mIu1Lq6Cqi0RkAd66sg7vswcvCX7q9iIFbzAm8Eb0a+PKpnF8e1Da9+WPV/Gaqua7bUQu3pmcE4Cz3OuuxvuNl+RFvG3EPBE5BhzDO1miSC0R+QHvD8c15QnMerkNYeINrnSPqmYGO5Zw4poTPlfVD4MdS3UhInVV9YCINMD7MzNAVbedxOv9E+84xtRy1pmoqtPcnvjNqvqHisZQU4nXLJyhqjsqUj8c9jSMMcH3uXjHw2KAv55MwnAewzv7sDyWquo0AFVdine2oalitqdhjDHGb2F1INwYY0xwWdIwxhjjN0saxhhj/GZJwxhjjN8saRhjjPHb/we7McnBi7pgmgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "## influence of hub height on diameter for fixed power\n",
    "cp_assume = 0.46 # assuem a likely Cp max\n",
    "Prated=10000000.0 # rated power in watts\n",
    "Pturb=Prated/(0.965*0.985) # power needed\n",
    "Hi = np.zeros(50)\n",
    "V = np.zeros(50)\n",
    "Di = np.zeros(50)\n",
    "\n",
    "for ic in range(50):\n",
    "   Hi[ic]=50.0+ic*(250.0-50.0)/50.0\n",
    "   A=0.14* np.log(Hi[ic]/10.0)\n",
    "   V[ic] = Uten*np.exp(A)\n",
    "   Di[ic] = np.sqrt(Pturb/(cp_assume*0.125*Density*np.power(V[ic],3)*np.pi))\n",
    "\n",
    "plt.figure()\n",
    "plt.xlabel('Wind Speed at Hub Height (m/s)')\n",
    "plt.ylabel('Hub height above sea level (m) ')\n",
    "plt.plot(V,Hi)\n",
    "\n",
    "plt.figure()\n",
    "plt.ylabel('Hub Height (m)')\n",
    "plt.xlabel('Diameter required for rated power at hub height (m), for assumed Cp ')\n",
    "plt.plot(Di,Hi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Need to input in the next box the input values for the selected turbine design"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Current Diameter  362.5  [m]\n",
      "Number of Blades 3 \n",
      "Current tip chord (m) 10\n",
      "Current pitch angle[deg] -1\n",
      "Current design [rpm] 3.1\n",
      "Current Hub Height 206.8  [m]\n",
      "Current rated speed  7.061761459072345  m/s\n",
      "Current clearance 25.55000000000001 [m]\n",
      "Current Rated Power 10.0  [MW]\n"
     ]
    }
   ],
   "source": [
    "print('Current Diameter ',D,' [m]')\n",
    "print('Number of Blades',Nblades,'')\n",
    "print('Current tip chord (m)',actual_chord)\n",
    "print('Current pitch angle[deg]',pitch)\n",
    "print('Current design [rpm]',RPM)\n",
    "print('Current Hub Height',HH,' [m]')\n",
    "print('Current rated speed ',vrated,' m/s')\n",
    "print('Current clearance', clearance,'[m]')\n",
    "print('Current Rated Power', Prated/1000000.0,' [MW]')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  Wind :   Pitch :  rpm:  TSR:  Cp:  Ct  Power (MW): \n",
      "  speed    [deg]          coeff     coeff    [N]     power \n",
      "  [m/s]\n",
      "   0.510  -1.000     3.1  115.325  -0.0000  0.0000  -0.0000\n",
      "   1.020  -1.000     3.1  57.663  -0.8617  0.0704  -0.0578\n",
      "   1.531  -1.000     3.1  38.442  -0.3395  0.1788  -0.0768\n",
      "   2.041  -1.000     3.1  28.831  -0.1155  0.3216  -0.0619\n",
      "   2.551  -1.000     3.1  23.065  0.0839  0.4033  0.0880\n",
      "   3.061  -1.000     3.1  19.221  0.1680  0.4535  0.3042\n",
      "   3.571  -1.000     3.1  16.475  0.2171  0.4995  0.6242\n",
      "   4.082  -1.000     3.1  14.416  0.2887  0.6498  1.2387\n",
      "   4.592  -1.000     3.1  12.814  0.3678  0.7684  2.2475\n",
      "   5.102  -1.000     3.1  11.533  0.4146  0.7791  3.4751\n",
      "   5.612  -1.000     3.1  10.484  0.4375  0.7695  4.8806\n",
      "   6.122  -1.000     3.1   9.610  0.4488  0.7535  6.5002\n",
      "   6.633  -1.000     3.1   8.871  0.4533  0.7357  8.3476\n",
      "   7.143  -1.000     3.1   8.238  0.4536  0.7190  10.4337\n",
      "   7.653  -1.000     3.1   7.688  0.4463  0.6949  12.6255\n",
      "   8.163  -1.000     3.1   7.208  0.4439  0.6817  15.2405\n",
      "   8.673  -1.000     3.1   6.784  0.4407  0.6682  18.1468\n",
      "   9.184  -1.000     3.1   6.407  0.4363  0.6541  21.3268\n",
      "   9.694  -1.000     3.1   6.070  0.4106  0.6029  23.6071\n",
      "  10.204  -1.000     3.1   5.766  0.4055  0.5888  27.1919\n",
      "  10.714  -1.000     3.1   5.492  0.3862  0.5455  29.9749\n",
      "  11.224  -1.000     3.1   5.242  0.3636  0.5124  32.4489\n",
      "  11.735  -1.000     3.1   5.014  0.3468  0.4757  35.3703\n",
      "  12.245  -1.000     3.1   4.805  0.3118  0.4274  36.1321\n",
      "  12.755  -1.000     3.1   4.613  0.3044  0.4115  39.8614\n",
      "  13.265  -1.000     3.1   4.436  0.2733  0.3596  40.2583\n",
      "  13.776  -1.000     3.1   4.271  0.2426  0.3183  40.0267\n",
      "  14.286  -1.000     3.1   4.119  0.2105  0.2700  38.7393\n",
      "  14.796  -1.000     3.1   3.977  0.1955  0.2479  39.9557\n",
      "  15.306  -1.000     3.1   3.844  0.1581  0.1983  35.7774\n",
      "  15.816  -1.000     3.1   3.720  0.1367  0.1695  34.1303\n",
      "  16.327  -1.000     3.1   3.604  0.1197  0.1456  32.8794\n",
      "  16.837  -1.000     3.1   3.495  0.0879  0.1063  26.4893\n",
      "  17.347  -1.000     3.1   3.392  0.0737  0.0895  24.2938\n",
      "  17.857  -1.000     3.1   3.295  0.0642  0.0774  23.0584\n",
      "  18.367  -1.000     3.1   3.203  0.0454  0.0543  17.7658\n",
      "  18.878  -1.000     3.1   3.117  0.0345  0.0417  14.6311\n",
      "  19.388  -1.000     3.1   3.035  0.0300  0.0364  13.7869\n",
      "  19.898  -1.000     3.1   2.957  0.0264  0.0322  13.1344\n",
      "  20.408  -1.000     3.1   2.883  0.0199  0.0249  10.6724\n",
      "  20.918  -1.000     3.1   2.813  0.0173  0.0220  10.0082\n",
      "  21.429  -1.000     3.1   2.746  0.0156  0.0199  9.6680\n",
      "  21.939  -1.000     3.1   2.682  0.0153  0.0196  10.2257\n",
      "  22.449  -1.000     3.1   2.621  0.0146  0.0187  10.4566\n",
      "  22.959  -1.000     3.1   2.563  0.0139  0.0177  10.5875\n",
      "  23.469  -1.000     3.1   2.507  0.0133  0.0173  10.8430\n",
      "  23.980  -1.000     3.1   2.454  0.0127  0.0163  11.0750\n",
      "  24.490  -1.000     3.1   2.403  0.0122  0.0156  11.3072\n",
      "  25.000  -1.000     3.1   2.354  0.0117  0.0150  11.5309\n",
      "  25.510  -1.000     3.1   2.307  0.0113  0.0144  11.7966\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAzjElEQVR4nO3dd3wUdf7H8dcnjZACoYSa0HsvkSLVjhW7AlZU5Kzn6d15zfNOz5/1LAcKqIgNET0LKgqHVKWGXgMBIQmhJEAaISHl+/tjBlniJiQxm9nyeT4eeWR3Z3bms5PNvne+35nviDEGpZRSqqwgpwtQSinlnTQglFJKuaUBoZRSyi0NCKWUUm5pQCillHJLA0IppZRbGhDKa4lIUxFZKiK5IvKSWN4RkWMislpEholIUiWWM05E5tdGzZ5S2dfqrUTkDhH5wcPr2CoiIz25jkCjAeHlRGSviJwQkTwROWR/QEY5XZcrEbnE5YM8Q0SWiMhVNbDoCUAmUM8Y8ygwFLgIiDPGDDDGLDPGdD7bQowxHxpjLq6BehARIyIdamJZVVlPZV9roBCRGSLytOtjxpjuxpjFDpXklzQgfMOVxpgooB9wDvBXJ4oQkWA3j10PfAK8B8QBTYEngCtrYJWtgW3m9NmcrYG9xpjjNbBs5aVEJMTpGpTNGKM/XvwD7AUudLn/AvC1ffsqYCuQBSwGutqP3wl85fKcZGC2y/1UoI99uwvwP+AokATc6DLfDOANYC5w3LUOe7oAKcDvK6g/CCvQ9gGHsYKkvsv0QcBy+zVsBEa6rLsIOAnkAfcCBUCJff8fwEggzWVZ8cBnQAZwBJhkP34H8IPLfGd7zZOBb4BcYBXQ3p62FDD2tsgDbirzWuvYr6OHy2OxwAmgCdAY+Nqe5yiwDAhys81+sR43r3Uv8CdgG3AMeAcIL+dv0AFYAmRj7ZF97DLNAA8Be+xpL7jWBIwHttvrmAe0ruR2bATMAXKA1cBTrn+DMvW1seu4C+v9tNR+/BPgoF33UqC7/fiEMu+Nr8r+r9h/i1eAdPvnFaCO0//PvvbjeAH6c5Y/0Jlv+nisQHgK6GR/gFwEhAJ/wAqCMKCd/SEUBDTH+nDeby+jnf3PHgREYoXFnUAI1h5Kpss/4gz7n3OIPX94mdq62P/YbSuof7xdVzsgCusD/H17WkusD/LL7OVfZN+PdVn/0y7LuoMzP+hHYn9oAsFYAfOy/brCgaFln1fJ13wUGGBP/xCY5bJOA3So4PVOB/7lcv9+4Dv79v8BU+y/VygwDJBylnPGenAfEFvs90RD4EfXbVVmWR8Bfzn1Nzy1XVzWs8heRitgJ3C3Pe1q+2/X1d4WfwWWV3I7zgJm2/P1APZz9oB4z56/rst7J5rTH/YbXJ5zxnvDzf/KP4GVWMEci/Ul5Cmn/5997cfxAvTnLH8g602fh/WBvw94HagL/I0z9wqC7H/Ckfb9VPuf9mZgGta3uC72P/Qce56bgGVl1jcV+Lt9ewbwXgW1DbH/sd1+c7Xn+R64z+V+Z6xvfyHAH7HDwmX6POB2l/VXNiAGY+05hLip4efnVfI1v+Uy7TJgh8v9swXEhcAel/s/ArfZt/8JfFnR88tbD+4DYmKZOneXs6z37PdAXDnrGeVy/z7ge/v2t8BdZd5j+VhNfeVuR6ywLgK6uEx7hrMHRLsKtkeMPU99d+8Nl21yKiB2A5e5TLsEq3nS8f9pX/rRPgjfcLUxJsYY09oYc58x5gTQAiswADDGlGKFQkv7oSVYHyrD7duLgRH2zxJ7ntbAQBHJOvUDjAOauaw7tYK6jti/m1cwzxl12rdDsPoqWgM3lFn/0LMsrzzxwD5jTPFZ5qvMaz7ocjsfa8+nshYCdUVkoIi0BvoAn9vTXsD6Rj5fRPaIyONVWK47rn+bfVjb2p0/YDUHrraP9BlfyeW0Bl512U5H7eW0pOLtGIv1Ny673Eq/HhEJFpFnRWS3iORgffiD1UxXGe7ed+VtH1UO7QzyXelAz1N3RESwPiT32w8tweoobov17S0L6x94MDDJnicVWGKMuaiC9ZgKpiXZy7gOeLGCOlu73G8FFAOH7Oe+b4y5p4J1VFYq0EpEQs4SEpV5zdVmjCkVkdnAGKzX+LUxJteelgs8CjwqIt2BRSKyxhjzfTVXF+9yuxXWtnZX00HgHgARGQosEJGlxphkl+VsdbOcVKzmsg/LLtMOP7fb0T6Yodhe7g6X5Z6N63ttLDAaa49sL1Afq2lU3Mzrzqn3nbvXpSpJ9yB812zgchG5QERCsT54CrHaWsEKiPOw2nPTsDpER2F1Hq635/ka6CQit4pIqP1zjoh0rUwBxtp3/x3wNxG5U0TqiUiQiAwVkWn2bB8Bj4hIW/vw3GewOkmLgQ+AK+3DZINFJFxERopIXDW2x2rgAPCsiETayxriZr5f9ZqxPvTbnWWemVhNMOPs2wCIyBUi0sEO8xysDveSX7Ge+0UkTkQaAn8GPnY3k4jc4LJNj2F9uLqu9/ci0kBE4oGHXZYzBfiTHWaISH0RucGeVu52NMaUYPU1PSkiESLSDbj9LK+lrGis9/MRIALrfePqbNvnI+CvIhIrIo2xjqz7oIo1BDwNCB9ljEkCbgH+g9U5eCXW4bAn7ek7sfoultn3c7COVPnR/gc+9Y32Yqx+inSsppXnsDoFK1vHp1gfhuPtZRwCnsZqawer0/Z9rKNQfsI6EulB+7mpWN8S/4zVf5AK/J5qvC/t13Ql1hE7KUCaXVfZ+X7ta34SeNduVrmxnFpWYR1A0AKrHf+UjsACrL/LCuB1U/5x+2ddD1b4zMf6u+7B2u7unAOsEpE8rCOLHjbG/OQy/UtgLbAB6+itt+3X8TnWtpllN/NsAS61p51tOz6A1TR3EKu/4J1yaivPe9gHV2AdqbWyzPS3gW729vnCzfOfBhKBTcBmYB3lbx9VDrE7cJRSPkRE9mIdbbTgVy7HAB1dmpuU+pnuQSillHLLowEhIqNEJElEkss7YsNuc95gH12xxN08Simlap/HmpjsIxl2Yp38lAasAcYYY7a5zBOD1ak6yhiTIiJNjDGHPVKQUkqpKvHkHsQAINkYs8fuOJ2F1SHpaizwmTEmBUDDQSmlvIcnz4NoyZknyqQBA8vM0wkIFZHFWIe1vWqMea/sgkRkAtb4K0RGRvbv0qWLRwpWSil/tXbt2kxjTGxVnuPJgBA3j5VtzwoB+gMXYA0fsUJEVtqHaJ5+kjHTsIYKICEhwSQmJnqgXKWU8l8iUpmz2c/gyYBI48wzPeP45ZmMaUCmsYZvPi4iS4HeWH0XSimlHOTJPog1QEf7DNowrBNq5pSZ50tgmIiEiEgEVhPUdg/WpJRSqpI8tgdhjCkWkQewRucMBqYbY7aKyER7+hRjzHYR+Q7rbMdSrFE0t3iqJqWUUpXnc2dSax+EUkpVnYisNcYkVOU5eia1UkoptzQglFJKuaUBoZRSyi0NCKWUUm5pQCillHJLA0IppZRbGhBKKaXc0oBQSinllgaEUkoptzQglFJKuaUBoZRSyi0NCKWUUm5pQCillHJLA0IppZRbGhBKKaXc0oBQSinllgaEUkoptzQglFJKuaUBoZRSyi0NCKWUUm5pQCillHJLA0IppZRbGhBKKaXc0oBQSinllgaEUkoptzQglFJKuaUBoZRSyi2PBoSIjBKRJBFJFpHH3UwfKSLZIrLB/nnCk/UopZSqvBBPLVhEgoHJwEVAGrBGROYYY7aVmXWZMeYKT9WhlFKqejy5BzEASDbG7DHGnARmAaM9uD6llFI1yJMB0RJIdbmfZj9W1mAR2Sgi34pIdw/Wo5RSqgo81sQEiJvHTJn764DWxpg8EbkM+ALo+IsFiUwAJgC0atWqhstUSinljif3INKAeJf7cUC66wzGmBxjTJ59ey4QKiKNyy7IGDPNGJNgjEmIjY31YMlKKaVO8WRArAE6ikhbEQkDbgbmuM4gIs1EROzbA+x6jniwJqWUUpXksSYmY0yxiDwAzAOCgenGmK0iMtGePgW4HviNiBQDJ4CbjTFlm6GUUko5QHzt8zghIcEkJiY6XYZSSvkUEVlrjEmoynP0TGqllFJuaUAopZRySwNCKaWUWxoQSiml3NKAUEop5ZYGhFJKKbc0IJRSSrmlAaGUUsotDQillFJuaUAopZRyy5PDfSsFQHZ+Ect3Z1JcaggJEoKC5IzfwUFCsAiNosJo3SiS0GD93qKUN9CAUB5xsriURUmH+XzdfhbuOMzJktJKPS80WGjbOJKOTaLp2DTq599tGkUSFqLBoVRt0oBQNcYYw7qUY3y2bj/fbD5AVn4RjaPCGDeoFVf0ak79uqEUlxpK7J8zbpcYDucWsOtwHrsO5bElPZu5Ww5waizJkCChTeNIOjaJYkDbhlzeqzlNosOdfcFK+TkNCFVlxhiy8ovYdzSflKP5pBw5TsrRfFbuOUrK0XzCQ4O4pHszru7bkmEdGhNSzSajgqISdmdYgbHrcC67DuWx7UAO3245yFNfb+Pc9o25qk8LRvVoRr3w0Bp+lUopHe5bVai4pJRdh/PYlJbFxrRstuzP5qeM4+QWFp8xX2x0Hbo2r8dVva0P7Kg6nvvusetQLnM2pvPlhnRSjuYTFhLE+Z2bcFWfFpzfpQnhocEeW7dSvqo6w31rQKgzGGPYvD+b+VsPsXLPEbam53CiqASA6Doh9IyrT4cmUbRqGEHrRpG0ahhBfMO6RITV/s6oMYaNadl8uWE/X208QGZeIVF1QrikezNG92nBue0bVXvvRSl/owGhqqW4pJTVPx1l3taDzN92iAPZBQQHCb3j6tM7PobecTH0iqtPm0aRBAWJ0+W6VVJqWLH7CHM27ufbLQfJLSimcVQYl/dszt3D2hHfMMLpEpVylAaEqpa7Zqzh+x2HqRMSxPBOsVzSvRkXdGlCg8gwp0urloKiEhYnZfDVxnT+t/0QAvxmZHsmjmivzU8qYFUnILSTOsAl7j3K9zsOM3FEex66oIMjTUU1LTw0mFE9mjGqRzPSs07wr7nbeWXBLj5dm8ZfL+/GJd2bIuKde0JKeRNtoA1wkxYl0zAyzG/CoawWMXWZPLYfM+8ZSGRYCBM/WMtt01eTfDjP6dKU8noaEAFsU1oWi5MyuHtYW78MB1fntm/MNw8N5e9XdmNDahajXlnKM3O3k1tQ5HRpSnktDYgANmlhMvXCQ7h1UGunS6kVIcFB3DmkLYseG8m1/Voybekezn9pCZ+tS8PX+uKUqg0aEAFqx8Ec5m87xJ1D2hIdYCeZNY6qw/PX9+aL+4fQon44v5u9keunrGDL/mynS1PKq2hABKhJC5OJqhPCnUPaOF2KY/rEx/D5fUN4/rpe7M08zpWTfuAvn2/m2PGTTpemlFfQgAhAuzPy+GbzAW4d3JqYCN88lLWmBAUJN54Tz8LHRnL74DbMWpPKeS8t5v2V+ygp1WYnFdg0IALQ5EXJ1AkJ4q6hbZ0uxWvUrxvKk1d1Z+5Dw+jSLJq/fbGF3368geJKjkKrlD/SgAggBUUlfLQ6hS83pDNuYGsaR9VxuiSv07lZNB/dM4g/jOrMVxvTefjjDRRpSKgA5d/HNioAsvJP8sHKfcxYvo/MvEJ6tqzPxBHtnS7La4kI943sQEiQ8MzcHZSWGl4b01cvZKQCjkcDQkRGAa8CwcBbxphny5nvHGAlcJMx5lNP1hRo1u47xq1vryL/ZAkjOsUyYXg7zm3fSM8kroQJw9sTJMLT32zngZnr+M+YfnrRIhVQPBYQIhIMTAYuAtKANSIyxxizzc18zwHzPFVLoDLG8K9vtlEvPJTP7juXLs3qOV2Sz7l7WDuCg4R/fLWN+2euY/JYDQkVODz5Th8AJBtj9hhjTgKzgNFu5nsQ+C9w2IO1BKSFOw6zLiWLhy/sqOHwK9w5pC3/HN2d/207xG8+WEthcYnTJSlVKzwZEC2BVJf7afZjPxORlsA1wJSKFiQiE0QkUUQSMzIyarxQf1RaanhhXhJtGkVwff84p8vxebcNbsPTV/ewBjZ8fy0FRRoSyv95MiDcNXKXPbD8FeCPxpgK/9uMMdOMMQnGmITY2Niaqs+vfb35ADsO5vLIRZ20c7WG3DKoNc9c05NFSRncqyGhAoAnPznSgHiX+3FAepl5EoBZIrIXuB54XUSu9mBNAaGopJR/z0+iS7NoruzVwuly/MrYga147rqeLN2VwT3vJWpIKL/myYBYA3QUkbYiEgbcDMxxncEY09YY08YY0wb4FLjPGPOFB2sKCP9dm8beI/k8dnFnr70CnC+76ZxWPH9dL35IzmT8jDWcOKkhofyTxwLCGFMMPIB1dNJ2YLYxZquITBSRiZ5ab6ArKinl1e930bdVDBd0beJ0OX7rhoR4Xry+Nyv2HOHOGavJP1nsdElK1TiPngdhjJkLzC3zmNsOaWPMHZ6sJVDM32pdU/pf1/TQcx087Lr+cQQHCb+bvYE73lnDO3ecQ2QdPfdU+Q/tvfQzH67aR8uYuozopHsPteHqvi155ea+rN13jDveWU1eoe5JKP+hAeFHdmfksXz3EcYObEWw9j3Umqt6t+C1m/uyLiWL26ev1qvUKb+hAeFHPlyZQmiwcGNC/NlnVjXq8l7NmTSmLxtTs7ht+mpyNCSUH9CA8BMFRSV8ujaVS7o3IzZaR2l1wqU9mzNpbD82p2Vz69uryT6hIaF8mwaEn/hqYzo5BcWMGxgY15f2VqN6NOONW/qzLT2bW99eRXa+hoTyXRoQfuLDVSm0j41kULuGTpcS8C7q1pQpt/Rnx4Fc7p+5Tq9Mp3yWBoQf2LI/mw2pWYwb2FoPbfUSF3RtylNXd+eH5ExeX5TsdDlKVYsGhB+YvCiZqDohXKeD8nmVGxPiubpPC15esJOVe444XY5SVaYB4eO2pefw7ZaDjB/Shvp1Q50uR7kQEZ6+pidtGkXy0EfrycwrdLokpapEA8LHvfr9TqLDQ7hraDunS1FuRNUJYfK4fmSfKOKRjzdQqv0RyodoQPiwrenZzNt6iPFD2lI/QvcevFXX5vX4+5XdWbYrkzeW7Ha6HKUqTQPCh726YBfR4SGMH9rW6VLUWYwZEM9VvVvw0vwkVv901OlylKoUDQgftTU9m/nbDnHX0Lba9+ADRIRnru1J60aRPPjROo5of4TyARoQPmreloMECdx5ru49+IqoOiFMGtuXY/lF/G72Ru2PUF5PA8JHrUvJokuzetr34GO6t6jPE1d0Y8nODKYs1f4I5d00IHxQSalhQ2oWfVvFOF2KqoZxA1txea/mvDR/J2v2an+E8l4aED4o+XAeeYXF9GvVwOlSVDWICM9e25O4BnV56KP1HD1+0umSlHJLA8IHrU85BqB7ED4sOjyUyWP7cSTvJI/O1vMjlHfSgPBB61KOERMRStvGkU6Xon6FHi3r89crurIoKYM3l+1xuhylfkEDwgetT8mib3yMDsznB24d1JrLejbj+XlJrN2n/RHKu2hA+JjsE0XsOpxHX+1/8AsiwrPX9aJlTF0enLmeY9ofobyIBoSP2ZiaBaAd1H6kXngok8b2JSOvkMc+2Ygx2h+hvIMGhI9Zn5KFCPSOr+90KaoG9YqL4S+XdeX7HYd5a9lPTpejFKAB4XPWpRyjU5NoosP1BDl/c/u5bRjVvRnPfbeDdfaRako5SQPCh5TqCXJ+TUR47vpeNKsfzoMz15OVr/0RylkaED5kWXIm2SeKNCD8WP26oUwa24/DuQU89skm7Y9Qjqp0QIjItSLybxF5SUSu8WRR6pfWpxzjvg/W0rFJFKN6NHe6HOVBfeJjePzSrizYfojpP+51uhwVwCoVECLyOjAR2AxsAe4VkcmVeN4oEUkSkWQRedzN9NEisklENohIoogMreoLCATb0nO4ffpqGkfX4YO7B+rw3gFg/JA2XNStKc9+u50N9pFrStU2qcwurIhsBXoYe2YRCQI2G2O6V/CcYGAncBGQBqwBxhhjtrnMEwUcN8YYEekFzDbGdKmoloSEBJOYmHj2V+ZHrn39R/ZnneDTiecS3zDC6XJULcnOL+Ky15YhAt88OExH7lW/ioisNcYkVOU5lW1iSgJaudyPBzad5TkDgGRjzB5jzElgFjDadQZjTJ45nVCRgDa4llFQVMLm/dlc0zdOwyHA1I+wzo84mF3A7z/V8yNU7atsQDQCtovIYhFZDGwDYkVkjojMKec5LYFUl/tp9mNnEJFrRGQH8A0w3t2CRGSC3QSVmJGRUcmS/cP2AzkUlRj66HkPAalvqwY8fmkX5m87xIzle50uRwWYkIomikgHoCnwRJlJI4D9QHJFT3fz2C++AhljPgc+F5HhwFPAhW7mmQZMA6uJqaKa/c2pM6d7x8c4Wodyzl1D27JyzxGembud/q0b0CsuxumSVIA42x7EK0CuMWaJ6w8wF7ja5b47aVhNUafEAenlrcgYsxRoLyKNK1++/9uYlk2T6Do0qxfudCnKISLCizf0pkl0OPfPXEf2iSKnS1IB4mwB0cYY84u+BmNMItDmLM9dA3QUkbYiEgbcDJzRHCUiHcQeklRE+gFhwJFK1h4QNqZm0VtHbg14MRFhvDamLweyCvjDpxsp0etHqFpwtoCo6Gtr3YqeaIwpBh4A5gHbsY5Q2ioiE0Vkoj3bdcAWEdkATAZuMtoT97Ps/CL2ZB6njzYvKaB/6wb86bKuzNt6iIdnredkcanTJSk/V2EfBLBGRO4xxrzp+qCI3AWsPdvCjTFzsZqjXB+b4nL7OeC5ypcbWDbtzwKgt7Y5K9tdQ9tSUlrKM3N3cLywmDdu6U94aLDTZSk/dbaA+C1WB/I4TgdCAlZTkJ5N7WGnOqh7xukRTOq0CcPbE1UnlL98sZnbpq/m7dsTdPBG5REVNjEZYw4ZY84F/gHstX/+YYwZbIw56PnyAtuG1GzaxUbqmdPqF8YObMWrN/dl3b5jjHtrFUf1QkPKA862BwGAMWYRsMjDtSgXxhg2pmUxrIMe1KXcu6p3C6LqBPObD9Zx09QVfHD3QJrq0W6qBulorl7qYE4BGbmF9NLmJVWB87s0ZcadA0jPOsH1U5aTciTf6ZKUH9GA8FKb0rIB6KVHMKmzGNy+ETPvGURuQTHXT1nOzkO5Tpek/IQGhJfauj+b4CChW/N6TpeifEDv+Bg+njAYgJumrmBTWpazBSm/oAHhpTbvz6ZDbJQewqgqrXOzaD6ZOJjIOiGMfXMVK/foOafq19GA8ELGGDbvz6FHS+1/UFXTulEkn048l2b1w7l9+moW7TjsdEnKh2lAeKFDOYVk5hXSs6U2L6mqa1Y/nNn3DqZT02jueS+RrzaWOwSaUhXSgPBCW/ZbHdR6gpyqroaRYXx4z0D6tWrAQ7PWM2t1itMlKR+kAeGFNu/PJkigq3ZQq1+hXngo744fwIhOsTz+2WbeXLrH6ZKUj9GA8EJb9mfTPjaKiLBKnceoVLnqhgUz7dYELu/ZnH/N3c6/5yfplelUpeknkBfakp7NkPZ6BrWqGWEhQbw2pi9RdUJ4bWEyOQXFPHFFN4KCdAh5VTENCC9zOLeAQzmFdNcjmFQNCg4Snr2uJ1HhIbz9w0/kFRbz7LU9CQnWRgRVPg0IL7N1fw4APTUgVA0TEf56eVfqhYfy8oKd5BUU8+qYPtQJ0XNtlHv69cHLbDtgBUTX5tEOV6L8kYjw8IUdeeKKbny39SB3v5tI/slip8tSXkoDwsskHcylZUxdHd9fedT4oW15/rpe/JicyW1vr9brXCu3NCC8TNLBXLo0070H5Xk3nhPPpLH92JiWxZhpK8nMK3S6JOVlNCC8yMniUnZn5NFZA0LVkst6NufN2xLYk5nHjVNXkJ51wumSlBfRgPAiezLzKC41GhCqVo3s3IT3xg8kI6eQG6as4KfM406XpLyEBoQXSTpojePfpZmeQa1q14C2DflowiBOFJVww5QVbLcPllCBTQPCi+w4mEtosNAuNtLpUlQA6tGyPrPvHURIkHDT1BWsSznmdEnKYRoQXiTpYC7tY6MI1ZOXlEM6NLGuKdEgMoxb3lrF8uRMp0tSDtJPIi+y40AOnZpq/4NyVnzDCD65dzDxDSK4Y8Ya3l+xl8LiEqfLUg7QgPAS2SeKSM8u0A5q5RWa1Avn43sH0Tc+hr99uZWRLyxmxo8/UVCkQRFINCC8xKkLzes5EMpbxESEMWvCIN4bP4D4BhE8+dU2hj63kClLdpNXqGdfBwIdi8lLbEjJAvQaEMq7iAjDO8UyvFMsq/YcYdKiZJ79dgdvLN7N+CFtuePcNtSP0LP+/ZVH9yBEZJSIJIlIsog87mb6OBHZZP8sF5HenqzHm329KZ3uLerRIqau06Uo5dbAdo14/66BfHH/EM5p05CXF+xkyHMLee67HXoWtp/yWECISDAwGbgU6AaMEZFuZWb7CRhhjOkFPAVM81Q93mxv5nE2pmVzVe8WTpei1Fn1iY/hrdsTmPvQMEZ0jmXKkt0MfW4h//xqGwezC5wuT9UgTzYxDQCSjTF7AERkFjAa2HZqBmPMcpf5VwJxHqzHa526qPwVGhDKh3RrUY/JY/uRfDiP1xcn8+6KvXywch83JMQxcUR74htGOF2i+pU82cTUEkh1uZ9mP1aeu4Bv3U0QkQkikigiiRkZGTVYovOMMczZmM45bRrQUpuXlA/q0CSKf9/Yh8WPjeT6hDg+SUxj5IuLeeyTjezJyHO6PPUreDIg3F3P0O3FcEXkPKyA+KO76caYacaYBGNMQmxsbA2W6LykQ7nsOpynzUvK58U3jOCZa3qy5A8juW1wa77elM4F/17CAzPXseOgDt3hizwZEGlAvMv9OCC97Ewi0gt4CxhtjDniwXq80ua0bACGdNBrUCv/0Lx+Xf5+ZXeW/eF87h3enkU7DjPqlWXc814iG1OznC5PVYEnA2IN0FFE2opIGHAzMMd1BhFpBXwG3GqM2enBWrxW6rETiEBcA22vVf4lNroOj1/ahR8fP5/fXtiRVXuOMHryj9w2fTVr9h51ujxVCR4LCGNMMfAAMA/YDsw2xmwVkYkiMtGe7QmgEfC6iGwQkURP1eOtUo/m07xeOGEhes6i8k8xEWH89sJO/Pj4+fxxVBe27s/mhikruGnqCpbtysAYty3PyguIr/1xEhISTGKi/+TIdW8sJyRI+PjewU6XolStOHGyhI9WpzB16W4O5RTSOz6GB8/rwAVdmyDirutS1QQRWWuMSajKc/Rrq8NSj+bTSg8HVAGkblgw44e2ZekfzuNf1/TgSF4hd7+XyGWv/cA3mw5QUupbX1r9mQaEgwqKSjicW6jHi6uAVCckmHEDW7PosZG8dENvCotLuH/mOi5+eQmfrUujuKTU6RIDngaEg9KO5QPoHoQKaKHBQVzXP47/PTKCSWP7EhocxO9mb+S8lxYzc1WKDjXuIA0IB6UctQIivqGeIKdUcJBwRa8WzH1oGG/elkDDiDD+/PlmRjy/mHd+/IkTJzUoapsGhINSj54A0CYmpVwEBQkXdWvKF/cP4f27BtCqUQT/+Gobw57XocZrmwaEg1KO5hMeGkRsVB2nS1HK64gIwzrGMvvewXw8YRBdm9fj2W93MOTZhbz9w08UaR+Fx2lAOCj1aD7xDSL00D6lzsJ1qPFecfV56uttXP7aMpbv1mtme5IGhINS9BBXpaqkT3wM740fwNRb+5N/soSxb67i/pnrSM864XRpfkkDwiHGGNKOndD+B6WqSES4pHszFvxuBI9c2IkF2w5xwUtLmLwoWY94qmEaEA7JyCskr7BYA0KpagoPDebhCzuy4HcjGN6pMS/MS+Lil5eycMchp0vzGxoQDtmUao3i2rNlfYcrUcq3xTeMYOqtCbx/1wBCgoTxMxIZP2MNezOPO12az9OAcMiG1CyCg0QDQqkaMqxjLN8+PJw/X9aFVXuOcPHLS3lh3g7yT+phsdWlAeGQDalZdG4aTd2wYKdLUcpvhIUEMWF4exY+NpLLezVn8qLdXPDSEr7ZdEBHja0GDQgHlJYaNqZm0adVjNOlKOWXmtYL5+Wb+vDJxMHERIRx/8x1jH1zFTsP5Tpdmk/RgHDAnsw8cguL6RMf43QpSvm1c9o05OsHh/LU6O5sO5DDpa8u459fbSOnoMjp0nyCBoQD1qdkAdBXA0IpjwsOEm4d3IZFj43kxoR43ln+E+e/uJhPElMp1aHFK6QB4YCNaVlE1QmhfWyU06UoFTAaRobxf9f2ZM79Q4lvGMHvP93EdVOW/3xdePVLGhAOSNx7jD7xMQQF6RAbStW2nnH1+e/Ec3nxht6kHs3nqsk/8KfPNnP0+EmnS/M6GhC17Ojxk+w4mMugdg2dLkWpgBUUJFzfP46Fj41k/JC2zE5M5bwXF/P+ir16RTsXGhC1bNWeIwAMbt/I4UqUUvXCQ/nbFd349uFhdG9Rj799uZUr//MDa/Yedbo0r6ABUctW7jlC3dBgeraMcboUpZStU9NoPrx7IK+P60dW/klumLKCRz7ewOGcAqdLc5QGRC1bsecICW0aEBaim14pbyIiXNazOQseHcED53Xgm00HOO/FxUxbupuTxYF57Qn9lKpFmXmF7DyUp81LSnmxiLAQHrukM/MfGc6gdo14Zu4ORr26lKU7M5wurdZpQNSiDfb5DwPaaAe1Ut6uTeNI3r7jHKbfkUBpqeG26au59/1EUu1ryQcCDYhalGK/sdo2jnS4EqVUZZ3fpSnzHhnO7y/pzNKdmVz47yW8smAnBUX+f+0JDYhalHosn4iwYBpGhjldilKqCuqEBHP/eR34/tERXNStKa8s2MWF/17C/K0H/XoQQA2IWpR27ARxDerqNaiV8lEtYuoyaWw/Zt4zkIiwYCa8v5bb31nD7ow8p0vzCI8GhIiMEpEkEUkWkcfdTO8iIitEpFBEHvNkLd4g9Wg+8Q30CnJK+bpz2zfmm4eG8cQV3Vi/7xijXlnK/327nbxC/7r2hMcCQkSCgcnApUA3YIyIdCsz21HgIeBFT9XhLYwx7Lf3IJRSvi80OIjxQ9uy8LGRXN2nJVOX7OGClxbz5Yb9ftPs5Mk9iAFAsjFmjzHmJDALGO06gzHmsDFmDeD3Y+9mnygiV69BrZTfiY2uwws39Obz+86lab1wHp61gZumrmT7gRynS/vVPBkQLYFUl/tp9mMBKe3YCQDdg1DKT/Vt1YAv7hvCs9f2JDkjj8tfW8bfv9xCdr7vfv/1ZEC464mt1n6XiEwQkUQRSczI8M2TVU4dOx2nfRBK+a2gIOHmAa1Y9OhIbh3UmvdX7uO8lxYza3WKT157wpMBkQbEu9yPA9KrsyBjzDRjTIIxJiE2NrZGiqttp/YgtJNaKf9XPyKUf4zuwdcPDqNDbBSPf7aZq1//kfUpx5wurUo8GRBrgI4i0lZEwoCbgTkeXJ9XSz2WT3R4CPUjQp0uRSlVS7q1qMfH9w7i1Zv7cCingGteX87vP9lIZl6h06VVSoinFmyMKRaRB4B5QDAw3RizVUQm2tOniEgzIBGoB5SKyG+BbsYY3+/dKSP1aL42LykVgESE0X1ackHXpvxn4S6m//AT3209yCMXduK2wa0JCfbe09HE1w7HSkhIMImJiU6XUSUlpYb+T/+PC7o05aUbeztdjlLKQbsz8nhyzlaW7cqkc9Nonr6mB+fUwvhsIrLWGJNQled4b3T5kW3pOWTlFzGsY2OnS1FKOax9bBTvjR/A1Fv7c/xkMWOmreSrjdXqnvU4DYha8ENyJgDndtBhvpVSVrPTJd2bMffhYfRr1YCHZq3no9UpTpf1CxoQteCH5Ay6NIumSXS406UopbxIvfBQ3h0/gBGdYvnTZ5uZumS30yWdQQPCwwqKSliz9xhDOmjzklLql+qGBTPt1gSu6NWc//t2B89/t8Nrhurw2FFMyrJg+yFOFpdyfpcmTpeilPJSYSFBvHpzX6LDQ3h98W5yCor451U9CApyduRnDQgP+3RtGs3rhzOonfY/KKXKFxwkPHNNT+rVDWXqkj3kFRTzwg29CXXwMFgNCA86nFPA0p0Z/GZke4Id/iaglPJ+IsKfLu1K/bqhPP9dEnmFxUwa24/w0GBH6tE+CA/6butBSg1c2y/O6VKUUj7kvpEdeOrqHny/4zB3vLPasetMaEB40Nb9OTSMDKOdXoNaKVVFtw5qzcs39mHN3mOMe3Mlx46frPUaNCA8aMfBHLo2j9ZLjCqlquXqvi2Zekt/th/M5capKziYXVCr69eA8JCSUkPSoVy6NKvndClKKR92YbemvHvnANKzTnDD1OWkHMmvtXVrQHjIviPHKSgqpUuzaKdLUUr5uMHtGzHznkHkFhRz/ZTlJB3MrZX1akB4yPYD1h+wa3Pdg1BK/Xq942OYfe9gAG6cuqJWri2hAeEhOw7mEBwkdGgS5XQpSik/0alpNJ9OPJf6dUMZ99YqltvjvHmKBoSHrE/JomOTKMeOX1ZK+adWjSL4dOJg4htEcMeMNczfetBj69KA8IC8wmJW/XSE4Z188/KoSinv1qReOB/fO4iuzevxmw/X8dm6NI+sRwPCA37YlUlRieG8zjr+klLKM2Iiwvjw7oEMbNuQ383eyLvL99b4OjQgPGBx0mGiw0NIaNPA6VKUUn4sqk4I0+84h4u6NeXvc7YyaeGuGh0JVgOihu3NPM6XG9K5oEsTRwfZUkoFhvDQYN4Y149r+7bkxfk7eWbu9hoLCR2srwYVlZTy6CcbCQkW/nhpF6fLUUoFiJDgIF68oTfR4SG8uewnck4U88y1PX/1IKEaEDXo6a+3sXbfMV4b05fm9es6XY5SKoAEBQlPXtWdenVD+c/CZPIKi3n5pj6EhVS/JUMDoobszTzOeyv3cfvg1lzVu4XT5SilApCI8OjFnakXHsq/5m4nr7CYKbf0p25Y9Q6310byGvLByn0Ei3DfeR2cLkUpFeDuGd6O567rybJdGdz69iqyTxRVazm6B1ED1u47yszVKVzSvRlN64U7XY5SSnHTOa2IqhPKbz9ez5hpK6u1DA2IasotKOL77YdZl3KMWWtSaRlTl79d0c3pspRS6meX92pOZJ1gJn6wtlrP14CoJGMMyYfzWLorkyU7M1i55wgni0upGxrM+Z2b8PQ1PWgcVcfpMpVS6gwjOzfh/bsGMuDpqj9XA6KSJn6wlnlbDwHQLjaSWwa25vJezekbH0OQXm9aKeXFzmnTsFrP04CopGv7xXFe5yYM6dCY+IYRTpejlFIe59GAEJFRwKtAMPCWMebZMtPFnn4ZkA/cYYxZ58maquuS7s2cLkEppWqVxw5zFZFgYDJwKdANGCMiZXtxLwU62j8TgDc8VY9SSqmq8eR5EAOAZGPMHmPMSWAWMLrMPKOB94xlJRAjIs09WJNSSqlK8mQTU0sg1eV+GjCwEvO0BA64ziQiE7D2MAAKRWRLzZbqsxoDnr2klO/QbXGabovTdFuc1rmqT/BkQLg7tKfsEIOVmQdjzDRgGoCIJBpjEn59eb5Pt8Vpui1O021xmm6L00QksarP8WQTUxoQ73I/DkivxjxKKaUc4MmAWAN0FJG2IhIG3AzMKTPPHOA2sQwCso0xB8ouSCmlVO3zWBOTMaZYRB4A5mEd5jrdGLNVRCba06cAc7EOcU3GOsz1zkosepqHSvZFui1O021xmm6L03RbnFblbSE1eXk6pZRS/kOH+1ZKKeWWBoRSSim3fCogRGSUiCSJSLKIPO50PU4Skb0isllENlTn8DVfJiLTReSw6/kwItJQRP4nIrvs3w2crLG2lLMtnhSR/fZ7Y4OIXOZkjbVBROJFZJGIbBeRrSLysP14wL0vKtgWVX5f+EwfhD10x07gIqzDY9cAY4wx2xwtzCEishdIMMYE3ElAIjIcyMM6C7+H/djzwFFjzLP2l4cGxpg/OllnbShnWzwJ5BljXnSyttpkj8DQ3BizTkSigbXA1cAdBNj7ooJtcSNVfF/40h5EZYbuUAHAGLMUOFrm4dHAu/btd7H+IfxeOdsi4BhjDpwa6NMYkwtsxxqVIeDeFxVsiyrzpYAob1iOQGWA+SKy1h6KJNA1PXUOjf27icP1OO0BEdlkN0H5fbOKKxFpA/QFVhHg74sy2wKq+L7wpYCo1LAcAWSIMaYf1oi499tNDUqBNSpye6AP1rhmLzlaTS0SkSjgv8BvjTE5TtfjJDfbosrvC18KCB2Ww4UxJt3+fRj4HKsJLpAdOjUSsP37sMP1OMYYc8gYU2KMKQXeJEDeGyISivWB+KEx5jP74YB8X7jbFtV5X/hSQFRm6I6AICKRducTIhIJXAwE+gi3c4Db7du3A186WIujygyZfw0B8N6wLz72NrDdGPNvl0kB974ob1tU533hM0cxAdiHZb3C6aE7/uVsRc4QkXZYew1gDZcyM5C2hYh8BIzEGsr5EPB34AtgNtAKSAFuMMb4fedtOdtiJFYzggH2Avf6+xhnIjIUWAZsBkrth/+M1fYeUO+LCrbFGKr4vvCpgFBKKVV7fKmJSSmlVC3SgFBKKeWWBoRSSim3NCCUUkq5pQGhlFLKLQ0I5XNEpJHLiJQHXUaozBOR16u4rL/YI15uspcx0FN12+tbLCIJ5TyeJCIbRWSNiPQ5y3JiROQ+l/stRORTD5SsApjHLjmqlKcYY45gHc/9q0YuFZHBwBVAP2NMoYg0BsJqsNSqGmeMSRSRO4EXsEYuLk8McB/wOvx8Zv31Hq9QBRTdg1B+Q0RGisjX9u0nReR9EVloXwvgHjdPaQ5kGmMKAYwxmaeGMLGvt/GciKy2fzrYj8eKyH/tb/lrRGSI/XikPQDaGhFZLyKj7cfrisgsew/lY6BuJV7KCuyBKEUkSkS+F5F1Yl3/49QIxs8C7e29nhdEpI3Y14QQkXAReceef72InFfNTaoCnO5BKH/WCxgERALrReSbUwFgmw88ISI7gQXAx8aYJS7Tc4wxA0TkNqwz+K8AXgVeNsb8ICKtgHlAV+AvwEJjzHgRiQFWi8gC4F4g3xjTS0R6AesqUfcorDPDAQqAa4wxOfYezkoRmQM8DvQwxvSBn0ftPOV+AGNMTxHpgjXqbydjTEEl1q3UzzQglD/70hhzAjghIouwBif74tREY0yeiPQHhgHnAR+LyOPGmBn2LB+5/H7Zvn0h0M0a7gaAeva4WBcDV4nIY/bj4VjDOwwHXrPXt0lENlVQ74f22FrBQD/7MQGesUfrLcXas2h6ltc9FPiPvc4dIrIP6ARUtG6lfkEDQvmzsuPI/GJcGWNMCbAYWCwim7EGdJvhZv5Tt4OAwXbw/MweIO06Y0xSmcfdrrcc44CNWM1Hk4Fr7cdigf7GmCKxriQYfpbluBsaX6kq0z4I5c9G2+3xjbAGsFvjOlFEOotIR5eH+gD7XO7f5PJ7hX17PvCAyzL62DfnAQ/aQYGI9LUfX4r1IY+I9MBq9iqXMaYI+CswSES6AvWBw3Y4nAe0tmfNBaLLWYzrOjth7ckklTOvUuXSPQjlz1YD32B9QD5Vpv8BIAr4j91nUAwkA65X56sjIquwvkiNsR97CJhsNxWFYH0YTwSewuqn2GSHxF6sPos3gHfs+TfYNVXIGHNCRF4CHgP+CHwlIon283fY8xwRkR/tjulvsfY4TnkdmGLvERUDd5zqiFeqKnQ0V+WXfs3hr/bz9wIJxpjMmqxLKV+iTUxKKaXc0j0IpZRSbukehFJKKbc0IJRSSrmlAaGUUsotDQillFJuaUAopZRy6/8BFiAVt3MiYnMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Calculate non-dimensional turbine performance \n",
    "outputString = ''\n",
    "outputString += '  Wind : '\n",
    "outputString += '  Pitch :'\n",
    "outputString += '  rpm:'\n",
    "outputString += '  TSR: '\n",
    "outputString += ' Cp: '\n",
    "outputString += ' Ct '\n",
    "outputString += ' Power (MW): '\n",
    "outputString += '\\n'\n",
    "outputString += '  speed    [deg]          coeff     coeff    [N]     power \\n'\n",
    "outputString += '  [m/s]' \n",
    "outputString += '\\n'\n",
    "\n",
    "V_min = 0.0 # minimum wind or tide speed\n",
    "V_max = 25.0 # maximum wind or tide speed\n",
    "No_Int = 50 # assess performance at tshi number of values\n",
    "dV=(V_max-V_min)/(No_Int-1)\n",
    "TSRs=np.zeros(No_Int)\n",
    "Cps=np.zeros(No_Int)\n",
    "Cts=np.zeros(No_Int)\n",
    "Pws=np.zeros(No_Int)\n",
    "Thr = np.zeros(No_Int)\n",
    "Vs = np.zeros(No_Int)\n",
    "\n",
    "for ic in range(No_Int):\n",
    "    VWIND=V_min+(ic+1)*dV # find each w​for ic in range(No_Int):\n",
    "    Vs[ic] = VWIND\n",
    "    \n",
    "    TSRs[ic],Cts[ic],Cps[ic] = TurbinePerformance(dataBlade,dataAero,pitch*np.pi/180.0,ScaleFactor,Nblades,D,VWIND,RPM)\n",
    "    # dimensional power take off values\n",
    "    FCTOR1=Density*np.pi*D*D*VWIND*VWIND/8.0\n",
    "    Pws[ic] = Cps[ic]*FCTOR1*VWIND\n",
    "    Thr[ic] = Cts[ic]*FCTOR1\n",
    "    # output data\n",
    "    # outputString += '  Overall performance data\\n'\n",
    "    outputString += '  {:6.3f}'.format(VWIND)\n",
    "    outputString += '  {:6.3f}'.format(pitch)\n",
    "    outputString += '  {:6.1f}'.format(RPM)\n",
    "    outputString += '  {:6.3f}'.format(TSRs[ic])\n",
    "    outputString += '  {:.4f}'.format(Cps[ic])\n",
    "    outputString += '  {:.4f}'.format(Cts[ic])\n",
    "    outputString += '  {:.4f}'.format(Pws[ic]/1000000.0) # units are MW\n",
    "    outputString += '\\n'\n",
    " \n",
    "## now try plotting some stuff\n",
    "plt.figure()\n",
    "plt.axis([0, 25.0, 0 , 0.6])\n",
    "plt.plot(TSRs,Cps)\n",
    "plt.title('Power Coefficient vs tip speed ratio ')\n",
    "plt.xlabel('Tip Speed Ratio')\n",
    "plt.ylabel('Cp')\n",
    "##\n",
    "print(outputString) # switched off output\n",
    "with open(outputFilename,\"w\") as f:\n",
    "    f.write(outputString)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cell below was used to iteratively find the optimal RPM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8.332115542116904\n",
      "Max Cp is: 0.45364880577896627\n",
      "Cp at Umean is: 0.45359870044706124\n",
      "8.237517737100235\n",
      "3.064804473237398\n"
     ]
    }
   ],
   "source": [
    "TSRU = (np.pi*(RPM/60)*(D))/vrated\n",
    "print(TSRU)\n",
    "print(f'Max Cp is: {max(Cps)}')\n",
    "Cpsint = interp1d(TSRs, Cps)\n",
    "print(f'Cp at Umean is: {Cpsint(TSRU)}')\n",
    "TSRsint = interp1d(Cps, TSRs)\n",
    "TSRmax = TSRsint(max(Cps))\n",
    "print(TSRmax)\n",
    "Nnew = (TSRmax*vrated)/(np.pi*D)\n",
    "print(Nnew*60)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x214e1e1dfa0>]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEDCAYAAAAlRP8qAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAArIElEQVR4nO3deXxU5d338c9vJgsJ2QhZScgChDUsQgibK0plUbHVWrHaulSkVatd7tYuT73b2qf33T5dXFotVLRaFdtikSLW1pVFthB2whIIkA0SCNn3yfX8kcHGmMgEJnNm+b1fr3llZs6Vye8wyZcz17nOdYkxBqWUUr7PZnUBSiml3EMDXSml/IQGulJK+QkNdKWU8hMa6Eop5Sc00JVSyk9YGugislxEKkRkrwttfyMiO523QyJS7YESlVLKZ4iV49BF5HKgHnjBGJPdh+97ELjEGHN3vxWnlFI+xtIjdGPMOqCq63MiMlxE/iki20VkvYiM7uFbFwGveKRIpZTyEUFWF9CDpcASY8xhEZkG/B6YfW6jiKQDmcC7FtWnlFJeyasCXUQigJnAX0Xk3NOh3ZrdCvzNGOPwZG1KKeXtvCrQ6ewCqjbGTPqUNrcC93umHKWU8h1eNWzRGFMLFInI5wGk08Rz20VkFDAI2GRRiUop5bWsHrb4Cp3hPEpESkTkHuCLwD0isgvYByzs8i2LgBVGp4hUSqlPsHTYolJKKffxqi4XpZRSF86yk6JxcXEmIyPDqh+vlFI+afv27aeNMfE9bbMs0DMyMsjLy7PqxyullE8SkeO9bdMuF6WU8hMa6Eop5Sc00JVSyk9ooCullJ/QQFdKKT+hga6UUn7C5UAXEbuI7BCRNT1sExF5QkQKRWS3iEx2b5lKKaXOpy/j0B8CCoCoHrbNA7Kct2nA086vSvkkYwyHK+pZd6iSqLBgxiZHMSIhggHBdqtLU6pXLgW6iKQCC4CfAd/soclCOpeRM8BmEYkRkWRjTLn7SlWqfzk6DDtOnOVf+0/xr30nOXam8WPb7TZhePxAxiRHMTopigXjk0kbHG5RtUp9kqtH6L8FvgNE9rI9BSju8rjE+ZwGuvJ6zW0O/u/aAtbuOcnp+haC7cKM4XF85bJhXD0mgcZWBwXltRwor6OgvJZtRVW8vrOMpeuO8Mri6YxO6ulD638YY/jt24dZtbOU1++fRUx4iIf2TAWa8wa6iFwHVBhjtovIlb016+G5T0zjKCKLgcUAaWlprlepVD/62/YSXth0nPnjk5ibncyVo+KJGhD8sTbD4yO4bsJ/Hh+prOe2ZZv54rItrFg8nazEno91jDE89kYBz24oAuC9gxV89pLUftsXFdhcOSk6C7hBRI4BK4DZIvLnbm1KgKFdHqcCZd1fyBiz1BiTY4zJiY/vcW4ZpTxuze4yhscP5He3TeaGiUM+EeY9GR4fwSv3TsdmExYt20JhRf0n2nR0GH64ai/PbijizpkZxEeG8nZBRX/sglKAC4FujPmeMSbVGJNB5/Jv7xpjbu/WbDXwJedol+lAjfafK19QUdvMlqIqrpswhC7r2LpkmDPUAW5btpmjlf8J9XZHB//1t928tOUES64YzqPXj+Xq0Ql8cLCS1vYOt+6DUudc8Dh0EVkiIkucD9cCR4FCYBnwNTfUplS/W7unHGPg+onJF/T9IxIiePneaTg6DIuWbebY6QbaHB089OpOVuaX8M05I/nu3FGICNeMSaS+pZ0tRWfcvBdKderT9LnGmPeB9533n+nyvEEXblY+aM3uckYnRTIiobfz/ec3MjGSl+6dxqKlm1m0bDMjEyP54FAlP5g/hnsvH/ZRu1kj4ggNsvFOQQWXZWmXo3I/vVJUBayy6ibyjp/lugkXdnTe1eikKP78lWk0tjr44FAlP1047mNhDhAWYueyrDj+vf8UuvSj6g+WLXChlNXW7uk8zXPdhCFueb1xQ6L5+9dmcrK2mZnD43psc82YRN4uqODAyTrGJH/6cEel+kqP0FXA+sfucrJTosiIG+i21xwWH9FrmAPMHp0AwDsFp9z2M5U6RwNdBaTiqkZ2FVe77ejcVQlRA5g4NIZ/6/BF1Q800FVAWrO7s7tlwfiL7z/vqzljEthVXE1FbbPHf7bybxroKiCt2V3GpKExDI31/Fws14xNBODdA3qUrtxLA10FnKLTDewrq3XL6JYLMSoxkpSYMN7WfnTlZhroKuCs2dU5K8UCiwJdRJgzNpH1h0/T1OqwpAblnzTQVcBZs7ucqRmDSI4Os6yGa8Yk0tLewcbC05bVoPyPBroKKIdP1XHwVJ3HR7d0l5sZS2RokHa7KLfSQFcB5R+7y7EJzBufZGkdIUE2Lh8VzzsHKujo0KtGlXtooKuAYYxhze4ypmUOJiFygNXlMGdMIpV1LewurbG6FOUnNNBVwNhfXsvRygauu8CZFd3tylHx2G3C2/u120W5hwa68nvGGFbvKuPu57cRGmRj7jhru1vOiQkPISd9kPajK7fRQFd+7fCpOm5btoWvv7KD+MhQXr1vBoMjQq0u6yNzxiZy4GQdxVWN52+s1HlooCu/VN/Szs/XFjDv8fXsK6vhpzdm8/r9lzJpaIzVpX3MteOSCLHb+MGqvTj05Ki6SOcNdBEZICJbRWSXiOwTkR/30OZKEakRkZ3O24/6p1ylzu/DwtNc86sP+MO6o3xucgrvfftK7piejt3WtyXmPGFobDiP3jCWdYcqefLdw1aXo3ycK/OhtwCzjTH1IhIMbBCRN40xm7u1W2+Muc79JSrlOkeH4TsrdzMg2MbKr85kSvogq0s6r9ty09h+7CyPv3OYS9IGccVIXc1IXRhXFok2xphzq98GO2/62VB5pX/vP0XJ2Sa+O3e0T4Q5dE4F8LPPjmdkQiQPr9hBaXWT1SUpH+VSH7qI2EVkJ1AB/NsYs6WHZjOc3TJvisi4Xl5nsYjkiUheZWXlhVetVC+e21hESkwYc5wzGvqKsBA7T98+mTaH4f6X8mlt77C6JOWDXAp0Y4zDGDMJSAVyRSS7W5N8IN0YMxF4EljVy+ssNcbkGGNy4uP1Y6Vyr31lNWwpquJLM9IJsvve+f5h8RH84uYJ7Cyu5mdv7Le6HOWD+vRbb4ypBt4H5nZ7vvZct4wxZi0QLCK9r8OlVD94fuMxwoLt3Do1zepSLtj88cncc2kmf9p0nNXOWSGVcpUro1ziRSTGeT8MuAY40K1NkoiI836u83XPuL1apXpxpr6F13eV8bnJKUSHB1tdzkV5ZN5octIH8cjK3Rw6VWd1OcqHuHKEngy8JyK7gW109qGvEZElIrLE2eZmYK+I7AKeAG41xuiJU+UxL285QWt7B3fNyrC6lIsWbLfx1G2TCQ8J4ot/3EJhhYa6co1Ylbs5OTkmLy/Pkp+t/EtreweX/u+7jEqK5MV7plldjtscPlXHomVbAMPL905nZGKk1SUpLyAi240xOT1t870zR0p18+becirqWrh7VqbVpbhVVmIkKxZPxybCoqWbOXCy1uqSlJfTQFc+zRjD8g1FDIsb6JcX5IxIiGDF4ukE2YXblm1hf5mGuuqdBrryafknqtlVUsOXZ2Zg88JL+91hWHwEry6eQWiQjdv+uJm9On+66oUGuvJpz20sInJAEDdPSbW6lH6VETeQVxfPYKDzRKmGuuqJBrryWeU1Tby59yRfyBnKwFBXpiXybWmDw1mxeDoDgm08phceqR5ooCuf9eKm4xhj+PLMDKtL8ZihseHMHZfEnpIaXYtUfYIGuvJJdc1tvLz1BNeMSWRobLjV5XhUdko0Da0Ojp5usLoU5WU00JVPWr7hGNWNbTwwe4TVpXhcdko00Dl3jVJdaaArn1PT2MYfNxxlzthEJqTGWF2Ox2UlRBAaZGNPiQa6+jgNdOVzlq0/Sl1zO9+cM9LqUiwRZLcxJjmKPTrSRXWjga58ypn6Fp7bWMSCCcmMSY6yuhzLjE+JZl9ZrZ4YVR+jga58yh/WHaWpzcE3rsmyuhRLZadEUd/SzvGqRqtLUV5EA135jIq6Zl7YdIyFk1IYkRDYE1WdOzGq3S6qKw105TN+/94R2hyGh64O7KNzgJGJkYQE2fSKUfUxGujKJ5TXNPHylhPcNDmFjLiBVpdjuWC7jTFJkTrSRX2MKysWDRCRrc4FoPeJyI97aCMi8oSIFIrIbhGZ3D/lqkD11LuFGAwPztaj83PGpUSzt6wGXUtGnePKEXoLMNu5APQkYK6ITO/WZh6Q5bwtBp52Z5EqsBVXNfKXvGK+MHVowF0V+mnGp0RT19zOCT0xqpzOG+imU73zYbDz1v2QYCHwgrPtZiBGRJLdW6oKVE++exgR4YGr9Oi8q/F6YlR141IfuojYRWQnUEHnmqJbujVJAYq7PC5xPqfURSmuamRlfilfnJZGUvQAq8vxKiMTIwm2iwa6+ohLgW6McRhjJgGpQK6IZHdr0tPKAp/o2BORxSKSJyJ5lZWVfS5WBZ5nNxRhE7jv8uFWl+J1QoJsjEqK1JEu6iN9GuVijKkG3gfmdttUAgzt8jgVKOvh+5caY3KMMTnx8f63XJhyr7MNrby6rZgbJqbo0XkvxqdEs7e0Vk+MKsC1US7xIhLjvB8GXAMc6NZsNfAl52iX6UCNMabc3cWqwPLSluM0tTlYfPkwq0vxWtkp0dQ0tVFytsnqUpQXcGWZl2TgTyJip/M/gL8YY9aIyBIAY8wzwFpgPlAINAJ39VO9KkA0tzl4/sPjXDEynlFJgX1V6KfpemJURwCp8wa6MWY3cEkPzz/T5b4B7ndvaSqQrdpRyun6Fu7To/NPNTIxkiBb54nR+eN1YFmg0ytFldfp6DAsXX+U7JQoZgwfbHU5Xm1AsJ2RiXpiVHXSQFde550DFRytbODey4Yh0tMAKtVV54lRvWJUaaArL7Rs3VFSYsJYoF0ILslOjeZsYxul1XpiNNBpoCuvsuPEWbYeq+LuSzMJsuuvpyuyh3Qu9KHdLkr/YpRXWbb+KFEDgrh16tDzN1YAjEmOwm7TK0aVBrryIsfPNPDPvSe5fXo6A0NdGVGroPPEaFZCBHtKa60uRVlMA115jWc3FBFks3HnzAyrS/E541Oi2acnRgOeBrryClUNrfwlr5gbLxlCQpRe5t9X2SnRnGlopbym2epSlIU00JVX+PuOUprbOrjnUr2Q6ELoGqMKNNCVl3h9ZynZKVF6mf8FGpschU10pEug00BXljtaWc/ukhpunKRT6F+osBA7WQl6xWig00BXlnt9ZxkicN2EIVaX4tOyU6LZUVxNY2u71aUoi2igK0sZY1i9q4wZwwbrnOcX6dbcoVQ3tvGHD45aXYqyiAa6stTukhqKTjewcJIenV+sqRmxXDchmWc+OKLTAAQoDXRlqdd3lhFitzE3W+dtcYfvzR8DwM/XFlhcibKCBrqyjKPD8I/dZVw1Op7osGCry/ELKTFhLLliOGt2l7O1qMrqcpSHubIE3VAReU9ECkRkn4g81EObK0WkRkR2Om8/6p9ylT/58MhpKutadHSLmy25YjjJ0QP48T/24ejQK0cDiStH6O3At4wxY4DpwP0iMraHduuNMZOct5+4tUrll17fWUZkaBBXjU6wuhS/EhZi55F5o9lXVsvfthdbXY7yoPMGujGm3BiT77xfBxQAekilLkpzm4N/7j3J3OwkBgTbrS7H79wwcQg56YP45VsHqW1us7oc5SF96kMXkQw61xfd0sPmGSKyS0TeFJFxvXz/YhHJE5G8ysrKvler/Ma7Byqob2lnoXa39AsR4dHrx3GmoZWn3i20uhzlIS4HuohEACuBh40x3efpzAfSjTETgSeBVT29hjFmqTEmxxiTEx8ff4ElK3+wakcp8ZGhumZoPxqfGs3np6Ty3MYiik43WF2O8gCXAl1EgukM85eMMa91326MqTXG1DvvrwWCRSTOrZUqv1HT2Mb7Byu5fsIQ7DZdM7Q/ffvaUYQG2fnZG/utLkV5gCujXAR4Figwxvy6lzZJznaISK7zdc+4s1DlP97cW06ro4MbL9GLifpbQuQAHpw9grcLKvjt24d0vnQ/58qyMLOAO4A9IrLT+dz3gTQAY8wzwM3AV0WkHWgCbjX6m6N68frOMjLjBjLeOeWr6l9fuWwYh07V89u3D9PuMHzrMyNxHn8pP3PeQDfGbAA+9d03xjwFPOWuopT/OlnTzOaiMzx0dZaGiofYbcIvb55ASJDw1HuFtDo6+N680frv74d04UblUSvzSzAGHd3iYTab8LMbxxNst7F03VFa2zt49PqxGup+RgNdeUx1YytL1x3lipHxZMYNtLqcgGOzCT++YRzBdhvPbiiizdHBTxdmY9MT035DA115zJPvFlLX3Mb3nRNIKc8TEX64YAzBdhvPfHCENkcHP//cBB1t5Cc00JVHHD/TwAubjnFLzlBdZs5iIsJ3544ixC488W4hs0bEaReYn9DZFpVH/OKfBwm22/jmnJFWl6LoDPWHrxlJdFgwGwtPW12OchMNdNXvth+v4o095dx3+XASonRVIm9hswlTMwax7dhZq0tRbqKBrvqVMYbH3iggITKUey/PtLoc1c3UjFiKTjdQUddsdSnKDTTQVb9au+ckO05U8+3PjCI8RE/ZeJupmbEA5OlRul/QQFf9pqXdwf/8s4DRSZHcNCXV6nJUD7KHRBMWbNfVjfyEBrrqNy9uOk5xVRM/WDBGh8V5qZAgG5ekxWig+wkNdNUvqhtbeeKdw1wxMp7LsnSqZG82NSOWgpO1uhCGH9BAV/3iqXcLqW9p14uIfMC0zFiMge3HtR/d12mgK7drbG1nxbZibpg4RC8i8gGXpA0iyCba7eIHNNCV2/1jVxn1Le3cPj3d6lKUC8JC7GSnRLNNA93naaArt3t5azFZCRFMSR9kdSnKRdMyY9ldUkNzm8PqUtRFcGXFoqEi8p6IFIjIPhF5qIc2IiJPiEihiOwWkcn9U67ydvvKathVXM2i3DSdmtWHTM2IpdXRwc7iaqtLURfBlSP0duBbxpgxwHTgfhEZ263NPCDLeVsMPO3WKpXPWLG1mJAgG5+brJM9+ZKcjM5PU9rt4tvOG+jGmHJjTL7zfh1QAHT/a10IvGA6bQZiRCTZ7dUqr9bY2s6qHaUsGJ9MTHiI1eWoPogJD2FUYiRbj2mg+7I+9aGLSAZwCbCl26YUoLjL4xI+GfrKz63ZXU5dSzuLctOsLkVdgNzMWPKPn6Xd0WF1KeoCuRzoIhIBrAQeNsbUdt/cw7d8YpFoEVksInkikldZWdm3SpXXe2XrCYbHD2Rqhp4M9UVTM2NpaHWwv7z7n7fyFS4FuogE0xnmLxljXuuhSQkwtMvjVKCseyNjzFJjTI4xJic+Xq8e9CcF5bXsOKEnQ31ZbkbnRF06Ht13uTLKRYBngQJjzK97abYa+JJztMt0oMYYU+7GOpWXW7H1BCF2GzdN1km4fFVS9ADSYsPZpv3oPsuV+UxnAXcAe0Rkp/O57wNpAMaYZ4C1wHygEGgE7nJ7pcprNbU6eG1HKfPGJzFooJ4M9WVTM2J572AFxhj9pOWDzhvoxpgN9NxH3rWNAe53V1HKt7yxp5y6Zj0Z6g9yMwexMr+EI5X1jEjQaRt8jV4pqi7aK1tPMCxuINOciyUo35WbORiArUU6UZcv0kBXF+XgyTq2Hz+rJ0P9RMbgcOIiQrUf3UdpoKuL8sq5k6G6IpFfEBFyMwfpSBcfpYGuLlhjazuv5ZfwmXGJxOrJUL8xNSOW0uomSqubrC5F9ZEGurpgK/NLqW1u586ZGVaXotwo13kuROd18T0a6OqCdHQYnt9YxITUaJ0m18+MTooiOiyY9w5WWF2K6iMNdHVB1h2u5EhlA3fNytCToX7GbhMWThrCm3tPcrah1epyVB9ooKsLsnzjMeIjQ1kwfojVpah+cNu0NFrbO1iZX2J1KaoPNNBVnxVW1LHuUCV3TE8nJEh/hfzR6KQopqQP4uUtJ+i8blD5Av1rVH323MZjhATZuG2aXhnqz744LY2jpxvYdPSM1aUoF2mgqz6pbmzltfxSbpw0hLiIUKvLUf1o/vhkosOCeWnLCatLUS7SQFd9smJbMU1tDu6alWl1KaqfDQi2c/OUVP617ySVdS1Wl6NcoIGuXNbu6OCFD48xY9hgxiRHWV2O8oBFuWm0OQx/3V58/sbKchroymVv7TtFWU0zd1+qR+eBYkRCBNOHxfLK1hN0dOjJUW+nga5ctnxjEWmx4cwenWB1KcqDbpuWTnFVE+sLT1tdijoPDXTlkl3F1Ww/fpY7Z2Zgt+mFRIHk2nGJDB4Ywkubj1tdijoPV5agWy4iFSKyt5ftV4pIjYjsdN5+5P4yldWe21hERGgQn8/RWRUDTWiQnZtzUnnnQAUna5qtLkd9CleO0J8H5p6nzXpjzCTn7ScXX5byJidrmlmzu5zP56QSOSDY6nKUBW7LTcPRYXh1m54c9WbnDXRjzDpAp10LYEvXHcUAd+tQxYCVPnggl2XFsWLbCdodHVaXo3rhrj70GSKyS0TeFJFxvTUSkcUikicieZWVlW760ao/na5v4eWtx7lxUgpDY8OtLkdZ6IvT0iivaeb9g/q3663cEej5QLoxZiLwJLCqt4bGmKXGmBxjTE58fLwbfrTqb39cX0RLewdfu2q41aUoi109JpH4yFB+9e9DnDjTaHU5qgcXHejGmFpjTL3z/logWETiLroyZbnqxlZe3HSMBeOTGR4fYXU5ymLBdhuP3ZhNSVUj1/52HX/68JiOTfcyFx3oIpIkzgmxRSTX+Zo6m48feG7jMRpaHTwwe4TVpSgvce24JN76xuXkZsby6Op9LFq2WY/WvYgrwxZfATYBo0SkRETuEZElIrLE2eRmYK+I7AKeAG41Ot+mz6trbuO5jUV8Zmwio5P0Mn/1H0Niwnj+rqn84qYJ7C+r1aN1LxJ0vgbGmEXn2f4U8JTbKlJe4cXNx6ltbtejc9UjEeGWqUO5bGQcj6zcw6Or97FmdxmPzBvNlPRYq8sLWHqlqPqExtZ2/ri+iCtGxjMhNcbqcpQXS452Hq3fPIGi0w3c9PQm7npuK3tKaqwuLSBpoKtPeGVrMVUNrTyoR+fKBSLCLTlDWfedq3hk3mh2FFdz/VMbuO/FPA6crLW6vICiga4+prnNwdJ1R5g+LJacDP3orFwXHhLEkiuGs/47V/GNa0byYeEZ5j2+nodX7KCp1WF1eQFBA119zF+3l3CqtoUHZ2dZXYryUZEDgnnomizWf/cq7rt8OKt2lvH8h8esLisgaKCrj7Q5Onjm/SNckhbDzOGDrS5H+biY8BAemTeaq0bF88wHR6htbrO6JL+nga4+8pe8Ykqrm3hw9giclxYoddG+OWcUNU1tLN9QZHUpfk8DXQFQ09TGr/91iKkZg7hqlC5godxnfGo0c8cl8ez6Is42tFpdjl/TQFcAPPnOYaoaW3n0+nF6dK7c7htzRlLf2s7S9UetLsWvaaArjlTW8/yHx/hCzlCyU6KtLkf5oVFJkVw/YQjPbzzG6foWq8uxVJujg/qW9n55bQ10xWNr9hMWbOdbnxlldSnKjz18TRYt7Q6efv+I1aV4jDGGotMNrNpRyo//sY/P/X4j2Y++xbJ1/fNJ5byX/iv/9t6BCt47WMkP5o8hPjLU6nKUHxsWH8FNk1N5cfNx7r1sGEnRA6wuya2a2xwUVtRz4GQdB0/WcuBkHbtLaqhp6hzdExZsJzslijump/fbKDIN9ADW2t7BT9/YT2bcQL48M8PqclQA+PrVWazaWcpT7x3msRvHW13OBTHGUFbTTEFZLQXlncF94GQtRacbODc/WUiQjayECOaPT2ZiajQTh8aQlRBBkL1/O0U00APYC5uOcbSygeV35hASpL1vqv8NjQ3nlpyhvLqtmPsuH+71q2A1tzk4fKqegvJa9pf/J8DPHXUDpMWGMzopkgXjkxmVFMWopEgyBof3e3j3RAM9QJ2ub+Hxdw5z5ah4Zo9OtLocFUAemD2Cv24v4Yl3DvPLz0+0upyPnKlvYU9pjTO46ygo7zzqdjgPu8ND7IxKimTBhGTGJEcxNjmSUUlRRIR6T4x6TyXKo371r0M0tTr44YKxVpeiAkxydBi3T0vnT5uOsWhaGpPTBnn05xtjONPQSkF5LbtLathTUsOe0hpKq5s+apMSE8aY5CjmZycxJjmK0clRpMeGY7N595De8wa6iCwHrgMqjDHZPWwX4HFgPtAI3GmMyXd3ocp99pXVsGLbCe6elcmIBF1aTnneV68czupdpXzu9x9yw8QhfGPOSDLjBrrltY0xnKptoaC8lkOn6iivaaairplTtS2crGmmsq6FVkfHR+0zBodzSVoMX56ZTnZKNOOGRBMdFuyWWjzNlSP05+lcwOKFXrbPA7Kct2nA086vygt1dBgefX0fg8JD+PrVOgGXskZ8ZCjvfPNKlq4/wvINx3hjTzm35KTy4OwshsSE9em1ztS38O6Bio+6SQ6crOVs43/6uCNCg0iICiUpagC5mbEkRg0gMSqUUYmRjEvx3fDuiSsrFq0TkYxPabIQeMG57NxmEYkRkWRjTLm7ilTuszK/hLzjZ/nFzRP86hdZ+Z7o8GD+69rR3Dkzk9+9V8jLW06wMr+UO6anc+fMjPOeMD10qo7lG4r4+45SWto7GBBsY1RiJNeOc3aTJEUyOimK6PDA+T13Rx96ClDc5XGJ8zkNdC9T3djKz988wJT0Qdw8OdXqcpQCOo/W//uGcXzlskwef/swz20s4tkNRWTGDeTSEXHMGhHHjOGDiQ4LpqPD8MHhSpZvKGL94dOEBtn43ORUbp+exuikKOxe3sfd39wR6D39C/a4WqyILAYWA6SlpbnhR6u++MVbB6lpauOxG7O9/uSOCjypg8L55ecn8uDsLN4uOMWGwtOszC/hxc3HsQlMHBpDbVMbRyobSIgM5b+uHcWi3DRiB4ZYXbrXcEeglwBDuzxOBcp6amiMWQosBcjJydElwj1oZ3E1r2w9wV0zMxmTHGV1OUr1Km1wOHdfmsndl2bS2t7BzuJqNhyuZH3haaLDgvnNFyayYPwQvXaiB+4I9NXAAyKygs6ToTXaf+5dHB2GH67aQ3xEKN+YoydCle8ICbKRmxlLbmYs39S5hs7LlWGLrwBXAnEiUgI8CgQDGGOeAdbSOWSxkM5hi3f1V7Hqwry85Th7S2t5YtElRA4InBNESgUaV0a5LDrPdgPc77aKlFtV1rXwi7cOMmvEYK6fkGx1OUqpfqSdUH7u52sLaG5z8JOF2bpwhVJ+TgPdj20+eobXdpSy+PJhDI/XK0KV8nca6H6quc3BD/6+h5SYMB64Sk+EKhUIdHIuP/XLtw5ypLKBF+7OJSzEbnU5SikP0CN0P7TpyBmWbyzijunpXD4y3upylFIeooHuZ+qa2/j2X3eRHhvO9+aPtrocpZQHaZeLn3lsTQHlNU38dckMwkP07VUqkOgRuh95e/8pXs0rZskVw5mSHmt1OUopD9NA9xNVDa088toeRidF8tA1OqpFqUCkn8n9gDGdc7XUNLXy4j25hAbpqBalApEeofuB1bvKWLvnJN+YM1JnUlQqgGmg+7gTZxr54aq9TEkfxH2XD7e6HKWUhTTQfVhTq4P7/rwdmwi/uWVSwK/WolSg0z50H2WM4Qd/38OBk7Usv3MqaYM/ff1FpZT/0yN0H/Xnzcd5bUcpD12dxVWjEqwuRynlBTTQfdD242f5yZr9XDUqnq/P1iGKSqlOLgW6iMwVkYMiUigij/Sw/UoRqRGRnc7bj9xfqoLOBSu+9tJ2kqPD+O0XLtHFnpVSH3FlCTo78DtgDp0LQm8TkdXGmP3dmq43xlzXDzUqp3ZHBw+8nE91YxuvfW0q0eG6nJxS6j9cOULPBQqNMUeNMa3ACmBh/5alevK//zzAlqIq/u9nxzNuSLTV5SilvIwrgZ4CFHd5XOJ8rrsZIrJLRN4UkXE9vZCILBaRPBHJq6ysvIByA9df8opZtr5zStybpqRaXY5Sygu5Eug9ddKabo/zgXRjzETgSWBVTy9kjFlqjMkxxuTEx+s83a5698ApvvfaHi4dEcf/uW6s1eUopbyUK4FeAgzt8jgVKOvawBhTa4ypd95fCwSLSJzbqgxg+SfO8rWX8hmTHMkzd0whJEgHJimleuZKOmwDskQkU0RCgFuB1V0biEiSOJeUF5Fc5+uecXexgaawop67n99GYtQAnrszl4hQvQ5MKdW78yaEMaZdRB4A3gLswHJjzD4RWeLc/gxwM/BVEWkHmoBbjTHdu2VUH5ysaebLy7cSZBNeuDuX+MhQq0tSSnk5sSp3c3JyTF5eniU/29vVNLXxhT9soriqkVfvm0F2io5oUUp1EpHtxpicnrbpZ3gv09zm4N4X8jhSWc/yO6dqmCulXKaB7kWa2xwsfnE7W4uqePzWSVyWpSOBlFKu00D3Eo2t7dzzfB6bi87wvzeNZ+Gknob6K6VU7zTQvUBdcxt3PbeN/BNn+fUtE/nsJXrhkFKq7zTQLVbT2MaXntvKvtIanlw0mQUTkq0uSSnlozTQLVTV0Modz27h8Kl6nr59CnPGJlpdklLKh2mgW6SyroXb/7iFY2caWPqlKVypi1QopS6SBroF9pbWcN+L26lqaOW5O6cyc4TOkqCUung6MYiHvZZfwk1Pf0iHMaxYPF3DXCnlNnqE7iFtjg5+9kYBz394jGmZsfzui5OJi9DL+ZVS7qOB7gGn61v42kv5bC2q4u5ZmXxv/miC7frhSCnlXhro/Sz/xFnufymfqoZWfvMFHWOulOo/Guj9ZFdxNU+9V8i/958iJSaMlV+dqfOyKKX6lQa6m207VsWT7xay7lAl0WHBPHxNFnfNyiQ6TBd0Vkr1Lw10N2huc7Dp6Bmeef8IW4qqGDwwhO/OHc3t09OIHKBBrpTyDJcCXUTmAo/TucDFH40x/9Ntuzi3zwcagTuNMflurtVrtLQ72FVcw6YjZ9h89Az5J87S0t5BYlQoP7puLIty0wgLsVtdplIqwJw30EXEDvwOmEPn+qLbRGS1MWZ/l2bzgCznbRrwtPOrz2tsbefwqXoOnarj0Kk69pXVkn/iLM1tHYjA2OQobp+ezoxhg7k0K44BwRrkSilruHKEngsUGmOOAojICmAh0DXQFwIvOJed2ywiMSKSbIwpd3fBZxtaOVHVSHiInfDQIAaG2AkLsRNit+Fc1rRPjDGcaWilrLqJsuomSs42UVbdzImqBg6dqqf4bCPnFnUKCbIxMjGCRblpTB82mGmZscSEh7h5D5VS6sK4EugpQHGXxyV88ui7pzYpgNsD/cMjZ7j/5U/25gTZ5KNgt9vko1uQ86sBOjoMDmPo6ACH835tUxst7R0fe62wYDtDY8OYkBrNzVNSGZkYycjECNIHD8Ru6/t/Gkop5QmuBHpPCdZ9IVJX2iAii4HFAGlpaS786E+amjGI5Xfm0NDioKnVQUNrO42tDhpaOr+2d3Tg6DC0O8xHod3uMCBgl85wt4lgt4HdJkSEBpESE8YQ5y11UBjRYcEXdLSvlFJWciXQS4ChXR6nAmUX0AZjzFJgKXQuEt2nSp0SogYwO2rAhXyrUkr5NVeuP98GZIlIpoiEALcCq7u1WQ18STpNB2r6o/9cKaVU7857hG6MaReRB4C36By2uNwYs09Elji3PwOspXPIYiGdwxbv6r+SlVJK9cSlcejGmLV0hnbX557pct8A97u3NKWUUn2hU/4ppZSf0EBXSik/oYGulFJ+QgNdKaX8hAa6Ukr5CTHmgq7vufgfLFIJHD9PszjgtAfK8Sa6z4FB9zkw9Mc+pxtj4nvaYFmgu0JE8owxOVbX4Um6z4FB9zkweHqftctFKaX8hAa6Ukr5CW8P9KVWF2AB3efAoPscGDy6z17dh66UUsp13n6ErpRSykUa6Eop5Se8NtBFZK6IHBSRQhF5xOp6PEFEjonIHhHZKSJ5VtfTH0RkuYhUiMjeLs/Fisi/ReSw8+sgK2t0t172+b9FpNT5Xu8UkflW1uhOIjJURN4TkQIR2SciDzmf99v3+VP22aPvs1f2oYuIHTgEzKFzNaRtwCJjzP5P/UYfJyLHgBxjjN9efCEilwP1dC4qnu187hdAlTHmf5z/eQ8yxnzXyjrdqZd9/m+g3hjz/6ysrT+ISDKQbIzJF5FIYDtwI3Anfvo+f8o+34IH32dvPULPBQqNMUeNMa3ACmChxTUpNzDGrAOquj29EPiT8/6f6PxD8Bu97LPfMsaUG2PynffrgAI6F4332/f5U/bZo7w10FOA4i6PS7DgH8cCBviXiGx3LqgdKBLPLVno/JpgcT2e8oCI7HZ2yfhN90NXIpIBXAJsIUDe5277DB58n7010KWH57yvb8j9ZhljJgPzgPudH9WVf3oaGA5MAsqBX1laTT8QkQhgJfCwMabW6no8oYd99uj77K2BXgIM7fI4FSizqBaPMcaUOb9WAH+ns+spEJxy9kGe64ussLiefmeMOWWMcRhjOoBl+Nl7LSLBdAbbS8aY15xP+/X73NM+e/p99tZA3wZkiUimiIQAtwKrLa6pX4nIQOfJFERkIPAZYO+nf5ffWA182Xn/y8DrFtbiEeeCzemz+NF7LSICPAsUGGN+3WWT377Pve2zp99nrxzlAuAc3vNbwA4sN8b8zNqK+peIDKPzqBw6F+9+2R/3WUReAa6kc1rRU8CjwCrgL0AacAL4vDHGb04i9rLPV9L5MdwAx4D7zvUv+zoRuRRYD+wBOpxPf5/OPmW/fJ8/ZZ8X4cH32WsDXSmlVN94a5eLUkqpPtJAV0opP6GBrpRSfkIDXSml/IQGulJK+QkNdKWU8hMa6Eop5Sf+P5GXf1yhTv91AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(Vs, Pws)\n",
    "#plt.axis([0, 5, -2 , 10])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pitch Angle for Rated Speed\n",
    "The following is an ammendment to the initial code which iteratively searches for the required pitch angle to achieve rated power above rated speed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate non-dimensional turbine performance \n",
    "outputString = ''\n",
    "outputString += '  Wind : '\n",
    "outputString += '  Pitch :'\n",
    "outputString += '  rpm:'\n",
    "outputString += '  TSR: '\n",
    "outputString += ' Cp: '\n",
    "outputString += ' Ct '\n",
    "outputString += ' Power (MW): '\n",
    "outputString += '\\n'\n",
    "outputString += '  speed    [deg]          coeff     coeff    [N]     power \\n'\n",
    "outputString += '  [m/s]' \n",
    "outputString += '\\n'\n",
    "\n",
    "V_min = 0.0 # minimum wind or tide speed\n",
    "V_max = 20.0 # maximum wind or tide speed\n",
    "No_Int = 100 # assess performance at tshi number of values\n",
    "dV=(V_max-V_min)/(No_Int-1)\n",
    "TSRs=np.zeros(No_Int)\n",
    "Cps=np.zeros(No_Int)\n",
    "Cts=np.zeros(No_Int)\n",
    "Pws=np.zeros(No_Int)\n",
    "pitchNEW = 0\n",
    "TSRnew = 0\n",
    "Ctnew = 0\n",
    "Cpnew = 0\n",
    "gen_eff = 0.965\n",
    "gbox_eff = 0.985\n",
    "\n",
    "Prated = 10000000 # W\n",
    "Pturb = Prated/(gen_eff*gbox_eff)\n",
    "pitches = []\n",
    "pitches2 = []\n",
    "Vs = []\n",
    "\n",
    "for ic in range(No_Int):\n",
    "    VWIND=V_min+(ic+1)*dV # find each w​for ic in range(No_Int):\n",
    "   \n",
    "    TSRs[ic],Cts[ic],Cps[ic] = TurbinePerformance(dataBlade,dataAero,pitch*np.pi/180.0,ScaleFactor,Nblades,D,VWIND,RPM)\n",
    "    # dimensional power take off values\n",
    "    FCTOR1=Density*np.pi*D*D*VWIND*VWIND/8.0\n",
    "    Pws[ic] = Cps[ic]*FCTOR1*VWIND\n",
    "    # output data\n",
    "    # outputString += '  Overall performance data\\n'\n",
    "    outputString += '  {:6.3f}'.format(VWIND)\n",
    "    outputString += '  {:6.3f}'.format(pitch)\n",
    "    outputString += '  {:6.1f}'.format(RPM)\n",
    "    outputString += '  {:6.3f}'.format(TSRs[ic])\n",
    "    outputString += '  {:.4f}'.format(Cps[ic])\n",
    "    outputString += '  {:.4f}'.format(Cts[ic])\n",
    "    outputString += '  {:.4f}'.format(Pws[ic]/1000000.0) # units are MW\n",
    "    outputString += '\\n'\n",
    "    \n",
    "    pitchNEW = pitch\n",
    "    \n",
    "    \n",
    "    while Pws[ic] > Pturb:\n",
    "        pitchNEW = pitchNEW + 0.01\n",
    "        TSRnew,Cts[ic],Cpnew = TurbinePerformance(dataBlade,dataAero,pitchNEW*np.pi/180.0,ScaleFactor,Nblades,D,VWIND,RPM)\n",
    "        Pws[ic] = Cpnew*FCTOR1*VWIND\n",
    "        \n",
    "    if Pws[ic]<10400000 and VWIND>7.5:\n",
    "        pitchint = interp1d(Vs,pitches2,fill_value='extrapolate')\n",
    "        pitchNEW = pitchint(VWIND)\n",
    "        TSRnew,Cts[ic],Cpnew = TurbinePerformance(dataBlade,dataAero,pitchNEW*np.pi/180.0,ScaleFactor,Nblades,D,VWIND,RPM)\n",
    "        Pws[ic] = Cpnew*FCTOR1*VWIND\n",
    "        while Pws[ic] < 10400000 or Pws[ic] > 10600000:\n",
    "            pslim = pitchNEW\n",
    "            ps = np.arange(0,90,0.001)\n",
    "            for i in ps:\n",
    "                pitchNEW = i\n",
    "                TSRnew,Cts[ic],Cpnew = TurbinePerformance(dataBlade,dataAero,pitchNEW*np.pi/180.0,ScaleFactor,Nblades,D,VWIND,RPM)\n",
    "                Pws[ic] = Cpnew*FCTOR1*VWIND\n",
    "                if Pws[ic] > 10400000 and Pws[ic] < 10600000:\n",
    "                    break\n",
    "            '''\n",
    "            if Pws[ic] < 10400000:\n",
    "                pitchNEW = pitchNEW - 0.0001\n",
    "                TSRnew,Ctnew,Cpnew = TurbinePerformance(dataBlade,dataAero,pitchNEW*np.pi/180.0,ScaleFactor,Nblades,D,VWIND,RPM)\n",
    "                Pws[ic] = Cpnew*FCTOR1*VWIND\n",
    "            if Pws[ic] > 10600000:\n",
    "                pitchNEW = pitchNEW + 0.00011\n",
    "                TSRnew,Ctnew,Cpnew = TurbinePerformance(dataBlade,dataAero,pitchNEW*np.pi/180.0,ScaleFactor,Nblades,D,VWIND,RPM)\n",
    "                Pws[ic] = Cpnew*FCTOR1*VWIND\n",
    "                print(pitchNEW)\n",
    "                '''\n",
    "    print(f'Pitch at V = {VWIND} is: {pitchNEW}')\n",
    "    pitches.append(pitchNEW)\n",
    "    if VWIND>7.3:\n",
    "        Vs.append(VWIND)\n",
    "        pitches2.append(pitchNEW)\n",
    "        \n",
    "## now try plotting some stuff\n",
    "plt.figure()\n",
    "plt.axis([0, 25.0, 0 , 0.6])\n",
    "plt.plot(TSRs,Cps)\n",
    "plt.title('Power Coefficient vs tip speed ratio ')\n",
    "plt.xlabel('Tip Speed Ratio')\n",
    "plt.ylabel('Cp')\n",
    "##\n",
    "print(outputString) # switched off output\n",
    "with open(outputFilename,\"w\") as f:\n",
    "    f.write(outputString)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Pshift2(V_min, V_max, No_Int, cut_in=2.6, P=Pws, Pturb=Pturb, vrated=vrated):\n",
    "    dV=(V_max-V_min)/(No_Int-1)\n",
    "    vc=np.zeros(No_Int)\n",
    "    Pws_new=np.zeros(No_Int)\n",
    "    for ic in range(No_Int):\n",
    "        vc[ic] = V_min+(ic+1)*dV\n",
    "        if vc[ic]<=cut_in:\n",
    "            Pws_new[ic]=0\n",
    "        else:\n",
    "            Pws_new[ic] = P[ic-1]\n",
    "        \n",
    "    return(Pws_new)\n",
    "\n",
    "powers = Pshift(0,20,100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell outputs the power profile and pitch profile graph, indicating cut-in, cut-out and rated speed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax1 = plt.subplots()\n",
    "\n",
    "ax2 = ax1.twinx()\n",
    "ax1.axvline(x=vrated, color='coral')\n",
    "ax1.axvline(x=2.6, color='coral')\n",
    "ax1.axvline(x=20, color='coral')\n",
    "#ax1.axhline(y=Pturb/1e6)\n",
    "ax1.plot(V+dV, Pws/1e6, label='Power', color='blue')\n",
    "ax2.plot(V+dV, pitches, 'g-', label='Pitch angle')\n",
    "\n",
    "ax1.set_xlabel('Wind Speed (m/s)')\n",
    "ax1.set_ylabel('Power (MW)', color='blue')\n",
    "ax2.set_ylabel('Pitch angle (degrees)', color='g')\n",
    "plt.text(1.5, 28.5, 'Cut-in \\nspeed', fontsize=10, bbox=dict(facecolor='red', alpha=0.3))\n",
    "plt.text(5.9, 28.5, 'Rated \\nspeed', fontsize=10, bbox=dict(facecolor='red', alpha=0.3))\n",
    "plt.text(18.5, 28.5, 'Cut-out \\nspeed', fontsize=10, bbox=dict(facecolor='red', alpha=0.3))\n",
    "\n",
    "plt.title('')\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "V = np.linspace(0,20,100)\n",
    "plt.axis([0, 25.0, 0 , 1.1e7])\n",
    "plt.plot(V,Pws)\n",
    "plt.axvline(x=7.061761459072345  )\n",
    "#print(Pws)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(V,pitches)\n",
    "plt.axis([2.6, 20.0, -2 , 30])\n",
    "plt.xlabel('Wind Speed (m/s)')\n",
    "plt.ylabel('Pitch Angle (degrees)')\n",
    "plt.title('Pitch Angle vs Wind Speed between cut-in and cut-out speeds')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PwrR = []\n",
    "count = 0\n",
    "for i in V:\n",
    "    if i >= vrated and i<=19.7:\n",
    "        PwrR.append(Pws[count])\n",
    "    count += 1\n",
    "#print(PwrR)\n",
    "std = statistics.stdev(PwrR, Pturb)\n",
    "print(std)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Examine power capture behaviour\n",
    "# Define Cut in speed\n",
    "# Define cut out speed\n",
    "# operational conditions, change these value for tidal turbine and to see how they \n",
    "\n",
    "V_cutin = 2.6 # m/s    \n",
    "V_cutout = 20.0 #m/s\n",
    "\n",
    "#  generator + gearbox _efficiency\n",
    "gen_eff = 0.965\n",
    "gbox_eff = 0.985\n",
    "\n",
    "# rated power\n",
    "Prated = 10000000 # W\n",
    "Pturb = Prated/(gen_eff*gbox_eff)\n",
    "\n",
    "vc=np.zeros(No_Int)\n",
    "rc=np.zeros(No_Int)\n",
    "for ic in range(No_Int):\n",
    "    vc[ic]=V_min+(ic+1)*dV # find each w for ic in range(No_Int):\n",
    "    rc[ic]=Rayleigh(vc[ic],vrated)\n",
    "    if vc[ic] < V_cutin:\n",
    "        Pws[ic] = 0.0\n",
    "    elif vc[ic] >= V_cutout:\n",
    "        Pws[ic] = 0.0\n",
    "    elif Pws[ic] > Pturb:\n",
    "        Pws[ic] = Pturb\n",
    "\n",
    "plt.figure()\n",
    "plt.axis([0, V_max, 0 , 1.2*Pturb])\n",
    "plt.plot(vc,Pws)\n",
    "plt.title('Wind speed vs Regulated Power output ')\n",
    "plt.xlabel('Wind Speed (m/s)')\n",
    "plt.ylabel('MW')   \n",
    "\n",
    "plt.figure()\n",
    "plt.axis([0, V_max, 0 , 0.1])\n",
    "plt.plot(vc,rc)\n",
    "plt.title('Wind speed vs Rayleight probability ')\n",
    "plt.xlabel('Wind Speed (m/s)')\n",
    "plt.ylabel('f(V)')   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function alters the power values slightly to ensure below cut-in the power is zero and above rated speed the power is rated power"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Pshift(V_min, V_max, No_Int, cut_in=2.6, P=Pws, Pturb=Pturb, vrated=vrated):\n",
    "    dV=(V_max-V_min)/(No_Int-1)\n",
    "    vc=np.zeros(No_Int)\n",
    "    Pws_new=np.zeros(No_Int)\n",
    "    for ic in range(No_Int):\n",
    "        vc[ic] = V_min+(ic+1)*dV\n",
    "        if vc[ic]<cut_in:\n",
    "            Pws_new[ic]=0\n",
    "        elif vc[ic]>=cut_in and vc[ic]<vrated:\n",
    "            Pws_new[ic]=P[ic]\n",
    "        elif vc[ic]>=vrated:\n",
    "            Pws_new[ic] = Pturb\n",
    "    return(Pws_new)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell defines the Weibull distribution used for energy capture calculations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Weibull(U, Uhat=vrated):\n",
    "    exp = np.exp(-np.pi*(U/(2*Uhat))**2)\n",
    "    frac = (np.pi*U)/(2*Uhat**2)\n",
    "    f = frac*exp\n",
    "    return(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell determines the probability of a wind speed above 20m/s occuring"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "V_min = 0.0 # minimum wind or tide speed\n",
    "V_max = 10000 # maximum wind or tide speed\n",
    "No_Int = 10000 # assess performance at tshi number of values\n",
    "dV=(V_max-V_min)/(No_Int-1)\n",
    "vc=np.zeros(No_Int)\n",
    "W=np.zeros(No_Int)\n",
    "cum = 0\n",
    "for ic in range(No_Int):\n",
    "    vc[ic]=V_min+(ic+1)*dV # find each w for ic in range(No_Int):\n",
    "    W[ic]=Weibull(vc[ic],vrated)\n",
    "    if vc[ic] >=20:\n",
    "        cum = cum + W[ic]\n",
    "\n",
    "print(cum*100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell defines the function used for energy capture calculations using the Weibull distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def EnergyCap(V_max, V_min, Pws=Pws, T=1, No_Int=100, gen_eff=0.965, gbox_eff=0.985):\n",
    "    P = Pshift(V_min, V_max, No_Int, P=Pws)\n",
    "    Pcum = 0\n",
    "    dV=(V_max-V_min)/(No_Int-1)\n",
    "    vc=np.zeros(No_Int)\n",
    "    hrs = T*365*24\n",
    "    for ic in range(No_Int):\n",
    "        vc[ic]=V_min+(ic+1)*dV\n",
    "        Pu = (P[ic]*gen_eff*gbox_eff)*Weibull(vc[ic])*dV\n",
    "        Pcum = Pcum + Pu\n",
    "    E = hrs*Pcum\n",
    "    return(E)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "energy = EnergyCap(20, 0, No_Int=100)\n",
    "print(energy/1e9) # Print the expected energy capture in GWh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell reads the 2019 wind speed data and takes out only the component of wind speed in the downstream direction, with those +/- 90 degrees of the wind turbine orientation being set to 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lines = []\n",
    "windDir = 247.5\n",
    "upper = windDir + 90\n",
    "lower = windDir - 90\n",
    "speeds = []\n",
    "\n",
    "with open('2019Wind.csv', 'r') as wind:\n",
    "    lines = wind.readlines()\n",
    "for line in lines:\n",
    "    split = line.split(',')\n",
    "    Uwind = float(split[0])\n",
    "    Dir = float(split[1])\n",
    "    if Dir>upper or Dir<lower:\n",
    "        speeds.append(0)\n",
    "    else:\n",
    "        angle = ((Dir-windDir)*np.pi)/180\n",
    "        speeds.append(Uwind*np.cos(angle))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Computes the speeds at hub height"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "speedsHH = []\n",
    "for speed in speeds:\n",
    "    speed = speed*np.exp(A)\n",
    "    speedsHH.append(speed)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell finds the probability of 100 wind speed ranges for the 2019 wind speed data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hist, bins = np.histogram(speedsHH, bins=100)\n",
    "hist = hist/np.sum(hist)\n",
    "bin_width = bins[1]-bins[0]\n",
    "probs = []\n",
    "speedsHH = sorted(speedsHH)\n",
    "for speed in speedsHH:\n",
    "    bin_idx = int((speed-bins[0])/bin_width)\n",
    "    prob = hist[bin_idx-1]\n",
    "    probs.append(prob)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(speedsHH,probs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A function is defined to determine the probability of a given wind speed occuring within a given wind speed dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def probability(speed, data=speedsHH):\n",
    "    hist, bins = np.histogram(data, bins=100)\n",
    "    hist = hist/np.sum(hist)\n",
    "    bin_width = bins[1]-bins[0]\n",
    "    bin_idx = int((speed-bins[0])/bin_width)\n",
    "    prob = hist[bin_idx]\n",
    "    return(prob)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following defines a function to determine the energy capture of one turbine using the realistic probability distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def EnergyCap2(V_max, V_min, Pws=Pws, T=1, No_Int=100, gen_eff=0.965, gbox_eff=0.985, data=speedsHH):\n",
    "    P = Pshift(V_min, V_max, No_Int, P=Pws)\n",
    "    Pcum = 0\n",
    "    dV=(V_max-V_min)/(No_Int-1)\n",
    "    vc=np.zeros(No_Int)\n",
    "    hrs = T*365*24\n",
    "    for ic in range(No_Int):\n",
    "        vc[ic]=V_min+(ic+1)*dV\n",
    "        Pu = (P[ic]*gen_eff*gbox_eff)*probability(vc[ic],data=data)*dV\n",
    "        Pcum = Pcum + Pu\n",
    "    E = hrs*Pcum\n",
    "    return(E/1e9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "EnergyCap2(20, 0, No_Int=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Ct distribution against wind speed is determined, and from this a function is defined to determine the wind speed loss due to wake interaction for a given wind speed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "CtInt = interp1d(V,Cts,fill_value='extrapolate')\n",
    "def deltaV(U):\n",
    "    speedsNEW = []\n",
    "    for i in U:\n",
    "        Ct = np.abs(CtInt(i))\n",
    "        deltaV = i*(1-np.sqrt(1-Ct))*(1/((1+2*0.04*6))**2)\n",
    "        speednew = i-deltaV\n",
    "        speedsNEW.append(speednew)\n",
    "    return(speedsNEW)\n",
    "Ctrated = CtInt(vrated)\n",
    "Thrust = (1/8)*Ctrated*1.223*vrated*vrated*np.pi*D*D\n",
    "print(Thrust)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The speed datasets are modified for each layer of turbines in the array, i.e. turbines with more turbines in front of them will experience lower wind speeds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "speedsHH2 = deltaV(speedsHH)\n",
    "speedsHH3 = deltaV(speedsHH2)\n",
    "speedsHH4 = deltaV(speedsHH3)\n",
    "speedsHH5 = deltaV(speedsHH4)\n",
    "speedsHH6 = deltaV(speedsHH5)\n",
    "speedsHH7 = deltaV(speedsHH6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The energy capture is determined for turbines at each layer within the array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'EnergyCap2' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m~\\AppData\\Local\\Temp/ipykernel_28384/948480364.py\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mEnergy1\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mEnergyCap2\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m20\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mdata\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mspeedsHH\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      2\u001b[0m \u001b[0mEnergy2\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mEnergyCap2\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m20\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mdata\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mspeedsHH2\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[0mEnergy3\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mEnergyCap2\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m20\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mdata\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mspeedsHH3\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[0mEnergy4\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mEnergyCap2\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m20\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mdata\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mspeedsHH4\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[0mEnergy5\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mEnergyCap2\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m20\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mdata\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mspeedsHH5\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'EnergyCap2' is not defined"
     ]
    }
   ],
   "source": [
    "Energy1 = EnergyCap2(20,0,data=speedsHH)\n",
    "Energy2 = EnergyCap2(20,0,data=speedsHH2)\n",
    "Energy3 = EnergyCap2(20,0,data=speedsHH3)\n",
    "Energy4 = EnergyCap2(20,0,data=speedsHH4)\n",
    "Energy5 = EnergyCap2(20,0,data=speedsHH5)\n",
    "Energy6 = EnergyCap2(20,0,data=speedsHH6)\n",
    "Energy7 = EnergyCap2(20,0,data=speedsHH7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "HHs = np.arange(1,8,1)\n",
    "energies = [Energy1,Energy2,Energy3,Energy4,Energy5,Energy6,Energy7]\n",
    "plt.plot(HHs,energies)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The number of turbines at each layer is counted and from this the total energy capture is determined"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "block0=22\n",
    "block1=21\n",
    "block2=20\n",
    "block3=18\n",
    "block4=18\n",
    "block5=11\n",
    "block6=0\n",
    "energytot = block0*Energy1+block1*Energy2+block2*Energy3+block3*Energy4+block4*Energy5+block5*Energy6\n",
    "print(f'The total energy of the array over one year is: {energytot} GWh')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell was used to find the LCOE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "years = 27\n",
    "ran = np.arange(-5,years+2,1)\n",
    "r=0.06\n",
    "DevCost = 120000*1100\n",
    "PlantCost = 600000*1100\n",
    "InstCost = 650000*1100\n",
    "OMCost = 75000*1100\n",
    "DecommCost = 310000*1100\n",
    "towercost = 1240000*110\n",
    "othercost = 100000\n",
    "E=energytot*1000\n",
    "cumcost = 0\n",
    "cumE = 0\n",
    "for year in ran:\n",
    "    if year == 1:\n",
    "        I = towercost  + DevCost + PlantCost + InstCost\n",
    "    elif year == max(ran):\n",
    "        I = DecommCost\n",
    "    else:\n",
    "        I = 0\n",
    "    Mcost = OMCost\n",
    "    comp = (1+r)**year\n",
    "    cumcost = ((I+Mcost)/comp) + cumcost\n",
    "    cumE = E/comp + cumE\n",
    "\n",
    "LCOE = cumcost/cumE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(LCOE)"
   ]
  }
 ],
 "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}