{
 "metadata": {
  "name": "",
  "signature": "sha256:e90b1b15a82d2add026467aa7aa99dbb6e4a09bf9870c7b9bcee8126de0b5def"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Chapter4- Aircraft Gas turbine Engines "
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Ex1-pg133"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import math\n",
      "print(\"Example 4.1\");\n",
      "M0=0.85\n",
      "p0=10000. ##ambient static pressure in Pa\n",
      "pt2=15.88*10**3. ##total pressure at the engine face in Pa\n",
      "gm=1.4 ##gamma\n",
      "pt0=p0*((1.+((gm-1.)*(M0)**2.)/2.)**(gm/(gm-1.)))\n",
      "Pr=pt2/pt0 ##Pr=total pressure recovery\n",
      "ie=((pt2/p0)**((gm-1.)/gm)-1.)/(((gm-1.)/2)*M0**2.) ##inlet adiabatic efficiency.\n",
      "de=-math.log(Pr)\n",
      "print'%s %.3f %s'%(\"(a)The inlet total pressure recovery:\",Pr,\"\")\n",
      "print'%s %.3f %s'%(\"(b)The inlet adiabatic efficiency:\",ie,\"\")\n",
      "print'%s %.4f %s'%(\"(c)The nondimensional entropy rise caused by the inlet:\",de,\"\")\n",
      "\n",
      " "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Example 4.1\n",
        "(a)The inlet total pressure recovery: 0.990 \n",
        "(b)The inlet adiabatic efficiency: 0.978 \n",
        "(c)The nondimensional entropy rise caused by the inlet: 0.0099 \n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Ex2-pg138"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#calculate compressor exit total temperature and adiabatic efficency and compressor shaft power\n",
      "print(\"Example 4.2\")\n",
      "m=50. ##mass flow rate in kg/s\n",
      "ec=0.9 ##compressore polytropic efficiency\n",
      "Tt2=288. ##inlet total temp in K.\n",
      "pt2=100000. ## inlet total pressure in Pa\n",
      "gm=1.4 ##gama\n",
      "cp=1004. ##specific heat in J/kg.K\n",
      "p=35. ##total pressure ratio\n",
      "tr=p**((gm-1.)/(gm*ec)) ##relation between total pressure and temp ratios\n",
      "Tt3=Tt2*tr ##Total exit temp\n",
      "cae=(p**((gm-1.)/gm)-1.)/(tr-1.) ##compressor adiabatic efficiency\n",
      "pc=m*cp*(Tt3-Tt2)/10**6. ## compressor shaft power\n",
      "print'%s %.1f %s'%(\"(a)Compressor exit total temperature in\",Tt3,\" K :\")\n",
      "print'%s %.2f %s'%(\"(b)Compressor adiabatic efficiency:\",cae,\"\")\n",
      "print'%s %.1f %s'%(\"(c)Comprssor shaft power in\",pc,\" MW :\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Example 4.2\n",
        "(a)Compressor exit total temperature in 890.4  K :\n",
        "(b)Compressor adiabatic efficiency: 0.84 \n",
        "(c)Comprssor shaft power in 30.2  MW :\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Ex3-pg142"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import math\n",
      "#calcualte fuel to air ratio and combustor exit temperature\n",
      "print(\"Example 4.3\")\n",
      "Tt3=800.##in K\n",
      "pt3=2*10**6. ## in Pa\n",
      "m=50. ##air mass flow rate in kg/s\n",
      "gm=1.4 ##gamma\n",
      "cp3=1004. ##specific heat at inlet in j/kg.K.\n",
      "Qr=42000. ##heating valuein kJ/kg\n",
      "mf=1. ##fuel flow rate in kg/s\n",
      "be=0.995 ##burner efficiency\n",
      "p=0.96 ##p=pt4/pt3\n",
      "cp4=1156. ##specific heat at exit in J/kg.K\n",
      "f=mf/m ## fuel-to-air ratio\n",
      "Tt4=(((cp3/cp4)*Tt3)+((f*Qr*be*1000.)/cp4))/(1.+f)\n",
      "pt4=p*pt3/10**6.\n",
      "print'%s %.3f %s'%(\"(a)Fuel-to-air ratio :\",f,\"\")\n",
      "print'%s %.1f %s'%(\"b(1) combustor exit total temperature in\",Tt4,\" K:\")\n",
      "print'%s %.2f %s'%(\"b(2)combustor exit total pressure in\",pt4,\" MPa\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Example 4.3\n",
        "(a)Fuel-to-air ratio : 0.020 \n",
        "b(1) combustor exit total temperature in 1390.0  K:\n",
        "b(2)combustor exit total pressure in 1.92  MPa\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Ex4-pg149"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import math\n",
      "#calculate turbine exit temperature and turbine polytropic efficency and turbine exit total pressure and turbine shaft power\n",
      "print(\"Example 4.4\")\n",
      "m=50. ##air mass flow in kg/s\n",
      "mf=1. ## fuel mass flow in kg/s\n",
      "tae=0.88 ##turbine adiabatic efficiency\n",
      "pe=45*10**6 ##shaft power in Watt\n",
      "cp4=1156 ## in J/kg.K\n",
      "Tt4=1390.0197 ## in K\n",
      "pt4=1.92 ##units in MPa\n",
      "cp5=cp4##specific heat\n",
      "mt=m+mf##total mass\n",
      "gm=1.33 ##gamma\n",
      "ht5=cp4*Tt4/1000.-(pe/(mt*1000.)) \n",
      "##print'%s %.1f %s'%(ht5)\n",
      "Tt5=ht5/(cp5/1000.)\n",
      "y=Tt5/Tt4 ##turbine expansion parameter\n",
      "tpe=math.log(y)/math.log(1.-(1.-y)/tae)\n",
      "pr=y**(gm/((gm-1.)*tpe))\n",
      "pt5=pr*pt4*1000. ## turbine total exit pressure\n",
      "pt=mt*cp5*(Tt4-Tt5)/10**6.\n",
      "print'%s %.1f %s'%(\"(a)Turbine exit total temperature in\",Tt5,\" K :\")\n",
      "print'%s %.1f %s'%(\"(b)Turbine polytropic efficiency:\",tpe,\"\")\n",
      "print'%s %.1f %s'%(\"(c)Turbine exit total pressure in \",pt5,\"kPa :\")\n",
      "print'%s %.1f %s'%(\"(d)Turbine shaft power based on turbine expansion delta(Tt) in \",pt,\"MW:\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Example 4.4\n",
        "(a)Turbine exit total temperature in 626.7  K :\n",
        "(b)Turbine polytropic efficiency: 0.8 \n",
        "(c)Turbine exit total pressure in  37.3 kPa :\n",
        "(d)Turbine shaft power based on turbine expansion delta(Tt) in  45.0 MW:\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Ex5-pg150"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import math\n",
      "#calculate the mixed out total temperature \n",
      "print(\"Example 4.5\")\n",
      "mc=0.5 ##mass flow rate of coolant in kg/s\n",
      "mg=50. ##mass flow rate of hot gas in kg/s\n",
      "htg=1850. ## total enthalpy of gas in kJ/kg\n",
      "htc=904. ##total enthalpy of coolant in kJ/kg\n",
      "Cpmixout=1594. ##in j/kg.K\n",
      "##Energy equation between mixed out state and mixed out state and the hot and cold stream solves this problem:\n",
      "Htmixout=(mc*htc+mg*htg)/(mc+mg)\n",
      "Ttmixout=Htmixout/(Cpmixout/1000.)\n",
      "print'%s %.1f %s'%(\"Mixed-out total enthalpy after the nozzle in \",Htmixout,\"kJ/kg :\")\n",
      "print'%s %.1f %s'%(\"Mixed out temperature in\",Ttmixout,\" K :\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Example 4.5\n",
        "Mixed-out total enthalpy after the nozzle in  1840.6 kJ/kg :\n",
        "Mixed out temperature in 1154.7  K :\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Ex6-pg150"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import math\n",
      "#calculate the entropy change across turbine nozzle blade row\n",
      "print(\"Example 4.6\")\n",
      "Cpg=1156. ##in J/kg.K\n",
      "Pt4=1.92 ##in MPa\n",
      "gm=1.33 ##gamma\n",
      "htg=1850. ##from example 4.5 in kJ/kg\n",
      "htc=904. ##from example 4.5 in kJ/kg\n",
      "Cpc=1.04 ##in kJ/kg.K\n",
      "pl=.02 ##total pressure loss ratio\n",
      "Ttmixout=1154.7 ##from example 4.5 in K.\n",
      "##Calculations:\n",
      "Ttg=htg/(Cpg/1000.) ##hotgas total temp in K.\n",
      "Tt4=Ttg ##same as nozzle entrance temp.\n",
      "Ttc=htc/Cpc ##coolant total temp.\n",
      "Ptmixout=(1.-pl)*Pt4 ##mixed-out total temp.\n",
      "##using gibbs equation\n",
      "de=((gm/(gm-1))*math.log((Ttmixout/Tt4)))-math.log(Ptmixout/Pt4)\n",
      "print'%s %.1f %s'%(\"Entropy change across the turbine nozzle blade row:\",de,\"\")\n",
      "print(\"The negative sign of entropy change is due to cooling.\")\n",
      "print(\"*Ans in book is incorrect as Ptmixout is calculated wrong!\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Example 4.6\n",
        "Entropy change across the turbine nozzle blade row: -1.3 \n",
        "The negative sign of entropy change is due to cooling.\n",
        "*Ans in book is incorrect as Ptmixout is calculated wrong!\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Ex7-pg157"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import math\n",
      "#calculate nozzle totalpressure ratio and nozzle area ratio and nozzle exit mach number\n",
      "print(\"Example 4.7\")\n",
      "NPR=10. ##Pressure ratio\n",
      "gm=1.33 ##gamma\n",
      "Cp=1156. ## in J/kg.K\n",
      "ae=0.94 ##adiabatic efficiency\n",
      "tpr=((NPR)**((gm-1.)/gm)-(ae*((NPR)**((gm-1.)/gm)-1.)))**((-1)*(gm/(gm-1.)))\n",
      "print'%s %.1f %s'%(\"(a)Nozzle total pressure ratio:\",tpr,\"\")\n",
      "de=-math.log(tpr) ##entropy rise inadiabatic nozzle\n",
      "##let p=pt9/p9\n",
      "p=tpr*NPR*1 ##p=pt9/p9; p0=p9 foe expanded nozzle\n",
      "M9=((2/(gm-1))*((p)**(((gm-1)/gm))-1))**(1/2.)\n",
      "print'%s %.3f %s'%(\"(c)Nozzle exit Mach no. M9 (perfectly expanded)\",M9,\"\")\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Example 4.7\n",
        "(a)Nozzle total pressure ratio: 0.8 \n",
        "(c)Nozzle exit Mach no. M9 (perfectly expanded) 2.048 \n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Ex10-pg167"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import math\n",
      "#calculate propulsive efficency of turbojet engine \n",
      "print(\"Example4.10\")\n",
      "Vt0=160. ##takeoff velocity in m/s\n",
      "Vt9=1000. ##takeoff velocity in m/s\n",
      "Vc0=800. ##cruise velocity in m/s\n",
      "Vc9=1000. ##cruise velocity in m/s\n",
      "##using approximation: engine propulsive efficiencfy(pe)=2/(1+V9/V0)\n",
      "pet=2./(1.+(Vt9/Vt0)) ##takeoff\n",
      "pec=2./(1.+(Vc9/Vc0)) ##cruise\n",
      "print'%s %.3f %s'%(\"Engine propulsive efficiency while takeoff:\",pet,\"\")\n",
      "print'%s %.3f %s'%(\"Engine propulsive efficiency while cruise:\",pec,\"\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Example4.10\n",
        "Engine propulsive efficiency while takeoff: 0.276 \n",
        "Engine propulsive efficiency while cruise: 0.889 \n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Ex11-pg176"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import math\n",
      "#calculate total pressure and temperature throughout the engine aswell as fuel to air ratio and non-dimensional specific thrust and thurst specifi fuel consumption and thermal and propulsive efficency\n",
      "print(\"Example 4.11\")\n",
      "M0=2.0 ##Mach no.\n",
      "p0=10.##units in kPa\n",
      "T0=228. ##in K\n",
      "gmc=1.4 ##gamma compressor\n",
      "Cpc=1004. ##J/kg.K specific heat of compressor\n",
      "pd=0.88 ##compression ratio of diffuser\n",
      "pc=12. ## compression ratio of compressor\n",
      "ec=0.9 ##adiabatic efficiency of compressor\n",
      "tl=8.    ##enthalpy ratio\n",
      "Qr=42000. ##kJ/kg\n",
      "eb=0.98 ##adiabatic efficiency of burner\n",
      "pb=0.95 ##compression ratio of burner\n",
      "gmt=1.33 ##gamma turbne\n",
      "Cpt=1156. ##J/kg.K specific heat turbine\n",
      "et=0.82 ##adiabatic efficiency of turbine\n",
      "em=0.995 \n",
      "tlAB=11. ##enthalpy ratio of afterburner (AB==AfterBurner)\n",
      "QrAB=42000. ##kJ/kg\n",
      "eAB=0.98\n",
      "pAB=0.93\n",
      "gmAB=1.3 ## gama AB\n",
      "CpAB=1243. ##J/kg.K\n",
      "pn=0.93\n",
      "a0=((gmc-1.)*Cpc*T0)**(1/2.)\n",
      "V0=M0*a0\n",
      "pt0=p0*(1.+(((gmc-1.)*(M0)**2.)/2.))**(gmc/(gmc-1.)) ##total flight pressure\n",
      "Tt0=T0*(1.+(((gmc-1.)*(M0)**2)/2.)) ##total flight temp\n",
      "Tt2=Tt0 ##Adiabatic inlets\n",
      "pt2=pt0*pd ## in kPa\n",
      "pt3=pt2*pc ##compressor exit total pressure\n",
      "k2=((gmc-1.)/(gmc*ec))\n",
      "##print'%s %.1f %s'%(k2)\n",
      "tc=pc**k2 ##relation between temp and pressure ratios\n",
      "##print'%s %.1f %s'%(tc)\n",
      "Tt3=Tt2*tc ##total temp at compressor exit\n",
      "Tt4=Cpc*T0*tl/Cpt ##combustor exit total temp.\n",
      "pt4=pt3*pb ##combustor exit pressure\n",
      "f=(Cpt*Tt4-Cpc*Tt3)/(Qr*eb*1000.-Cpt*Tt4) ##fuel-to-air ratio in burner\n",
      "##print'%s %.1f %s'%(f)\n",
      "Tt5=Tt4-(Cpc*((Tt3-Tt2)/(Cpt*em*(1.+f)))) ## turbine exit total temp\n",
      "tt=Tt5/Tt4 ##temp ratio in turbine\n",
      "pt=tt**(gmt/(et*(gmt-1.)))\n",
      "pt5=pt4*pt ##in kPa\n",
      "pt7=pt5*pAB\n",
      "Tt7=Cpc*T0*tlAB/CpAB ##afterburner exit\n",
      "fAB=(1+f)*((CpAB*Tt7)-(Cpt*Tt5))/((QrAB*eAB*1000.)-(CpAB*Tt7))\n",
      "##print'%s %.1f %s'%(fAB)\n",
      "pt9=pt7*pn ##in kPA\n",
      "Tt9=Tt7 ##adiabatic flow in nozzle\n",
      "p9=p0\n",
      "M9=((2./(gmAB-1.))*((pt9/p9)**(((gmAB-1)/gmAB))-1))**(1/2.) ##nozzle exit\n",
      "##print'%s %.1f %s'%(M9)\n",
      "T9=Tt9/(1.+((gmAB-1)*(M9)**2)/2.)\n",
      "a9=((gmAB-1.)*CpAB*T9)**(1/2.)\n",
      "##print'%s %.1f %s'%(a9)\n",
      "V9=M9*a9\n",
      "##Performance parameters:\n",
      "st=(1.+f+fAB)*V9-V0 ##st=Fn/m0; specific thrust when nozzle is perfectly expanded\n",
      "ndst=((1.+f+fAB)*V9/a0)-M0 ##ndst=Fn/m0*ao ; nondimensional specific thrust\n",
      "TSFC=((f+fAB)/st)*10**6. ##units mg/s/N\n",
      "eth=(((1.+f+fAB)*((V9)**2)/2.)-((V0)**2.)/2.)/(f*Qr*1000.+fAB*QrAB*1000.) ##cycle thermal efficiency\n",
      "ep=st*V0/(((1.+f+fAB)*(((V9)**2)/2.))-((V0)**2)/2.) ##propulsive efficiency exact\n",
      "epa=2./(1.+V9/V0) ##approx\n",
      "print(\"a(1)Total temperatures across the engine in\"\" K:\")\n",
      "print'%s %.1f %s'%(\"Flight total temperaure:\",Tt0,\" \")\n",
      "\n",
      "print'%s %.1f %s'%(\"Toal temperature at compressor inlet:\",Tt2,\"\")\n",
      "print'%s %.1f %s'%(\"Total temperature at compressor exit:\",Tt3,\" \")\n",
      "print'%s %.1f %s'%(\"Total temperature at burner exit:\",Tt4,\"\")\n",
      "print'%s %.1f %s'%(\"Total temperature at turbine exit:\",Tt5,\"\")\n",
      "print'%s %.1f %s'%(\"Total temperature at afterburner exit:\",Tt7,\"\")\n",
      "print'%s %.1f %s'%(\"Total temperature at nozzle exit:\",T9,\"\")\n",
      "print'%s %.1f %s'%(\"Nozzle exit static temperature:\",T9,\"\")\n",
      "print(\"a(2)Total pressures across the engine in kPa:\")\n",
      "print'%s %.1f %s'%(\"Flight total pressure:\",pt0,\"\")\n",
      "\n",
      "print'%s %.1f %s'%(\"Toal pressure at compressor inlet:\",pt2,\"\")\n",
      "print'%s %.1f %s'%(\"Total pressure at compressor exit:\",pt3,\" \")\n",
      "print'%s %.1f %s'%(\"Total pressure at burner exit:\",pt4,\"\")\n",
      "print'%s %.1f %s'%(\"Total pressure at turbine exit:\",pt5,\"\")\n",
      "print'%s %.1f %s'%(\"Total pressure at afterburner exit:\",pt7,\"\")\n",
      "print'%s %.1f %s'%(\"Total pressure at nozzle exit:\",pt9,\"\")\n",
      "print'%s %.1f %s'%(\"Nozzle exit static pressure:\",p9,\"\")\n",
      "print'%s %.1f %s'%(\"(b)Nondimensional specific thrust:\",ndst,\"\")\n",
      "print'%s %.1f %s'%(\"(c)Thrust specific fuel consumption TSFC \",TSFC,\"(in mg/s/N):\")\n",
      "print'%s %.1f %s'%(\"d(1)Themal efficiency:\",eth,\"\")\n",
      "print'%s %.1f %s'%(\"d(2)Exact propulsive efficiency:\",ep,\"\")\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Example 4.11\n",
        "a(1)Total temperatures across the engine in K:\n",
        "Flight total temperaure: 410.4  \n",
        "Toal temperature at compressor inlet: 410.4 \n",
        "Total temperature at compressor exit: 903.2  \n",
        "Total temperature at burner exit: 1584.2 \n",
        "Total temperature at turbine exit: 1163.9 \n",
        "Total temperature at afterburner exit: 2025.8 \n",
        "Total temperature at nozzle exit: 1085.8 \n",
        "Nozzle exit static temperature: 1085.8 \n",
        "a(2)Total pressures across the engine in kPa:\n",
        "Flight total pressure: 78.2 \n",
        "Toal pressure at compressor inlet: 68.9 \n",
        "Total pressure at compressor exit: 826.3  \n",
        "Total pressure at burner exit: 784.9 \n",
        "Total pressure at turbine exit: 172.5 \n",
        "Total pressure at afterburner exit: 160.4 \n",
        "Total pressure at nozzle exit: 149.2 \n",
        "Nozzle exit static pressure: 10.0 \n",
        "(b)Nondimensional specific thrust: 3.3 \n",
        "(c)Thrust specific fuel consumption TSFC  54.2 (in mg/s/N):\n",
        "d(1)Themal efficiency: 0.5 \n",
        "d(2)Exact propulsive efficiency: 0.6 \n"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Ex13-pg187"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#calculate nozzle exit static pressure and actual and effective nozzle exit velocties and ratio of fan to core thrust and non-dimensional specific thrust and TSFC and all engine effeciences \n",
      "print(\"Example4.13\")\n",
      "M0=0.88 ##Mach no.\n",
      "p0=15 ## pressure in kPa\n",
      "T0=233 ##temperatue in K\n",
      "gmc=1.4 ##gamma compressor\n",
      "Cpc=1004 ##specific heat of compressor in J/kg.K\n",
      "pd=0.995 ## pressure compression ratio of diffuser\n",
      "pf=1.6 ##pressure compression ratio of fan\n",
      "ef=0.9 ##fan efficiency\n",
      "alfa=8\n",
      "pfn=0.95 ##compression ratio of convergent fan nozzle\n",
      "pc=40 ##compression ratio of compressor\n",
      "ec=0.9 ##compressor efficiency\n",
      "tl=8 ##temp. ratio\n",
      "Cpt=1152 ##in J/kg.K of turbine\n",
      "gmt=1.33 ##gamma turbine\n",
      "Qr=42000000 ##in J/kg\n",
      "pb=0.95 ##burner compression ratio\n",
      "eb=0.992 ##burner efficiency\n",
      "em=0.95\n",
      "et=0.85\n",
      "pn=0.98 ##primary nozzle\n",
      "a0=((gmc-1)*Cpc*T0)**(1/2.);\n",
      "V0=M0*a0;\n",
      "pt0=p0*(1.+((gmc-1.)*(M0)**2)/2.)**(gmc/(gmc-1.))\n",
      "Tt0=T0*(1+((gmc-1.)*(M0)**2)/2.)\n",
      "Tt2=Tt0\n",
      "pt2=pt0*pd\n",
      "##fan stream:\n",
      "pt13=pt2*pf\n",
      "tf=pf**((gmc-1.)/(ef*gmc))\n",
      "Tt13=Tt2*tf\n",
      "pt19=pt13*pfn\n",
      "p19=pt19/(1.+(gmc-1)/2.)**(gmc/(gmc-1.))\n",
      "M19=1.\n",
      "T19=Tt13/1.2\n",
      "a19=((gmc-1)*Cpc*T19)**(1/2.)\n",
      "V19=a19\n",
      "##V19eff=V19+((gmc*p19)/r19)*((1-p0/p19)/(gmc*V19)) i.e V19+a19**2\n",
      "V19eff=V19+(a19**2.)*((1.-p0/p19)/(gmc*V19))\n",
      "##Core stream\n",
      "pt3=pt2*pc\n",
      "tc=pc**((gmc-1.)/(ec*gmc))\n",
      "##print'%s %.1f %s'%(tc)\n",
      "Tt3=Tt2*tc\n",
      "pt4=pt3*pb\n",
      "Tt4=Cpc*T0*tl/Cpt\n",
      "##print'%s %.1f %s'%(Tt4)\n",
      "f=(Cpt*Tt4-Cpc*Tt3)/(Qr*eb-Cpt*Tt4)\n",
      "##print'%s %.1f %s'%(f)\n",
      "Tt5=Tt4-((Cpc*(Tt3-Tt2)+alfa*Cpc*(Tt13-Tt2)))/((1+f)*Cpt*em)\n",
      "##print'%s %.1f %s'%(Tt5)\n",
      "tt=Tt5/Tt4\n",
      "pt=tt**(gmt/(et*(gmt-1)))\n",
      "pt5=pt4*pt\n",
      "pt9=pt5*pn\n",
      "p9=pt9/((gmt+1)/2.)**(gmt/(gmt-1.))\n",
      "M9=1.\n",
      "T9=Tt5/((gmt+1)/2.)\n",
      "a9=((gmt-1)*Cpt*T9)**(1/2.)\n",
      "V9=a9\n",
      "V9eff=V9+(((a9)**2)*(1-(p0/p9)))/(gmt*V9)\n",
      "ndsft=alfa*(V19eff-V0)/((1+alfa)*a0)\n",
      "ndsct=((1+f)*V9eff-V0)/((1+alfa)*a0)\n",
      "ndst=ndsft+ndsct\n",
      "rfct=ndsft/ndsct\n",
      "fc=ndsft*100./(ndsft+ndsct)\n",
      "cc=ndsct*100./(ndsft+ndsct)\n",
      "TSFC=f/((1.+alfa)*a0*(ndsft+ndsct))*10**6.\n",
      "eth=(alfa*V19eff**2+(1+f)*V9eff**2-(1+alfa)*V0**2.)/(2.*f*Qr)\n",
      "ep=(2.*(ndsft+ndsct)*(1+alfa)*a0*V0)/(alfa*V19eff**2.+(1.+f)*V9eff**2.-(1.+alfa)*V0**2.)\n",
      "eo=eth*ep\n",
      "##Pressures\n",
      "print(\"a(1)Total pressures throughout the engine in kPa:\")\n",
      "print'%s %.1f %s'%(\"Total pressure of flight:\",pt0,\"\")\n",
      "print'%s %.1f %s'%(\"Total pressure at engine face:\",pt2,\"\")\n",
      "print'%s %.1f %s'%(\"Total pressure at fan exit:\",pt13,\"\")\n",
      "\n",
      "print'%s %.1f %s'%(\"Static pressure at nozzle exit:\",p19,\"\")\n",
      "print'%s %.1f %s'%(\"Total pressure at compressor exit:\",pt3,\"\")\n",
      "print'%s %.1f %s'%(\"Total pressure at burner exit:\",pt4,\"\")\n",
      "print'%s %.1f %s'%(\"Total pressure at turbine exit:\",pt5,\"\")\n",
      "print'%s %.1f %s'%(\"Total pressure at nozzle exit:\",pt9,\"\")\n",
      "\n",
      "##Temperatures\n",
      "print(\"a(2)Total temperatures across the engine in K:\")\n",
      "print'%s %.1f %s'%(\"Total temperature of flight:\",Tt0,\"\")\n",
      "print'%s %.1f %s'%(\"Total temperature at engine face:\",Tt2,\"\") ##Tt0=Tt2, since adiabatic!\n",
      "print'%s %.1f %s'%(\"Total temperature at fan exit:\",Tt13,\"\")\n",
      "print'%s %.1f %s'%(\"Static  temperature at fan nozzle exit:\",T19,\"\")\n",
      "print'%s %.1f %s'%(\"Total temperature at compressor exit:\",Tt3,\"\")\n",
      "print'%s %.1f %s'%(\"Total temperature at burner exit:\",Tt4,\"\")\n",
      "print'%s %.1f %s'%(\"Total temperature at turbine exit:\",Tt5,\"\")\n",
      "print'%s %.1f %s'%(\"Static temperature at nozzle exit:\",T9,\"\")\n",
      "print'%s %.1f %s'%(\"(b{1})Total pressure at fan nozzle exit:\",pt19,\"\")\n",
      "print'%s %.1f %s'%(\"(b{2})Static pressure at nozzle exit:\",p9,\"\")\n",
      "\n",
      "\n",
      "##Remaining results\n",
      "print'%s %.1f %s'%(\"(c{1}Actual fan nozzle exit velocity in\",V19,\" m/s:)\")\n",
      "print'%s %.1f %s'%(\"(c{2}Effective fan nozzle exit velocity in\",V9eff,\" m/s:)\")\n",
      "print'%s %.1f %s'%(\"(c{3})Actual core nozzle exit velocity in\",V9,\" m/s:\")\n",
      "print'%s %.1f %s'%(\"(c{4})Effective nozzle exit velocity in \",V9eff,\"m/s:\")\n",
      "print'%s %.1f %s'%(\"(d)Ratio of fan-tocore thrust:\",rfct,\"\")\n",
      "print'%s %.1f %s'%(\"(e)Nondimensional specific thrust:\",ndst,\"\")\n",
      "print'%s %.1f %s'%(\"(f)TSFC in mg/s/N:\",TSFC,\"\")\n",
      "print(\"(g)Engine efficiencies:\")\n",
      "print'%s %.1f %s'%(\"Thermal efficiency:\",eth,\"\")\n",
      "print'%s %.1f %s'%(\"Propulsion effciency:\",ep,\"\")\n",
      "print'%s %.1f %s'%(\"Overall efficiency:\",eo,\"\")\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Example4.13\n",
        "a(1)Total pressures throughout the engine in kPa:\n",
        "Total pressure of flight: 24.8 \n",
        "Total pressure at engine face: 24.7 \n",
        "Total pressure at fan exit: 39.5 \n",
        "Static pressure at nozzle exit: 19.8 \n",
        "Total pressure at compressor exit: 988.2 \n",
        "Total pressure at burner exit: 938.8 \n",
        "Total pressure at turbine exit: 28.7 \n",
        "Total pressure at nozzle exit: 28.1 \n",
        "a(2)Total temperatures across the engine in K:\n",
        "Total temperature of flight: 269.1 \n",
        "Total temperature at engine face: 269.1 \n",
        "Total temperature at fan exit: 312.4 \n",
        "Static  temperature at fan nozzle exit: 260.3 \n",
        "Total temperature at compressor exit: 867.9 \n",
        "Total temperature at burner exit: 1624.0 \n",
        "Total temperature at turbine exit: 778.1 \n",
        "Static temperature at nozzle exit: 667.9 \n",
        "(b{1})Total pressure at fan nozzle exit: 37.6 \n",
        "(b{2})Static pressure at nozzle exit: 15.2 \n",
        "(c{1}Actual fan nozzle exit velocity in 323.3  m/s:)\n",
        "(c{2}Effective fan nozzle exit velocity in 508.4  m/s:)\n",
        "(c{3})Actual core nozzle exit velocity in 503.9  m/s:\n",
        "(c{4})Effective nozzle exit velocity in  508.4 m/s:\n",
        "(d)Ratio of fan-tocore thrust: 3.5 \n",
        "(e)Nondimensional specific thrust: 0.4 \n",
        "(f)TSFC in mg/s/N: 22.1 \n",
        "(g)Engine efficiencies:\n",
        "Thermal efficiency: 0.4 \n",
        "Propulsion effciency: 0.8 \n",
        "Overall efficiency: 0.3 \n"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Ex15-pg199"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import math\n",
      "print(\"Example 4.15\")\n",
      "M0=2. ##Mach no.\n",
      "p0=10. ## in kPa\n",
      "T0=223. ##in K\n",
      "##the engine inlet total pressure loss is characterized by \n",
      "pd=0.9\n",
      "##The fan pressure ratio is\n",
      "pf=1.9\n",
      "##and polytropic efficiency of the fan is\n",
      "ef=0.9\n",
      "##The flow in the fan duct suffers 1% total pressure loss i.e.\n",
      "pfd=0.99\n",
      "##The compressor pressure ratio and polytropic efficiency are \n",
      "pc=13.\n",
      "ec=0.9 ##respectively\n",
      "##The combustor exit temperature is \n",
      "Tt4=1600. ##in K\n",
      "Qr=42000000. ##fuel heating value in J/kg\n",
      "pb=0.95 ##total pressure ratio\n",
      "eb=0.98 ##burner efficiency\n",
      "et=0.8 ##turbine polytropic efficiency\n",
      "em=0.95 ##mechanical efficiency of turbine\n",
      "M5=0.5 ##Mach no at turbine exit\n",
      "pmf=0.98 ##total pressure loss due to friction in mixer\n",
      "Tt7=2000. ##afterburner total temp in K\n",
      "QrAB=42000000. ##in J/kg\n",
      "pABon=0.92\n",
      "eAB=0.98\n",
      "pn=0.95 ##total pressure ratio at nozzle\n",
      "p=3.8 ##p=p9/p0\n",
      "gmc=1.4 ##gamma compressor\n",
      "Cpc=1004. ##specofic heat compressor in J/kg.K\n",
      "gmt=1.33 ##gamma turbine\n",
      "Cpt=1152. ##turbine\n",
      "gmAB=1.3 ##afterburner\n",
      "CpAB=1241. ##afterburner\n",
      "pt0=p0*(1.+((gmc-1.)*(M0)**2)/2.)**(gmc/(gmc-1.))\n",
      "Tt0=T0*(1.+((gmc-1.)*(M0)**2)/2.)\n",
      "pr=pt0/p0\n",
      "tr=Tt0/T0\n",
      "pt=pfd*pf/(pb*pc)\n",
      "a0=((gmc-1.)*Cpc*T0)**(1/2.);\n",
      "V0=a0*M0\n",
      "Tt2=Tt0\n",
      "pt2=pt0*pd\n",
      "pt13=pt2*pf\n",
      "tf=pf**((gmc-1.)/(ec*gmc))\n",
      "##print'%s %.1f %s'%(tf)\n",
      "Tt13=Tt0*tf\n",
      "Tt15=Tt13 ##adiabatic\n",
      "pt15=pt13*pfd\n",
      "pt3=pt2*pc\n",
      "tc=pc**((gmc-1.)/(ec*gmc))\n",
      "Tt3=Tt2*tc\n",
      "pt4=pt3*pb\n",
      "f=(Cpt*Tt4-Cpc*Tt3)/(Qr*eb-Cpt*Tt4)\n",
      "##print'%s %.1f %s'%(f)\n",
      "pt5=pt15 ##assumption\n",
      "pt=(pfd*pf)/(pb*pc)\n",
      "##print'%s %.1f %s'%(pt)\n",
      "tt=pt**(et*(gmt-1.)/(gmt))\n",
      "##print'%s %.1f %s'%(tt)\n",
      "Tt5=Tt4*tt\n",
      "tl=(Cpt*Tt4)/(Cpc*T0)\n",
      "tr=(1.+((gmc-1.)*(M0**2)/2.))\n",
      "alfa=((em*(1.+f)*tl*(1.-tt))-(tr*(tc-1.)))/(tr*(tf-1.))\n",
      "ht6M=Cpc*T0*((1.+f)*tt*tl+alfa*tf*tr)/(1.+alfa+f) ## mixed-out total enthalpy in J/kg\n",
      "Cp6M=(((1.+f)/alfa)*Cpt+Cpc)/(((1.+f)/alfa)+1.)\n",
      "gm6M=(((1+f)/alfa)*Cpt+Cpc)/(((1+f)/alfa)*(Cpt/gmt)+(Cpc/gmc))\n",
      "M15=((2./(gmc-1.))*((((1.+((gmt-1.)*(M5**2)/2.))**(gmt/(gmt-1.)))**((gmc-1.)/gmc))-1.))**(1/2.)\n",
      "T15=Tt15/(1.+((gmc-1.)*(M15)**2)/2.)\n",
      "p15=pt15/(1.+((gmc-1.)*(M15)**2)/2.)**(gmc/(gmc-1.))\n",
      "T5=Tt5/(1.+((gmt-1.)*(M5)**2)/2.)\n",
      "p5=pt5/(1.+((gmt-1.)*(M5)**2)/2.)**(gmt/(gmt-1.))\n",
      "a15=((gm6M-1.)*Cp6M*T15)**(1/2.)\n",
      "a5=((gm6M-1.)*Cp6M*T5)**(1/2.)\n",
      "A=((alfa/(1.+f))*(gmt/gmc)*((T15/T5)**(1/2.))*(M5/M15))\n",
      "C1=((1.+gmt*M5**2.)+(A*(1.+gmc*M15**2.)))/(1.+A)\n",
      "Tt6M=ht6M/Cp6M\n",
      "C2=((gmt/gm6M)*(M5/a5)+(gmc/gm6M)*(M15*A/a15))*(((gm6M-1.)*Cp6M*(Tt6M))**(1/2.))/(1.+A)\n",
      "C=(C1/C2)**2.\n",
      "M6M=((C-2*gm6M-((C-2.*gm6M)**2-4.*(gm6M**2.-(C*(gm6M-1))/2.))**(1/2.))/(2*(gm6M)**2.-C*(gm6M-1.)))**(1/2.)\n",
      "p6M=p5*(C1/(1.+gm6M*(M6M)**2.))\n",
      "pt6Mi=131.23\n",
      "pmi=0.9907\n",
      "pM=0.9709\n",
      "pt6M=pt6Mi*pmf\n",
      "Tt7=2000.\n",
      "pABon=0.92\n",
      "pt7=118.32\n",
      "fAB=(CpAB*Tt7-ht6M)/(QrAB*eAB-CpAB*Tt7)\n",
      "pt9=pt7*pn\n",
      "p9=p0*p\n",
      "M9=1.377\n",
      "T9=1557.2\n",
      "a9=761.4\n",
      "V9=a9*M9\n",
      "V9eff=V9+a9**2.*(1.-p0/p9)/(gmAB*V9)\n",
      "ndst=((1+alfa+f+fAB)/(1.+alfa))*(V9eff/a0)-M0\n",
      "TSFC=((f+fAB)/((1+alfa)*a0))*10**6./(ndst)\n",
      "eth=(((1+alfa+f+fAB)*((V9eff)**2.))-((1+alfa)*V0**2.))/(2.*(f*Qr+fAB*QrAB))\n",
      "ep=(2.*ndst*V0*a0*(1.+alfa))/((1.+alfa+f+fAB)*V9eff**2-(1.+alfa)*V0**2)\n",
      "e0=ep*eth\n",
      "print(\"a(1)Total pressures throughout the engine in kPa:\")\n",
      "print'%s %.1f %s'%(\"Total pressure of flight:\",pt0,\"\")\n",
      "print'%s %.1f %s'%(\"Total pressure at engine face:\",pt2,\"\")\n",
      "print'%s %.1f %s'%(\"Total pressure at fan exit:\",pt15,\"\")\n",
      "##print'%s %.1f %s'%(p19,\"Static pressure at nozzle exit:\")\n",
      "print'%s %.1f %s'%(\"Total pressure at compressor exit:\",pt3,\"\")\n",
      "print'%s %.1f %s'%(\"Total pressure at burner exit:\",pt4,\"\")\n",
      "print'%s %.1f %s'%(\"Total pressure at turbine exit:\",pt5,\"\")\n",
      "print'%s %.1f %s'%(\"Total pressure at nozzle exit:\",pt9,\"\")\n",
      "\n",
      "\n",
      "print(\"a(2)Total temperatures across the engine in K:\")\n",
      "print'%s %.1f %s'%(\"Total temperature of flight:\",Tt0,\"\")\n",
      "print'%s %.1f %s'%(\"Total temperature at engine face:\",Tt2,\"\") ##Tt0=Tt2, since adiabatic!\n",
      "print'%s %.1f %s'%(\"Total temperature at fan exit:\",Tt13,\"\")\n",
      "print'%s %.1f %s'%(\"Total temperature at fan duct :\",Tt15,\"\")\n",
      "print'%s %.1f %s'%(\"Total temperature at compressor exit:\",Tt3,\"\")\n",
      "print'%s %.1f %s'%(\"Total temperature at burner exit:\",Tt4,\"\")\n",
      "print'%s %.1f %s'%(\"Total temperature at turbine exit:\",Tt5,\"\")\n",
      "print'%s %.1f %s'%(\"a(3)Fan bypass ratio :\",alfa,\"\")\n",
      "print'%s %.3f %s'%(\"a(4)fuel-to-air ratio in primary :\",f,\"\")\n",
      "print'%s %.3f %s'%(\"a(5)fuel-to-air ratio in afterburner :\",fAB,\"\")\n",
      "print'%s %.1f %s'%(\"b(1)TSFC in mg/s/N :\",TSFC,\"\")\n",
      "print'%s %.1f %s'%(\"b(2)Non-dimensional specific thrust :\",ndst,\"\")\n",
      "print'%s %.1f %s'%(\"b(3)Propulsive efficiency :\",ep,\"\")\n",
      "print'%s %.1f %s'%(\"b(4)Thermal efficiency :\",eth,\"\")\n",
      "print'%s %.1f %s'%(\"b(5)Overall efficiency :\",e0,\"\")\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Example 4.15\n",
        "a(1)Total pressures throughout the engine in kPa:\n",
        "Total pressure of flight: 78.2 \n",
        "Total pressure at engine face: 70.4 \n",
        "Total pressure at fan exit: 132.5 \n",
        "Total pressure at compressor exit: 915.5 \n",
        "Total pressure at burner exit: 869.7 \n",
        "Total pressure at turbine exit: 132.5 \n",
        "Total pressure at nozzle exit: 112.4 \n",
        "a(2)Total temperatures across the engine in K:\n",
        "Total temperature of flight: 401.4 \n",
        "Total temperature at engine face: 401.4 \n",
        "Total temperature at fan exit: 492.1 \n",
        "Total temperature at fan duct : 492.1 \n",
        "Total temperature at compressor exit: 906.2 \n",
        "Total temperature at burner exit: 1600.0 \n",
        "Total temperature at turbine exit: 1101.3 \n",
        "a(3)Fan bypass ratio : 0.6 \n",
        "a(4)fuel-to-air ratio in primary : 0.024 \n",
        "a(5)fuel-to-air ratio in afterburner : 0.039 \n",
        "b(1)TSFC in mg/s/N : 48.5 \n",
        "b(2)Non-dimensional specific thrust : 2.7 \n",
        "b(3)Propulsive efficiency : 0.6 \n",
        "b(4)Thermal efficiency : 0.5 \n",
        "b(5)Overall efficiency : 0.3 \n"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Example 4.17 - pg "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import math\n",
      "print(\"Example 4.17\")\n",
      "M0=0.7 ##Mach no.\n",
      "T0=228 ## in K\n",
      "p0=16 ##kPa\n",
      "eprop=0.85 ## prop efficiency\n",
      "m=10. ##Kg/s\n",
      "pd=0.98 ##diffuser pressure ratio\n",
      "pc=30. ##compressor pressurer ratio\n",
      "ec=0.92 ##thermal efficiency of compressor\n",
      "Tt4=1600. ##in K\n",
      "Qr=42000000. ##in kJ/kg\n",
      "eb=0.99 ##thermal efficiency of burner\n",
      "pb=0.96 ##burner pressure ratio\n",
      "etHPT=0.82\n",
      "emHPT=0.99\n",
      "alfa=0.85 \n",
      "emLPT=0.99\n",
      "eLPT=0.88\n",
      "egb=0.995\n",
      "en=0.95\n",
      "gmc=1.4 ##gamma of compressor\n",
      "Cpc=1004. ## in J/kg.K\n",
      "gmt=1.33 ##gamma of turbine\n",
      "Cpt=1152. ## in J/kg.K\n",
      "Tt0=T0*(1.+((gmc-1.)*(M0)**2)/2.)\n",
      "pt0=p0*(1.+((gmc-1.)*(M0)**2)/2.)**(gmc/(gmc-1.))\n",
      "a0=((gmc-1.)*Cpc*T0)**(1/2.);\n",
      "V0=a0*M0\n",
      "pt2=pt0*pd\n",
      "Tt2=Tt0 ##Adiabatic\n",
      "pt3=pt2*pc\n",
      "tc=pc**((gmc-1.)/(ec*gmc))\n",
      "Tt3=Tt2*tc\n",
      "f=(Cpt*Tt4-Cpc*Tt3)/(Qr*eb-Cpt*Tt4)\n",
      "pt4=pt3*pb\n",
      "ht45=Cpt*Tt4-(Cpc*Tt3-Cpc*Tt2)/((1.+f)*emHPT)\n",
      "Tt45=ht45/Cpt\n",
      "pt45=pt4*(Tt45/Tt4)**(gmt/((gmt-1.)*etHPT))\n",
      "m9=(1.+f)*m\n",
      "sp=(1.+f)*m*eLPT*alfa*ht45*(1.-(p0/pt45)**((gmt-1.)/gmt))/10**6.\n",
      "Tt5=(ht45-sp*10**6./((1.+f)*m))/Cpt\n",
      "tt=Tt5/Tt45\n",
      "et=math.log(tt)/(math.log(1-((1-tt)/eLPT)))\n",
      "pt=tt**(gmt/(et*(gmt-1.)))\n",
      "pt5=pt45*pt\n",
      "p9=p0 ##assumption\n",
      "pi=p9/pt5\n",
      "ti=pi**((gmt-1.)/gmt)\n",
      "T9i=Tt5*ti\n",
      "T9=Tt5-en*(Tt5-T9i)\n",
      "V9=(2.*Cpt*(Tt5-T9))**(1/2.)\n",
      "Fprop=eprop*egb*emLPT*sp*10**3/V0\n",
      "a9=((gmt-1.)*Cpt*T9)**(1/2.)\n",
      "M9=V9/a9\n",
      "pt9=p9*(1.+((gmt-1)*M9**2.)/2.)**(gmt/(gmt-1.))\n",
      "pn=pt9/pt5\n",
      "Fncore=m*((1.+f)*V9-V0)/1000.\n",
      "spp=egb*emLPT*sp\n",
      "Ft=Fprop+Fncore\n",
      "mp=((m9*V9**2.)/2.-m*(V0**2.)/2.)/10**3.\n",
      "mf=m9-m\n",
      "PSFC=mf*10**6./((spp*10**3.)+mp)\n",
      "TSFC=mf*10**3./(Ft)\n",
      "eth=(spp*10**3.+mp)*10**3./(mf*Qr)\n",
      "ep=(Ft*V0)/(spp*10**3+mp)\n",
      "eo=eth*ep\n",
      "print(\"a(1)Total pressures throughout the engine in kPa:\")\n",
      "print'%s %.1f %s'%(\"Total pressure of flight:\",pt0,\"\")\n",
      "print'%s %.1f %s'%(\"Total pressure at engine face:\",pt2,\"\")\n",
      "##print'%s %.1f %s'%(p19,\"Static pressure at nozzle exit:\")\n",
      "print'%s %.1f %s'%(\"Total pressure at compressor exit:\",pt3,\"\")\n",
      "print'%s %.1f %s'%(\"Total pressure at burner exit:\",pt4,\"\")\n",
      "print'%s %.1f %s'%(\"Total pressure across HPT :\",pt45,\"\")\n",
      "print'%s %.1f %s'%(\"Total pressure at turbine exit:\",pt5,\"\")\n",
      "print'%s %.1f %s'%(\"Total pressure at nozzle exit:\",pt9,\"\")\n",
      "\n",
      "print(\"a(2)Total temperatures across the engine in K:\")\n",
      "print'%s %.1f %s'%(\"Total temperature of flight:\",Tt0,\"\")\n",
      "print'%s %.1f %s'%(\"Total temperature at engine face:\",Tt2,\"\") ##Tt0=Tt2, since adiabatic!\n",
      "print'%s %.1f %s'%(\"Total temperature at compressor exit:\",Tt3,\"\")\n",
      "print'%s %.1f %s'%(\"Total temperature at burner exit:\",Tt4,\"\")\n",
      "print'%s %.1f %s'%(\"Total temperature across HPT :\",Tt45,\"\")\n",
      "print'%s %.1f %s'%(\"Total temperature at turbine exit:\",Tt5,\"\")\n",
      "print'%s %.1f %s'%(\"a(3)fuel-to-air ratio in burner :\",f,\"\")\n",
      "print'%s %.1f %s'%(\"(b)Engine core thrust in kN :\",Fncore,\"\")\n",
      "print'%s %.1f %s'%(\"(c)Propeller thrust in kN :\",Fprop,\"\")\n",
      "print'%s %.1f %s'%(\"(d)Power-specific fuel consumption in\",PSFC,\" mg/s/kW :\")\n",
      "print'%s %.1f %s'%(\"(e)TSFC in\",TSFC,\" mg/s/N :\")\n",
      "print'%s %.3f %s'%(\"f(1)Propulsive efficiency :\",ep,\"\")\n",
      "print'%s %.3f %s'%(\"f(2)Thermal efficiency :\",eth,\"\")\n",
      "print'%s %.3f %s'%(\"(g)Overall efficiency :\",eo,\"\")\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Example 4.17\n",
        "a(1)Total pressures throughout the engine in kPa:\n",
        "Total pressure of flight: 22.2 \n",
        "Total pressure at engine face: 21.7 \n",
        "Total pressure at compressor exit: 652.5 \n",
        "Total pressure at burner exit: 626.4 \n",
        "Total pressure across HPT : 151.1 \n",
        "Total pressure at turbine exit: 24.5 \n",
        "Total pressure at nozzle exit: 24.0 \n",
        "a(2)Total temperatures across the engine in K:\n",
        "Total temperature of flight: 250.3 \n",
        "Total temperature at engine face: 250.3 \n",
        "Total temperature at compressor exit: 719.9 \n",
        "Total temperature at burner exit: 1600.0 \n",
        "Total temperature across HPT : 1198.0 \n",
        "Total temperature at turbine exit: 815.2 \n",
        "a(3)fuel-to-air ratio in burner : 0.0 \n",
        "(b)Engine core thrust in kN : 2.2 \n",
        "(c)Propeller thrust in kN : 17.9 \n",
        "(d)Power-specific fuel consumption in 54.6  mg/s/kW :\n",
        "(e)TSFC in 14.0  mg/s/N :\n",
        "f(1)Propulsive efficiency : 0.827 \n",
        "f(2)Thermal efficiency : 0.436 \n",
        "(g)Overall efficiency : 0.361 \n"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Ex17-pg211"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print(\"Example 4.17\")\n",
      "%pylab inline\n",
      "import warnings\n",
      "warnings.filterwarnings('ignore')\n",
      "\n",
      "import numpy\n",
      "import math\n",
      "from math import log\n",
      "import matplotlib\n",
      "from matplotlib import pyplot\n",
      "M0=0.7\n",
      "T0=228. ##in K\n",
      "p0=16. ##kPa\n",
      "eprop=0.85 ##efficiency of prop\n",
      "m=10.##Kg/s\n",
      "pd=0.98\n",
      "pc=30.\n",
      "ec=0.92\n",
      "Tt4=1600.\n",
      "Qr=42000000.## in kJ/kg\n",
      "eb=0.99\n",
      "pb=0.96\n",
      "etHPT=0.82\n",
      "emHPT=0.99\n",
      "alfa=0.79\n",
      "emLPT=0.99\n",
      "eLPT=0.88\n",
      "egb=0.995\n",
      "en=0.95\n",
      "gmc=1.4\n",
      "Cpc=1004.\n",
      "gmt=1.33\n",
      "Cpt=1152.\n",
      "z0=numpy.linspace(0.79,0.97,19)\n",
      "leng=len(z0)\n",
      "g1=numpy.zeros(leng)\n",
      "gc1=0\n",
      "g2=numpy.zeros(leng)\n",
      "gc2=0\n",
      "g3=numpy.zeros(leng)\n",
      "gc3=0\n",
      "g4=numpy.zeros(leng)\n",
      "gc4=0\n",
      "for alfa in z0:\n",
      "    Tt0=T0*(1.+((gmc-1)*(M0)**2)/2.)\n",
      "    pt0=p0*(1.+((gmc-1.)*(M0)**2)/2.)**(gmc/(gmc-1.))\n",
      "    a0=((gmc-1.)*Cpc*T0)**(1/2.);\n",
      "    V0=a0*M0\n",
      "    pt2=pt0*pd\n",
      "    Tt2=Tt0 ##Adiabatic\n",
      "    pt3=pt2*pc\n",
      "    tc=pc**((gmc-1.)/(ec*gmc))\n",
      "    Tt3=Tt2*tc\n",
      "    f=(Cpt*Tt4-Cpc*Tt3)/(Qr*eb-Cpt*Tt4)\n",
      "    pt4=pt3*pb\n",
      "    ht45=Cpt*Tt4-(Cpc*Tt3-Cpc*Tt2)/((1+f)*emHPT)\n",
      "    Tt45=ht45/Cpt\n",
      "    pt45=pt4*(Tt45/Tt4)**(gmt/((gmt-1)*etHPT))\n",
      "    m9=(1+f)*m\n",
      "    sp=(1+f)*m*eLPT*alfa*ht45*(1-(p0/pt45)**((gmt-1)/gmt))/10**6\n",
      "    Tt5=(ht45-sp*10**6/((1+f)*m))/Cpt\n",
      "    tt=Tt5/Tt45\n",
      "    et=log(tt)/(log(1-((1-tt)/eLPT)))\n",
      "    pt=tt**(gmt/(et*(gmt-1)))\n",
      "    pt5=pt45*pt\n",
      "    p9=p0 ##assumption\n",
      "    pi=p9/pt5\n",
      "    ti=pi**((gmt-1)/gmt)\n",
      "    T9i=Tt5*ti\n",
      "    T9=Tt5-en*(Tt5-T9i)\n",
      "    V9=(2*Cpt*(Tt5-T9))**(1/2)\n",
      "    Fprop=eprop*egb*emLPT*sp*10**3/V0\n",
      "    a9=((gmt-1)*Cpt*T9)**(1/2)\n",
      "    M9=V9/a9\n",
      "    pt9=p9*(1+((gmt-1)*M9**2)/2)**(gmt/(gmt-1))\n",
      "    pn=pt9/pt5\n",
      "    Fncore=m*((1+f)*V9-V0)/1000\n",
      "    spp=egb*emLPT*sp\n",
      "    Ft=Fprop+Fncore\n",
      "    Fr=Fprop/Ft\n",
      "\n",
      "    mp=((m9*V9**2)/2-m*(V0**2)/2)/10**3\n",
      "    mf=m9-m\n",
      "    PSFC=mf*10**6/((spp*10**3)+mp)\n",
      "    TSFC=mf*10**3/(Ft)\n",
      "    eth=(spp*10**3+mp)*10**3/(mf*Qr)\n",
      "    ep=(Ft*V0)/(spp*10**3+mp)\n",
      "    eo=eth*ep\n",
      "    g1[gc1]=Ft;\n",
      "    gc1=gc1+1;\n",
      "    g2[gc2]=TSFC;\n",
      "    gc2=gc2+1\n",
      "    g3[gc3]=ep\n",
      "    gc3=gc3+1\n",
      "    g4[gc4]=Fr\n",
      "    gc4=gc4+1\n",
      "\n",
      "\n",
      "pyplot.plot(z0,g1)\n",
      "pyplot.title(\"Turboprop total thrust\")\n",
      "pyplot.xlabel(\"Power split(alfa)\")\n",
      "pyplot.ylabel(\"Fprop+Fcore(kN)\")\n",
      "pyplot.show()\n",
      "pyplot.plot(z0,g2)\n",
      "\n",
      "pyplot.title(\"TSFC in turboprop engine\")\n",
      "pyplot.xlabel(\"Power split(alfa)\")\n",
      "pyplot.ylabel(\"TSFC(mg/s/N)\")\n",
      "pyplot.show()\n",
      "pyplot.plot(z0,g3)\n",
      "pyplot.plot(z0,g4)\n",
      "\n",
      "pyplot.xlabel(\"Power split(alfa)\")\n",
      "pyplot.title(\"Propeller thrust as a fraction of total thrust and propulsive efficiency\")\n",
      "pyplot.legend(\"Prop efficiency\",\"Fprop/Ftotal\")\n",
      "pyplot.show()\n",
      "##plot2d(z0,g5,4)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Example 4.17\n",
        "Populating the interactive namespace from numpy and matplotlib\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "WARNING: pylab import has clobbered these variables: ['f', 'pi', 'log']\n",
        "`%matplotlib` prevents importing * from pylab and numpy\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEZCAYAAACAZ8KHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xm8XdP9//HXWySGEHMRP0SbRlvVftNWDEGusSGG+hYR\n81B8DdUq1a9WuV/VUjVVTTXELEENTQQRwyUSbdAgVUGMicQcQ4wZPr8/1r4c17n3njucu8859/18\nPO7DHtbe+3O3k/O5e62111JEYGZmVqrF8g7AzMyqixOHmZm1iROHmZm1iROHmZm1iROHmZm1iROH\nmZm1iROHVSxJDZIOyjuOSiOpXtLV5Spv1honDutUkuZJej/7WSTpw4L1EW08XWQ/VaWtCa8dCbLZ\neyKpTtLMUst3NkkvStqyq65n+Vg87wCstkTEMo3Lkl4ADoqIe9tyDkkC1NmxtXC9xSNiQSeesq1f\n1G0t39Z7U3L5TrgX0ZbrWXXyE4d1iabVJZL6ZU8ki2XrDZJOkTQJmAeskxXtL+mfkt6VdKukFQrO\nsZOkJyXNlXSfpG8U7HtR0v9m+9+WNFLSEtm+OkmzJB0naQ5wmaReks6R9Er2c7akXk3KHy/pDUkv\nSNqzmd/z98BmwHnZU9a52fZNJD0s6R1JUyRt3Er5P0t6Ofu9H5G0aQn3uDdwB9A3O9d7klYnfZn3\nknRltu3fkr7f5F4dJ+kJ4H1JPbL/N18tKHOFpN9lyytLui27729JekDJ1cBawNjs+se2FrNVJycO\n6yql/FW9N/ATYFngJdJfrvsCBwCrAwuAxi/WAcB1wFHAysDtpC+swqfoPYFtga8BA4ATCvatCqxA\n+qI7NNs3CPhu9jOoSPmVgL7AfsDFWQxf/CUjfgNMBI6IiGUj4ihJKwLjgHOAFYGzgHGSVihWPjvV\nlCyOFbLf88bGRNaciPgAGArMzs7VJyLmZPdxJ2AUsBwwBjivyeF7ANsBy0fEwmKn5/P/h8cAM0n3\n/SvA8ZHsA7wM7JBd/4yW4rXq5cRhXaW16osAroiIpyJiUVZdEsBVEfGfiPgQ+C2we/aUMhy4LSLu\nyb7ozgCWAjYpON95EfFKRMwFfg8UtrEsAk6KiPkR8TEpyZwcEW9GxJvA/wH7NInxt1n5B0iJYPcS\nf99hwNMRcW32u40GppO+zIven6zs3Kz8WcASwLotXK/oeQpMjIg7Iw1Odw0pKX12OeDc7F59UsI1\nPiUl8n4RsTAiJpVwjNUQJw6rJE0bdZtuexnoSfpLd/VsHYDsC3EmsEYLx/YtWH8jIj4tWO9Lespp\nrvzciPioYP2lJvubKnzC6lsYazPHf+GJTNKxkv6TVW3NJT0prNzC9VrzWsHyh8CSjdWEmWL3vqnG\npPQnYAZwl6TnJP2qA3FZFXLisK4yD1i6YH21ImWKVWet1WR5PvAGMBtYu3FH1qC+JvBKC8fObuFa\ns4F+LZRfQVJh/Gs3uVahpud+pTDWIsc3TRqbAb8EdouI5SNiBeBdSmt0LnYPS6kmbFrmQ774/6ux\nrYSImBcRx0bE10hPTb+QtEUbrmVVzonDuspjwOaS1pS0HHB8kTJNvxgF7C3pm9mX9snAjdnTxY3A\nMElbSupJqnf/GJhccOzhktbI2hh+A4xuIb5RwAlZw+/KwIlA03cf/k9Sz+yLfVgWQzGvkdpVGt0O\nDJA0QtLikoYD3wBua6b8sqT2nDezRvsTgT4txN702itJKizfnl5OjwF7ZQ3lQ4HNPzuZtIOk/lmy\nfg9YSKr6K/a7WA1y4rAuERF3A9cDTwAPA2P58l+nxdavAq4A5gC9SI3hRMTTpMb0v5CeQIYBOxZ0\nJQ1So/JdwHPAs8ApLVzrFOCRLL4nsuXC8q8Cc0lPIVcDh0bEM838un8Gds16c50TEW8DO5CS25vA\nsaQG5LeLlQfuzH6eAV4EPuKLVV3Nvt8SEdNJSfD57HyrN1O+tSeDnwE7Zr/znsAtBfv6AxOA90mJ\n+vyIuD/bdyopAc+V9ItWrmFVSuWayEnSSNI/5tcjYv1s2yBSb46epL+oDo+Ih4sc+yKf/yUzPyIG\nlSVIq1lq5zskzZyrDrg6ItbscGBmNaCcTxyXk7oGFjqd1DNlIKkq4PRmjg2gLiIGOmmYmVWWsiWO\niJhIeswtNIfUOwRgeZpvXAS/fWqVxY2+ZpmyVVVBejsYGFtQVbU28CDpH+FiwMYR8aVugJKeJ/Ui\nWQj8NSIuKVuQZmbWJl3dOH4ZcFRErAUcDYxsptzgrDprO+CIrBeLmZlVgK5+4ngvIvpkywLeiYjl\nmj8DSDoJmBcRZxbZ5+oDM7N2iIh2Nwd09RPHDElDsuUtSd0Nv0DS0pKWzZZ7k8YamtbcCSPCPxGc\ndNJJucdQKT++F74Xvhct/3RU2YZVlzQKGAKsrDQ/wInAIcD5SqOUfpStI6kvcElEDCO9UXxzeiBh\nceDaiLirXHGamVnblC1xRERzk/ZsWKTsbNI7H0TE88B/lSsuMzPrGL85XiPq6uryDqFi+F58zvfi\nc74XnaesjePlJimqOX4zszxIIqqocdzMzKqcE4eZWSd66628Iyg/Jw4zs04QASeeCDvumJZrWdl6\nVZmZdReffgoHHwzTp8PYsaAaH2nPicPMrAPefRd23RWWXhruvRd69847ovJzVZWZWTvNmgWbbw4D\nBsDNN3ePpAFOHGZm7TJtGmyyCey9N5x3HvTokXdEXcdVVWZmbXTPPTBiBJx7LuyxR97RdD0/cZiZ\ntcHVV8Oee8KNN3bPpAF+4jAzK0kE/OEPcMklcN998K1v5R1Rfpw4zMxasWABHHYYPPooPPQQrL56\n3hHly4nDzKwF8+bB7runJ44HHoBllsk7ovy5jcPMrBmvvgpDhsAaa8CYMU4ajZw4zMyKeOop2Hhj\n2GUXuPhi6Nkz74gqh6uqzMyaeOAB2G03OP102G+/vKOpPE4cZmYFRo+Go46C666DrbfOO5rK5MRh\nZkZq/D7jjPRS3913w3e+k3dElatsbRySRkp6TdK0gm2DJE2RNFXSw5I2aObYoZKmS3pW0q/KFaOZ\nGcDChfDTn6aX+x56yEmjNWWbOlbSZsA84KqIWD/b1gCcGhHjJW0HHBcRWzQ5rgfwNLA18ArwMDAi\nIp4qcg1PHWtmHfLhh+lN8Hnz4KabYLnl8o6o/Cp26tiImAjMbbJ5DtD4v2V5UmJoahAwIyJejIj5\nwGhg53LFaWbd1xtvwJZbQp8+cPvt3SNpdIau7o77v8CZkl4G/gQcX6TMGsDMgvVZ2TYzs07z7LOp\nu+0228CVV0KvXnlHVD26unH8MuCoiLhF0m7ASGCbJmXaVPdUX1//2XJdXR11dXUdDNHMat1DD6X3\nM373uzRzX61raGigoaGh085XtjYOAEn9gLEFbRzvRUSfbFnAOxGxXJNjNgLqI2Jotn48sCgi/ljk\n/G7jMLM2ueUWOOSQ9JSx/fZ5R5OPim3jaMYMSUOy5S2BZ4qUeQT4uqR+knoBw4ExXRWgmdWuc8+F\nI4+EO+/svkmjM5StqkrSKGAIsLKkmcCJwCHA+ZKWAD7K1pHUF7gkIoZFxAJJRwLjgR7AZcV6VJmZ\nlWrRIvjlL1MD+KRJ0K9f3hFVt7JWVZWbq6rMrDUffwz77psGLLz1Vlhxxbwjyl+1VVWZmXWZt95K\nw4ZIcNddThqdxYnDzGrSCy/A4MGwySYwahQsuWTeEdUOJw4zqzmPPJKSxk9/mka4XczfdJ3Kgxya\nWU257TY44AC49FLY2WNOlIXzsJnVjL/+Nb3Qd9ttThrl5CcOM6t6ixbBCSfAjTfCxInQv3/eEdU2\nJw4zq2qffgoHHgjPPQeTJ8Mqq+QdUe1zVZWZVa133oGhQ9PQ6Pfc46TRVZw4zKwqvfwybLoprL9+\nqqJaeum8I+o+nDjMrOo89ljqbnvggXDOOdCjR94RdS9u4zCzqnLXXbD33nD++bDbbnlH0z35icPM\nqsYVV6Rxp26+2UkjT37iMLOKFwEnn5zm0GhogG98I++IujcnDjOraPPnw6GHwrRpaea+VVfNOyJz\n4jCzivXee6lKqmfP9KTRu3feERm4jcPMKtTs2bD55rDOOmkeDSeNyuHEYWYV58knYeONYfhwuPBC\nWNx1IxXF/zvMrKLcd19KGGefDXvtlXc0VoyfOMysYlx3XUoao0c7aVSysj1xSBoJDANej4j1s22j\ngXWzIssD70TEwCLHvgi8BywE5kfEoHLFaWb5i4DTToOLLoJ774VvfzvviKwl5ayquhz4C3BV44aI\n2KNxWdIZwDvNHBtAXUS8Xcb4zKwCLFiQZuqbPDn9rLFG3hFZa8qWOCJioqR+xfZJErA7sEULp1AZ\nwjKzCvLBB7DHHvDJJ2kejT598o7ISpFXG8dmwGsR8Vwz+wO4W9Ijkg7uwrjMrIu89hrU1cHKK8O4\ncU4a1SSvXlUjgOta2D84IuZIWgWYIGl6REwsVrC+vv6z5bq6Ourq6jozTjMrg6efhu22g332gfp6\nkOsXyqqhoYGGhoZOO58iotNO9qWTp6qqsY2N49m2xYFZwPciYnYJ5zgJmBcRZxbZF+WM38w636RJ\n8OMfwx/+kIZFt64niYhod7rOo6pqa+Cp5pKGpKUlLZst9wa2BaZ1YXxmViZ/+xvssksarNBJo3qV\nLXFIGgVMBgZIminpgGzXcGBUk7J9JY3LVlcDJkp6DPgncFtE3FWuOM2sa5x9Nvz85zB+PPzwh3lH\nYx1R1qqqcnNVlVnlW7gQjjkGJkyAO+6AtdbKOyLraFWVhxwxs7L56KM0W99bb8GDD8IKK+QdkXUG\nDzliZmXx5puw1VawxBKpespJo3Y4cZhZp3vuOdhkExgyBK65JiUPqx1OHGbWqaZMgU03haOPhlNP\nhcX8LVNz3MZhZp1mzBg46CAYORJ23DHvaKxcnDjMrFNccAGccgrcfjtssEHe0Vg5OXGYWYcsWgTH\nH5+md33wQfjqV/OOyMrNicPM2u2TT2D//eHll9OQ6CutlHdE1hXcbGVm7TJ3bnoDfP58uPtuJ43u\nxInDzNrspZdg8GD43vfghhtgqaXyjsi6khOHmbXJv/6Vksahh8JZZ7m7bXfUYhuHpO+R5s7YHOhH\nmmDpJeAB4LqImFruAM2sctx5Z5pD46KL0tDo1j01O8ihpNuBucAYYAowhzSd6+rAIGBHYPmIGNY1\noRaN0YMcmnWRyy6D3/wGbr45vRVu1aujgxy2lDhWjYjXWrn4VyLi9fZevKOcOMzKLwJOOgmuvTaN\nbjtgQN4RWUeVbXTc1pJGVia3pGFm5ffpp3DwwfDUU/DQQ/CVr+QdkVWCZhOHpBdaOC4iwq/5mNWw\nd9+FXXdNPabuuw969847IqsULTWOFw4aEKQeWMOBY4F/lTMoM8vXrFmw/faw2WZw7rnQo0feEVkl\nabYjXUS8GRFvAm+TGsIbgI2B7SPC/SnMatQTT6TG7733hvPOc9KwL2upqqoXcCBwNPAgsHNEzOiq\nwMys691zD4wYAX/+c/qvWTEt9aqaBSwA/gy8TKqugtQlNyLi5hZPLI0EhgGvR8T62bbRwLpZkeWB\ndyJiYJFjhwLnAD2ASyPij81cw72qzDrJVVfBL3+Z3gQfMiTvaKycytkd94pssWiBiDiglcA2A+YB\nVzUmjib7zyAljlOabO8BPA1sDbwCPAyMiIinipzDicOsgyLg97+HSy9NQ6J/61t5R2TlVs7uuPtn\nF1gyIj5uctFWhzOLiImS+hXbJ0nA7sAWRXYPAmZExItZ2dHAzsCXEoeZdcyCBXDYYfDoo2l02759\n847IqkEpo8zcLKln44qk1YEJHbzuZsBrEfFckX1rADML1mdl28ysE82bBzvtlHpQ3X+/k4aVrpT5\nOG4BbpC0K7AmaQiSYzt43RHAdc3sa1PdU319/WfLdXV11NXVtTsos+7i1Vdh2DAYOBAuvBB69mz9\nGKteDQ0NNDQ0dNr5mm3j+EIh6UhgKLA28D8RMamkk6eqqrGFbRySFic9RXwvImYXOWYjoD4ihmbr\nxwOLijWQu43DrO2eeiq9o3HggXDCCaB213RbtSpbG4ekY7LFIPWkWhN4HNhI0oYRcVY7r7k18FSx\npJF5BPh6lnRmk146dMdAs07wwAOw225w+umw3355R2PVqqU2jmWBZQr+ewvwbMG2FkkaBUwGBkia\nKamxF9ZwYFSTsn0ljQOIiAXAkcB44D/A9cV6VJlZ21x/fRpC5JprnDSsY0qqqqpUrqoya10EnHlm\nGjrkttvgO9/JOyLLW0erqpp94pA0UtIGLezfUNLl7b2wmZXfwoVw1FFw5ZWpu62ThnWGlnpVnQ38\nMmusfprPJ3JajfT292TgjLJHaGbt8uGHsOeeqdvtgw/CcsvlHZHVilarqiQtAQwk9ahqnDr28aYv\nBebBVVVmxb3xBuy4Y5p06dJLoVevvCOySlK2IUeaXGRpYK2ImN7eC5WDE4fZlz37LGy3HeyxB/zu\nd+5ua19WtjaOggvsBEwF7szWB0oa094Lmln5/OMfsPnmcNxxcMopThpWHqUMOVIPbAjMBYiIqYBn\n/zOrMLfckqqnLrsMDjkk72islpUy5Mj8iHhHX/zTZVGZ4jGzdvjLX+C00+DOO+H73887Gqt1pSSO\nJyXtBSwu6evAUaQeVWaWs0WLUrXUuHGp59Q66+QdkXUHpVRVHQmsB3xCeuP7PeDn5QzKzFr38cep\nAXzKFJg0yUnDuk6LvaqyAQknRESxeTNy515V1l29/TbsvHMaCv3KK2HJJfOOyKpJWXtVZeNGLZK0\nfHsvYGad64UXYJNNYKONYNQoJw3reqW0cXwATJM0IVuGNOf4UeULy8yKeeSRNPnSr38NRx6ZdzTW\nXZWSOG7OfhrrhEQbJ1sys44bNw723x8uuQR+9KO8o7HurNQ3x5cABmSr0yNiflmjKpHbOKy7uPhi\nOOmk9K7GRhvlHY1Vu7JN5FRwgTrgStIYVQBrSdovIu5v70XNrDQRaZa+66+HiROhf/+8IzIrrarq\nLGDbiHgaQNIAYDTwvXIGZtbdffppmt51xgx46CFYZZW8IzJLSnmPY/HGpAEQEc9QWsIxs3Z65x0Y\nOjQNiX7vvU4aVllKSRyPSrpUUp2kLSRdSpoX3MzKYOZM2HRTWG89uOkmWHrpvCMy+6JS5uNYEjgC\nGJxtmghcEBGflDm2Vrlx3GrN44/DDjvAz38Ov/iFR7e18ij7fBySegMfR8TCbL0HsEREfNjKcSOB\nYcDrEbF+wfafAocDC4FxEfGrIse+SBraZCFpkMVBzVzDicNqxoQJsNdecN55sPvueUdjtazs83EA\n9wJLFawvDdxdwnGXA0MLN0jaAtgJ+E5EfJvmp54NoC4iBjaXNMxqyRVXwD77pKopJw2rdKU0ci8R\nEfMaVyLi/WxGwBZFxERJ/ZpsPgw4tfE9kIh4o4VT+CHdal4EnHxyGm+qoQG+8Y28IzJrXSlPHB9I\n+myEf0k/AD5q5/W+Dmwu6R+SGrJzFRPA3ZIekXRwO69lVtHmz4ef/ATGjoXJk500rHqU8sTxc+AG\nSXOy9dWB4R243goRsZGkDYAbKD6b4OCImCNpFWCCpOkRMbHYCevr6z9brquro66urp2hmXWd99+H\n3XaDHj3Sk8Yyy+QdkdWyhoYGGhoaOu18pQ450gtYN1t9OiI+LenkqapqbGPjuKQ7gNMa3zqXNAPY\nMCLeauEcJwHzIuLMIvvcOG5VZ/ZsGDYMBg2C88+Hxf1WlHWxsjWOS/pDweqQiJiW/ZSUNJpxK7Bl\ndv4BQK+mSUPS0pKWzZZ7A9sC0zpwTbOK8eSTaUj03XeHiy5y0rDq1FIbx3YFy6e39cSSRpGmmB0g\naaakA4CRwFclTSPNJrhvVravpHHZoasBEyU9BvwTuC0i7mrr9c0qTUMDbLklnHIKHH+839Gw6tVs\nVZWkqRExsOlyJXFVlVWL665LL/WNHp2Sh1meyjk67iqSfkHqFlu4DGkip7Pae1Gz7iIC/vhHuOCC\nNObUt7+dd0RmHddS4rgUWLbIspmVYMECOOoomDQpjW67xhp5R2TWOUrqVVWpXFVlleqDD2CPPeCT\nT+Bvf4M+ffKOyOxzXTHkSOOF9mvvRcy6k9degy22gJVWStO9OmlYrSk5cZBeBDSzFjz9NGy8MWy3\nHVx+OfTsmXdEZp3PvcjNOsmkSfDjH8Pvfw8HHZR3NGbl02LikHRfwWr/gvWICHcqNMvcdBMcdhhc\nfTX88Id5R2NWXi02jmdDhgSpG+44YPtsmYh4sezRtcKN41YJzj4bzjwzDVY4sOLedjL7snK+x/GF\n5CDp04h4qb0XMqs1CxfCMcfAXXelaqq11847IrOu0ZY2jhfLFYRZtfnoozTx0ptvpqSxwgp5R2TW\ndUruVRURu5QzELNq8eabsPXWqcfU+PFOGtb9tJo4lPxY0tmSzpK0i+Th2ax7eu65NLrtZpvBtdfC\nEkvkHZFZ12v1zXFJFwJfI41mK2B34PmIOLz84bXMjePWlaZMgR/9CE44AQ7P/dNv1n4dbRwvJXFM\nB74VEYuy9cWA/0RE7hNdOnFYVxkzJr2bcdllsNNOeUdj1jFdMeTIDGCtgvW1sm1m3cIFF8Chh6bh\nQ5w0zErrVdUHeErSFNI7HYOAhyWNJb0I6H9KVpMWLYJf/xpuvhkefBC+9rW8IzKrDKUkjhOz/zbW\nCRU+3rieyGrSJ5/AAQfASy/B5Mmw8sp5R2RWOUoaVl3SasAGpEQxJSJeL3dgpXAbh5XD3Lmwyy4p\nWVx9NSy1VN4RmXWusrdxSNqdNPf3bqQeVVMk7dbeC5pVspdegk03TUOHXH+9k4ZZMaU0jp8AbBAR\n+0bEvqQnj9+2dpCkkZJekzStyfafSnpK0r8l/bGZY4dKmi7pWUm/KuUXMeuoqVNh8GA4+OA0/lSP\nHnlHZFaZSkkcAt4oWH+LL7ZzNOdyYOgXTiRtAewEfCcivg2c8aWLST2A87JjvwWMkPTNEq5n1m7j\nx6dRbf/8Z/i5Z54xa1EpjeN3AuMlXUdKGMOBO1o7KCImZqPrFjoMODUi5mdl3mh6HKnX1ozGARYl\njQZ2Bp4qIVazNhs5MvWeuvXW9Fa4mbWstfk4BPyFVD01ONv814i4pZ3X+zqwuaQ/AB8Dx0bEI03K\nrAHMLFifBWzYzuuZNSsC6uvhmmvg/vth3XXzjsisOpTyxHF7Vq10Uyddb4WI2EjSBsANwFeblGlT\nN6n6+vrPluvq6qirq+tgiNYdzJ8PhxwCTz6ZutuuumreEZmVT0NDAw0NDZ12vlKGHLkSOD8iprT5\n5KmqamxErJ+t3wGcFhH3Z+szgA0j4q2CYzYC6iNiaLZ+PLAoIr7UkO7uuNYe770Hu+4KSy4Jo0ZB\n7955R2TWtbpiyJGNgIckPS9pWvbzRDuvdyuwJYCkAUCvwqSReQT4uqR+knqR2lTGtPN6Zl/wyitp\nZNv+/dMb4U4aZm1XSlVV4wzKxd4cb5akUcAQYCVJM0lvoI8ERmZddD8F9s3K9gUuiYhhEbFA0pHA\neKAHcFlEuGHcOmzaNNhhhzSy7XHHgScHMGufUt8c/z6wKbAImBQR/yp3YKVwVZWV6t57YY89Unfb\nESPyjsYsX13x5viJwBXAisAqwOWSWn0B0KxSXHNNShY33uikYdYZSmkcf4b0wt7H2fpSwOMRMaAL\n4muRnzisJRFw6qlw8cVpSPT11ss7IrPK0NEnjlLaOF4BliK9dwGwJOndCrOKtWABHHFEmrVv8mTo\n2zfviMxqRymJ4z3gSUl3ZevbkAY6/AtpPo6jyhadWTvMmwfDh8PChfDAA7DssnlHZFZbSqmq2r9g\nNUi9qj77b0RcWbboWuGqKmvq1VdTz6nvfhcuugh69sw7IrPKU7Y5xyWtFREvtzuyLuDEYYWmT4ft\nt4f994ff/tbdbc2aU85eVX8vuEhnDDdiVjYTJ8KQIXDSSXDiiU4aZuVUShsHfHk8KbOKccMNcOSR\ncO21sM02eUdjVvtKTRxmFScCzjoLzjkHJkxI7RpmVn4ttXEsBD7MVpcCPirYHRHRp8yxtcptHN3X\nwoVw9NFw331w++2w5pp5R2RWPcr2HkdEeOJMq0gffgh77ZVGuZ04EZZfPu+IzLqXUkbHNasYb7wB\nW20FyywDd9zhpGGWBycOqxozZqSpXbfaCq66Cnr1yjsis+7JicOqwj/+kebR+OUv4ZRT3N3WLE/u\nVWUV79Zb4eCD4cor0wt+ZpYvJw6raOedl0a4veMO+MEP8o7GzMCJwyrUokXwq1/B2LHw4IOwzjp5\nR2RmjZw4rOJ8/DHstx/Mnp2GRF9xxbwjMrNCbhy3ivL227Dttml5wgQnDbNKVLbEIWmkpNckTSvY\nVi9plqSp2c/QZo59UdITWZkp5YrRKsuLL8LgwbDhhjBqFCy5ZN4RmVkx5XziuBxomhgCOCsiBmY/\ndzZzbAB1WZlBZYzRKsSjj6akccQR8Kc/wWJ+FjarWGVr44iIiZL6FdlVag9899TvJm6/Pc2hcfHF\n8KMf5R2NmbUmj7/rfirpcUmXSWpuwIgA7pb0iKSDuzI461qXXAIHHQRjxjhpmFWLru5VdSFwcrb8\nO+BM4KAi5QZHxBxJqwATJE2PiInFTlhfX//Zcl1dHXV1dZ0asJVHRJql7/rr00CF/fvnHZFZ7Wpo\naKChoaHTztfqnOMdOnmqqhobEeu3ZV+TcicB8yLizCL7PKx6Ffr0U/jJT+CZZ9J7GquskndEZt1L\nOaeO7XSSVi9Y3QWYVqTM0pKWzZZ7A9sWK2fV6d1307Ah770H997rpGFWjcrZHXcUMBlYV9JMSQcC\nf8y62T4ODAGOzsr2lTQuO3Q1YKKkx4B/ArdFxF3litO6zqxZaaDCb34TbroJll4674jMrD3KWlVV\nbq6qqh5PPAE77AA/+xn84hce3dYsT2WbAdCss9x9N+y5ZxqwcPfd847GzDrKr1lZWV15ZZrm9aab\nnDTMaoWfOKwsItKESyNHQkNDatcws9rgxGGdbv58OPxwmDoVHnoIVlst74jMrDM5cVinev/9VCW1\n2GLpSWP+RsEmAAANEklEQVSZZfKOyMw6m9s4rNPMmQNDhsCaa8Lf/+6kYVarnDisU/znP7DxxrDr\nrvDXv8LifpY1q1n+520ddv/9qXrqzDNh773zjsbMys2Jwzpk1Kj0Ut/o0bDllnlHY2ZdwYnD2iUC\nTj8dzj8f7rkH1m9xqEozqyVOHNZmCxbAUUfBpEmpu+0aa+QdkZl1JScOa5MPPoARI+Djj9M8Gn36\n5B2RmXU196qykr3+emrHWHFFGDfOScOsu3LisJI880zqbvvDH8Lll0PPnnlHZGZ5cVWVtWryZPjv\n/4bf/z7ND25m3ZsTh7XoppvgsMPgqqtg6NC8ozGzSuDEYc065xw44wwYPx4GDsw7GjOrFE4c9iWL\nFsExx6SEMWkSrL123hGZWSUp55zjIyW9JmlawbZ6SbMkTc1+ilZ+SBoqabqkZyX9qlwx2pd99FEa\nPmTqVCcNMyuunL2qLgeaJoYAzoqIgdnPnU0PktQDOC879lvACEmeBqgLvPUWbL116jE1fjyssELe\nEZlZJSpb4oiIicDcIrtamyB9EDAjIl6MiPnAaGDnzo7Pvuj552GTTWCzzeDaa2GJJfKOyMwqVR7v\ncfxU0uOSLpO0fJH9awAzC9ZnZdusTB5+GDbdNA1WeNppaRImM7PmdPVXxIXAOsB/AXOAM4uUiS6N\nqJu77TYYNgwuuihN92pm1pou7VUVEa83Lku6FBhbpNgrwJoF62uSnjqKqq+v/2y5rq6Ourq6jobZ\nbVx0EZx8ckoegwblHY2ZlUtDQwMNDQ2ddj5FlO8PfEn9gLERsX62vnpEzMmWjwY2iIg9mxyzOPA0\nsBUwG5gCjIiIp4qcP8oZf61atAh+85v0ct8dd8DXvpZ3RGbWlSQREa21NzerbE8ckkYBQ4CVJc0E\nTgLqJP0XqTrqBeDQrGxf4JKIGBYRCyQdCYwHegCXFUsa1j6ffAIHHggvvJCGEll55bwjMrNqU9Yn\njnLzE0fbvPMO7LJLGt32mmtgqaXyjsjM8tDRJw73n+kmXn4ZBg+G734XbrjBScPM2s+Joxt47LH0\njsbBB6fxp3r0yDsiM6tmHquqxo0fD/vsAxdeCD/+cd7RmFkt8BNHDRs5EvbbD265xUnDzDqPnzhq\nUATU16cG8Pvvh3XXzTsiM6slThw1Zv58OOQQePLJ1N121VXzjsjMao0TRw157z3YdVdYckm47z7o\n3TvviMysFrmNo0a88koa2bZ/f7j5ZicNMysfJ44aMG1a6m67555w/vmwuJ8jzayM/BVT5e69F/bY\nI72fseeerZc3M+soJ44qds01aW7wG2+EIUPyjsbMugsnjioUAaeeChdfnJ441lsv74jMrDtx4qgy\nCxbAEUfAlCmpu23fvnlHZGbdjRNHFZk3D4YPh4UL4YEHYNll847IzLoj96qqEq++CnV1sNpqMHas\nk4aZ5ceJowpMn5662+60E1x6KfTsmXdEZtaduaqqwk2cmN4GP/30NGChmVnenDgq2I03pobwa6+F\nbbbJOxozs8SJowJFwFlnpZf6JkxIs/aZmVWKsrVxSBop6TVJ04rsO0bSIkkrNnPsi5KekDRV0pRy\nxViJFi6En/0Mrrgidbd10jCzSlPOxvHLgaFNN0paE9gGeKmFYwOoi4iBETGoTPFVnA8/hN12g3//\nO7VtrLlm6cc2NDSULa5q43vxOd+Lz/ledJ6yJY6ImAjMLbLrLOC4Ek6hzo2o8r30EnzlK3DnnbD8\n8m071v8oPud78Tnfi8/5XnSeLu2OK2lnYFZEPNFK0QDulvSIpIO7ILSK8M1vwkUXQa9eeUdiZta8\nLmscl7Q08GtSNdVnm5spPjgi5khaBZggaXr2BGNmZjlTRJTv5FI/YGxErC9pfeBu4MNs9/8DXgEG\nRcTrLZzjJGBeRJxZZF/5gjczq2ER0e7mgC574oiIacBnM2BLegH4fkS8XVguezLpERHvS+oNbAv8\nXzPn7HbtIGZmeStnd9xRwGRggKSZkg5oUiQKyvaVNC5bXQ2YKOkx4J/AbRFxV7niNDOztilrVZWZ\nmdWeihzkUNJQSdMlPSvpV0X2H5u9HDhV0jRJCyQtn+2rqZcHS7gXK0u6U9Jjkv4taf9Sj602HbwX\n3e1zsYKkWyQ9LumfktYr9dhq08F7UTOfi5Zeui4oc252nx6XNLBge9s+ExFRUT9AD2AG0A/oCTwG\nfLOF8jsAdxesvwCsmPfv0VX3AqgHTs2WVwbeIrVdtek+VvpPR+5FN/1c/An4bba8buO/kW76uSh6\nL2rwc7EZMBCY1sz+7YHbs+UNgX+09zNRiU8cg4AZEfFiRMwHRgM7t1B+T2BUk2210mheyr2YA/TJ\nlvsAb0XEghKPrSYduReNutPn4pvAfQAR8TTQT9JXSjy2mrT3XqxSsL8mPhfR/EvXjXYCrszK/hNY\nXtJqtOMzUYmJYw1gZsH6rGzbl2Q9sH4I3FSwuZZeHizlXlwCrCdpNvA48LM2HFtNOnIvoPt9Lh4H\n/htA0iBgbVIX+O74uWjuXkBtfS5a09y96tvM9mZV4ui4bWmt3xF4MCLeKdhWSy8PlnIvfg08FhF1\nkr5G+p1rcWjEdt+LiHif7ve5OA34s6SpwDRgKrCwxGOrSUfuBcCmETG7Rj4XpeiUp6tKfOJ4BSgc\n3m9NUgYsZg+aVFNFxJzsv28At5Aew6pVKfdiE+BGgIh4jlRnu25WrtT7WA06ci+63eciIt6PiAMj\nDRS6L7AK8Fwpx1aZ9t6L57N9s7P/1sLnojVN79X/I92rNn8mKjFxPAJ8XVI/Sb2A4cCYpoUkLQds\nDvy9YNvSkpbNlhtfHmy2h0EVKOVeTAe2BpC0KumL8vkSj60m7b4X3fFzIWm5bB9ZFcz9ETGvlGOr\nTLvvRQ1+LlozBtgXQNJGwDsR8Rrt+ExUXFVVRCyQdCQwntTaf1lEPCXp0Gz/X7OiPwLGR8RHBYev\nCtwiCdLvdm1U8cuDJd6LPwCXS3qc9IfAcZG9jV/s2Dx+j87QkXsh6avAzd3sc/Et4AqlYXn+DRzU\n0rF5/B6doSP3ghr7vlB66XoIsLKkmcBJpF5SRMRfI+J2SdtLmgF8AByQ7WvzZ8IvAJqZWZtUYlWV\nmZlVMCcOMzNrEycOMzNrEycOMzNrEycOMzNrEycOMzNrEycOq1qSFurzofVvkLRU3jG1RNL+kv6S\nLR8qaZ+C7as3KXt9NmxKS+drkPT9bHk3Sf+RdE8L5VeVdHvHfxPr7pw4rJp9mA0jsT7wKfA/5byY\npI6+MPvZS1PZC1lXZ6v7kQaaa7xOf6B3NmxKa+drPOdBwE8iYqtmC6e3hOdK+l57gjdr5MRhteJB\noL/SpD23ZhPVPCRpfYBssp4+St4q+Gv/KklbSVpM0p8kTcmOPSTbXydpoqS/A08WXlBSD0lXZE88\nT0j6Wba9QdI5BU9DGzQeUnBsvaRjJP0Y+AFwraR/SVqSNAbbmIKyF0h6WGlyqvomv7cknQgMBkZK\nOl3S2pIekPRo9rNxQfkxwIgO3mvr5pw4rOplTwJDgSeAk4FHI+K7pNFyr8qKTQI2BdYjDfa3abZ9\nI2Ay8BPS2D2DSAPdHSypX1ZmIHBURKzb5NL/BfSNiPUj4jvA5dn2AJaKiIHA4cDIImEHEBFxE2ms\noD0j4nsR8TEpCTxSUPY3EbEB8F1gSGMybDxPRJxccI7jgNeBbSLi+6QkdG5B+SmkMd7M2q3ixqoy\na4OllIbKBniA9AX9T7K5FyLiPkkrZQPZTSR9Yb4EXAgcIqkvMDciPpK0LbC+pF2z8/UB+gMLgCkR\n8VKR6z8HfFXSucA4oHCco1FZDBOzJ53lWvldCoe7Xps0KVWj4dngfIsDq5MmJio2GF/jOXoB5ykN\nr78QGFBQZg5ppjezdnPisGr2UfZX/WeyAeuazjkQpMRyJPAi8BtgF2DXbHujIyNiQpPz1ZEGhPuS\niHhH0ndITzv/A+zO5wPoNbWold+l6aBxyq6/DnAM8IOIeFfS5cCSrZzjaGBOROwjqQfwcZPzeoA6\n6xBXVVmtmQjsBZ996b8REfMiYhZpHvL+EfECqU3kWD5PHOOBwxsbwCUNUJphslmSViLNaX4z8FtS\nlRakL+fhWZlNSVVg7zc9nM8T3Pt8PuUtpKeixl5WfUiJ6z2loeK3K+Ee9AFezZb3JY142mj17Pxm\n7eYnDqtmxf5yric1Ej9O+sLdr2DfP/j8j6UHScOwP5itX0qqwvmX0mPL66SnksKeS02tQRrGvfGc\n/1sQ18eS/kX6N3ZgwfYosnwFcJGkD0mTUT1IajB/NCIez6rjppOm92yMtyUXADdJ2he4E5hXsG8Q\nX3zKMmszD6tu1skk3QccExH/aufxXwX+EhHDOjcykHQtcEZETG21sFkzXFVlVmEi4nng/dZeAGwr\nSV8BlnfSsI7yE4eZmbWJnzjMzKxNnDjMzKxNnDjMzKxNnDjMzKxNnDjMzKxNnDjMzKxN/j9Wa4hN\nQwruUQAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x2d8afb0>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEZCAYAAACAZ8KHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xm8ndPZ//HPN5M5MSeCSg3xM0sRsxxzaqwxDapmrdIJ\nTWnVeQyPtmhVaZ8iiaGEeoyRlAgOQdEQQxWPqJAQaUrMVIbr98e6j+xu+8xn73vvc77v12u/cu97\nvPZt29dZa91rLUUEZmZmrdUj7wDMzKy2OHGYmVmbOHGYmVmbOHGYmVmbOHGYmVmbOHGYmVmbOHFY\nzZH0gaRBOV6/QdKxeV2/2kj6UvbfRHnHYpXhxGGtIunD7MfhA0mLJH1c8H6kpOUljZE0W9L7kl6S\nNKrg+EVF53inYFtfSZdIei3bNl3SryWtVCqWiFguIma04zMMyuLo6Pc+spcBEfF69t/E96SbcOKw\nVomIZbMfh+WA14B9Gt9HxDjg18DSwP+LiL7AfsD0otNsWnDMigCS+gD3ARsAe2bn3xb4FzC0TB+n\nXX8ZK6nY/zOSelXqWmZt4cRhnWVLYFxEvAcQES9FxC2tOO5IYE3ggIh4MTt2bkScHxF/LnVAVmpY\nO1u+WtLlku7KSjqPNW4r4aHs33ezfbeRVC/puoJz/0epJKuWOk/SI8CHwJezXdeV9Lik9yTdLmmF\ngnPsJ+l5SfMkPSDp/xVsmyHpx9n2d7JS2hLZtjpJsyT9SNJsYLSkPllp7I3s9ess2Rbuf4akuZJe\nlXRYUzdaUj9JoyW9mR13bsHnPErSw5IuzOL6h6ThBcd+WdJD2X27N7vn1zVzz87Jzve+pHsKS4/Z\nfX80uz9PSxrWVMxWnZw4rLM8Bpyf/QCt18Q+pf7S3w34c0R83IFrjwDqgRVIpZzzm9hvx+zffhHR\nNyIeo3VVTkcAxwGNpS2REt7RwGrAAuBSAEmDgRuA7wIrAxOB8UWlh8OAPYB1gMHATwu29c8+x5eA\nE7NtQ4HNstfQEvuvBAwEvglckcVQytXAZ9l1h2QxHFewfSjwYna+XwKjC7bdQPpvvCLpXh9B8/du\nJHAUsCrQBzgNQNLqwF3AORGxQrb+FkkrN3MuqzJOHNZZTgGuB04Gnpf0cuFfrJmnsr8y50m6JFu3\nIjC7A9cN4NaImBoRC7MYNm9i31KJq6VqqwCujogXImJRRCzI1l0bEX/PEt5ZwKHZX9wjgLsi4r4s\nnouApYDtCs53WUS8ERHzSEluZMH1FgFnR8T8iPiUlGTOiYh/RcS/gP8CvlEU41nZ/g8BE4BDv/Ah\npf7AV4EfRMQnETEXuAT4esFur0XE6Kyt4lpgNUmrSvoSqUT5s4hYEBGPAHc2c+8CGBsR07PP8CcW\n/zc5ApgYEXcDRMRkYCqwVxPnsirkOlTrFNkPxAXABZKWA34M3CxpzYh4N9ttSET8o+jQt0l/LXfE\nnILlT4BlO3i+YjNbWPc60JtUwlgtew9ARISkmcDqzRxb+PnnRsRnBe8Hkko5Te0/LyI+KXj/GqXv\n51pZjLO1+OGnHoWxAm8VxP1xtt+ypFLDO9l/48LPsGaJ63zhXPznf5O1gEMk7VuwvRdwfzPnsirj\nEod1uoj4gJRElmFxm0BTJgN7Slq67IGVrlr5kNSo32hAK4/7UtHyfGAu8CbpxxFIDeqkH9g3mjn2\nzWau9SYwqJn9Vyi6d2sVXavRTODfwEoRsUL26hcRm5TYt9hsYEVJSzXxGdrideC6ghhWyB6W+GU7\nz2c5cOKwTiHpLElbZo25SwLfA+YBL7Vw6HWkH7VbJK0vqYeklSSdKemrrbl0G8KcS6oKWqdg3dPA\nTpLWlNQPOKMV1xBwhKQNsh/tc4Cbsyqem4G9Je0iqTdwKvAp8GjBsSdJWl3SisBPgBubiXkc8FNJ\nK2ftAD8j3bNC/yWpt6Qdgb2zGP5DRMwGJgG/krRcdp/XkbRTM9duPPY1UnVSfXadbYF9aL6No6n/\nLn8E9pW0h6SekpbMGvlXb2J/q0JOHNZZFgFjST/ObwC7AnsXNHqX/JHJqmV2IzXK3gu8BzxOavt4\nrIlrRdFy8bmbutbHpDaFR7J2lqFZHftNwLPAX4HxrThfYxvA1aS/xvuQGsOJiJdI9fi/Jd2LvYF9\ns7aRxmNvIP2IvwK8DJzXzLXOI/1oP5u9phbt/xYpQb9JSignRsT/lfr8pAb9PsDfgXdICaaxhNXS\nfTyc9Jj028C5pHv2WRP7Fr///NwRMQvYHzgT+CepBHIq/i2qKSpXnx1JY0j/0/yzsTgsaTPgf0hV\nGDOAw7NqjeJjZwDvAwuB+RFRruf5zSpK0qvAsRHR4Tp9SXWkap/m2hrKQtJNwN8j4r8qfW3LXzmz\n/Fig+Kmaq4AfRcSmwG3A6U0cG0BdRAxx0jDLX1YNuU5WxfVVUgfP2/OOy/JRtsQREVNIRehC62Xr\nITWKHtTMKTzujVnLKjXMxwDgAeAD0igB34qIZyp0basylX4c93lJ+0fEHcAhNP04XwCTJS0E/hAR\nV1YsQrMyioiWnjJry7kaaP/TTW291l2kjntmFW+QOob0RMlU0nPdnzWx3/YRMYTUYek72dMiZmZW\nBSpa4sieONkTPh+aYe8m9pud/TtX0m2koRCmFO8nyaNxmpm1Q0S0uzmgoiUOSatk//Ygjbfz+xL7\nLJ31PEbSMqTxdJ5r6pwR4VcEZ599du4xVMvL98L3wvei+VdHlS1xSBpH6vS0vqSZko4BRkp6CXgB\nmBURV2f7DpQ0ITt0ADBF0tOk5/nviohJ5YrTzMzapmxVVRExsolNl5bY902yaqtIYxk1NUidmZnl\nzL01q0gEvPVWy/uVUldX16mx1DLfi8V8Lxbzveg8Zes5XgmSopbjLzZtGuyzD9x/P6y/ft7RmFlX\nJYmolcZxa96QIXD++bDbbjC9eNJVM7Mq4fk4qsxRR8Fnn8Guu8KDD8KgQXlHZGb2n5w4qtAJJ8D8\n+bDLLil5rFnxIezMzJrmxFGlvvOdVPLYZRdoaIDVPVuBmVUJJ44q9oMfLK62amiAAaXmpjMzqzAn\njio3ahT8+9+pwfyBB2CVVfKOyMy6OyeOGnDWWanksfvu6VHdFVfMOyIz687cj6NGRKTSx/33w+TJ\nsPzyeUdkZrWqo/04nDhqSERq93jsMZg0Cfr2zTsiM6tFThw1HH97RKQnrp59Fu6+G5ZdNu+IzKzW\nOHHUcPzttWhR6uvxyiswYQIsvXTeEZlZLXHiqOH4O2LhQjj6aJg9G8aPhyWXzDsiM6sVThw1HH9H\nLVwIRxwB770Ht90GSyyRd0RmVgs8yGE31rMnXHttqqo66KDU38PMrNzKOQPgGElzJD1XsG4zSX+R\n9KykOxuniC1x7HBJL0p6WdKocsXYFfTuDePGwVJLOXmYWWWUs8QxFhhetO4q4EcRsSlwG3B68UGS\negKXZcduSJpudoMyxlnzeveGG25IJY8DD4RPP807IjPrysqWOCJiCjCvaPV62XqAycBBJQ4dCkyP\niBkRMR+4Edi/XHF2Fb17w/XXwzLLOHmYWXlVuo3jeUmNSeAQoNSA4asDMwvez8rWWQsak8eyy8IB\nBzh5mFl5VHqsqmOASyWdBdwJfFZinzY9JlVfX//5cl1dXbefV7ix2urww1PyuO02P6pr1t01NDTQ\n0NDQaecr6+O4kgYB4yNikxLbBgPXRcTWReu3AeojYnj2/gxgUUT8osQ5uvXjuM1ZsCA9qvvuu3D7\n7U4eZrZYTT2OK2mV7N8ewE+B35fYbSqwnqRBkvoAI0ilE2uDXr3gj3+EFVaAr33N1VZm1nnK+Tju\nOOBRYH1JMyUdQ3pC6iXgBWBWRFyd7TtQ0gSAiFgAnAzcA/wduCkiXihXnF1Zr15w3XVpGPb994dP\nPsk7IjPrCtxzvBtYsACOPBL+9S+4447U58PMuq+aqqqyfPTqlXqYr7KKSx5m1nFOHN1Er15wzTUp\neey3H3z8cd4RmVmtcuLoRhpLHv37p5KHk4eZtYcTRzfTs2cqeQwYAPvsAx99lHdEZlZrnDi6oZ49\n4eqrYdAgGD4c3n8/74jMrJY4cXRTPXvCVVfBxhvD7rvDvOJRxczMmuDE0Y316AG/+x1suy3suiu8\n/XbeEZlZLXDi6OYk+PWvYY89YOedYc6cvCMys2pX6UEOrQpJcMEFqWNgXR3cdx8MHJh3VGZWrZw4\nDEjJ4+yz07zlw4al5PGlL+UdlZlVIycO+w8//nEaSbcxeay9dt4RmVm1ceKwL/j+91PJo64OJk+G\nwYPzjsjMqokTh5X07W+n5LHzzjBpEmy0Ud4RmVm1cOKwJh1zTEoeu+0Gd98Nm22Wd0RmVg2cOKxZ\nhx+ekseee8Jdd8GWW+YdkZnlzYnDWnTwwSl57LVXmoZ2u+3yjsjM8lTOGQDHSJoj6bmCdUMlPSFp\nmqS/StqqiWNnSHo22++JcsVorbfvvmk2wf33T09bmVn3VbYZACXtCHwIXBsRm2TrGoALIuIeSV8F\nfhQRO5c49lVgi4h4p4VreAbACnvwQTjkEBg9OiUTM6s9VTsDYERMAYqHzpsN9MuWlwfeaOYU7f5Q\nVj7DhsGECXD88XDDDXlHY2Z5qHQbx4+BhyVdREpa2zaxXwCTJS0E/hARV1YqQGvZVlul/h3Dh8MH\nH8CJJ+YdkZlVUqUTx2jguxFxm6RDgDHA7iX22z4iZktaBbhX0otZCeYL6uvrP1+uq6ujrq6u86O2\nL9h441RttdtuaT6P00/POyIza0pDQwMNDQ2ddr6ytXEASBoEjC9o43g/IvpmywLejYh+TZ8BJJ0N\nfBgRF5fY5jaOnM2alebzOOggOPfcNOaVmVW3qm3jaMJ0ScOy5V2A/yveQdLSkpbLlpcB9gCeK97P\nqsMaa8BDD8HEifC978GiRXlHZGblVs6nqsYBw4CVgTnAz0gJ4HJgCeAT4KSImCZpIHBlROwtaW3g\n1uw0vYDrI+KCJq7hEkeVePfdNIf5uuummQV7uYeQWdXqaImjrFVV5ebEUV0++ggOPBCWWw6uvz51\nGjSz6lNrVVXWhS2zDNx5J0TAfvulRGJmXY8Th3WqJZaAm26CAQPS+FbvvZd3RGbW2Zw4rNP16gVj\nx8KQIWlY9rlz847IzDqTE4eVRY8ecOmlaWDEnXZKj+2aWdfgZ1+sbCQ47zzo1w922AHuuQfWXz/v\nqMyso5w4rOxOPx1WXjlNRTt+vOf0MKt1rqqyijj6aPif/0lVVx6W3ay2OXFYxey/P9x8M4wcCf/7\nv3lHY2bt5aoqq6hhw2DSJNh7b3j7bY+sa1aLnDis4jbfPI1vtcce6VHdn/zEgyOa1RIPOWK5mT07\nzelRVwe//nV6hNfMys9jVdVw/JYGR9x3X1hrrdRpsHfvvCMy6/o8VpXVtOWXT/073nsvNZ57fCuz\n6ufEYblbemm49VZYZZU0KdQ77+QdkZk1x4nDqkLv3qmqarvt0hAlb7yRd0Rm1hQ/VWVVo0cPuPDC\nVPJoHKJk8OC8ozKzYmUrcUgaI2mOpOcK1g2V9ISkaZL+KmmrJo4dLulFSS9LGlWuGK36SDBqFPz0\np6nPxxNP5B2RmRUrZ1XVWGB40bpfAmdFxBDSVLK/LD5IUk/gsuzYDYGRkjYoY5xWhY49Fq68MnUU\nnDAh72jMrFDZEkdETAHmFa2eDfTLlpcHStVkDwWmR8SMiJgP3AjsX644rXrtsw/cdRccdxyMHp13\nNGbWqNJtHD8GHpZ0ESlpbVtin9WBmQXvZwFbVyA2q0Jbbw0PPpg6Cr7xBpx1lnuZm+Wt0oljNPDd\niLhN0iHAGGD3on3a1KOvvr7+8+W6ujrq6uo6GKJVm8GD4dFH08i6b7wBl1+eZhk0s9ZpaGigoaGh\n085X1p7jkgYB4yNik+z9+xHRN1sW8G5E9Cs6ZhugPiKGZ+/PABZFxC9KnN89x7uRDz6Agw6CJZeE\nG29M/T/MrO1qref4dEnDsuVdgP8rsc9UYD1JgyT1AUYAd1YqQKteyy2X2jz69YNdd4V//SvviMy6\np3I+jjsOeBRYX9JMSUcDJwC/lPQ0cF72HkkDJU0AiIgFwMnAPcDfgZsi4oVyxWm1pU8fuPbaNDDi\n9tvDq6/mHZFZ9+NBDq1mXXYZXHBBKoUMGZJ3NGa1oyKj42b9KAYBi4DXIuLF9l6wMzlx2C23wLe/\nDddfn8a5MrOWlS1xSPoy8ANgL1J/izcBAasBawB3Ab+OiBntvXhHOXEYwJQpcPDBcPHFcMQReUdj\nVv3KmTj+BFwJNGQd8Qq39QZ2Bo6LiEPbe/GOcuKwRs8/nx7X/c534PTT3dfDrDmeyKmG47fONWtW\nSh477gi/+Y37epg1pZwljp2KVxUsR0Q81N6LdhYnDiv23ntwyCGwxBIwbhwsu2zeEZlVn3Imjrso\n3Yt7U2CNiOjZ3ot2FicOK2X+fPjWt+Dpp9MTV6utlndEZtWlYlVVkrYHziINTnh+RIxv70U7ixOH\nNSUCzj8frroqja670UZ5R2RWPcqeOCTtBvw0e3t+RNzb3ot1NicOa8kf/winnpqqrXbZJe9ozKpD\nOauq9gF+ArwL/Hc2THpVceKw1mhogBEj0uyCRx6ZdzRm+Stn4lhEGtL8mRKbIyL2a+9FO4sTh7XW\nCy+kSaG++U342c/8uK51b+VMHMMK3xZtjoh4sL0X7SxOHNYWc+akyaE22giuuCKNe2XWHZUzcVwB\n/BmYHBEftPcC5eTEYW310Udw+OFpiPZbboHll887IrPKK+ew6mOAzYGJku6XNErSZu29kFk1WGaZ\nlDA23hh22AFeey3viMxqT2sHOVwZ2AMYTurHMQ34c0T8qbzhtRiXSxzWbpdcAhddBHfeCV/5St7R\nmFVOLkOOSNoS2DMizm/vhTuDE4d11G23wYknwpgxqf3DrDso+wyAkr4vqZ+S0ZKeAlZqKWlIGiNp\njqTnCtbdKGla9npV0rQmjp0h6dlsvyfa/KnMWumAA2D8eDjhhFQC8d8hZi1rTQfAZyNiU0l7At8i\n9R6/LiKanTpH0o7Ah8C1jXOOF22/iDTn+Hkltr0KbBER77RwDZc4rFO89hrsuy9stx389rfQu3fe\nEZmVTyXmHG88+d6khPG31pw46zA4r+QJJQGHAuNacV2zsltrLXj4YZg5M/X3ePfdvCMyq16tSRxP\nSppEmtDpbkl9STMBdsSOwJyIeKWJ7QFMljRV0vEdvJZZq/TtC3fcARtsANtuC6809e006+aaTBzZ\nZE0AxwBnAFtFxMdAb+DoDl53JHBDM9u3z6rCvgp8J6v2Miu7Xr3SXB6nnJIe13344bwjMqs+zU11\n8xdJb5A6Ad4dEfMAIuJt4O32XlBSL+AAoMkHICNidvbvXEm3AUOBkmNl1dfXf75cV1dHXV1de0Mz\n+9xJJ8E668CBB6Ypab/xjbwjMmu/hoYGGhoaOu18zTaOZ/OODwf2JM0zPoWUSB6MiH+3eHJpEDC+\nsHFc0nBgVETs3MQxSwM9I+IDScsAk4D/iohJJfZ147iV1fPPp0bzww6Dc86BHq2p3DWrcpWcj6MP\nqW1iODAMmBsRezez/7hsv5WAfwI/i4ixksYCf4mIKwr2HQhcGRF7S1obuDXb1Au4PiIuaOIaThxW\ndnPnwte+BquvDtdcA0stlXdEZh2T25zjklaPiDfae+HO4MRhlfLpp3DssTB9empAHzAg74jM2q8S\nEzk9V7QqgPeAqcB5WZtHLpw4rJIi4NxzYfTo1Glw003zjsisfSqROC4EFpCeghLwdWBp4C3S00/7\ntvfiHeXEYXm48cb01NXYsR6mxGpTJRLHtOJe4o3rJD1Xqld4pThxWF4eeyw9cfXDH6apaT0xlNWS\nSvQc7ylp64ILDi04bkF7L2xWy7bZJiWP66+Ho45KbSBm3UVrShxbAWOBZbNVHwDHAs8De+c5tLpL\nHJa3jz5KiWPWrDTSrhvNrRZU8nHcfgAR8V57L9bZnDisGixatLjR/PbbPbeHVb9KtHGsABwJDGJx\nT/OIiO+296KdxYnDqsktt8C3vgWXXw6HHpp3NGZN62jiaG7IkUYTgb8Az5IGNxTpkVwzK3DQQWmY\nkv33h7/9Derr3dPcuqbWlDieioiqLHy7xGHVaM6c9MTVgAGpp/myy7Z8jFklVeKpqhsknSBpNUkr\nNr7ae0Gzrq5/f7j/fujXD7bfPk0SZdaVtCZxfApcCDwGPJm9ppYzKLNat8QSqbH8qKPS3B6PPJJ3\nRGadpzVVVa+S5uL4V2VCaj1XVVktuPtuOPJI+PnP4Zhj8o7GrDJVVS8Dn7T3Ambd3fDhMGVKShw/\n+AEscLdZq3GtKXHcDmwEPAA0zsHhx3HN2mjePBgxIg2WeOONsNJKeUdk3VUl+nEclS027ihS4rim\nvRftLE4cVmsWLIAzzoBbb009zT3CruUht/k4qoETh9WqG26A730Pfvc7OOSQvKOx7qZsbRySJkg6\nJJvKtXjb0pJGSJrYzPFjJM0pnM9D0o2SpmWvVyVNa+LY4ZJelPSypFFt/VBm1e6ww2DSJDj9dDjz\nTFi4MO+IzFqvyRKHpFWBk4GDgYXAbFI11QBSj/ObgMsjYm4Tx+8IfAhcW2rodUkXAe9GxHlF63sC\nLwG7AW8AfwVGRsQLJc7hEofVtLlzU7vHkkumkXZXWCHviKw7qEhVlaQBwFrZ2xkRMaeVwQ0Cxhcn\nDkkCXgN2johXirZtC5wdEcOz9z8GiIiflzi/E4fVvPnzU8ljwoQ0SOJGG+UdkXV15ayqWk/SDgAR\n8VZEPB4RjwPrSVqnvRfM7AjMKU4amdWBmQXvZ2XrzLqk3r3hkkvgrLNg551To7lZNWuuH8clwPsl\n1r+fbeuIkaSpaEtxEcK6pSOPhIkTU6P52Wen4drNqlFzo+P2j4hni1dGxLOSvtzeC0rqBRwANDVw\n4hvAmgXv1ySVOkqqr6//fLmuro66urr2hmaWuy23hL/+NT1pNW0aXHddGvPKrCMaGhpoaGjotPM1\n1zg+PSLWbeu2ov0GUdTGIWk4MCoidm7imF6kxvFdgTeBJ3DjuHUzn32W5jOfPBnuuAPWXz/viKwr\nKeeQI1MlnVDigseTBjpsKbBxwKPAYEkzJR2dbRoBjCvad6CkCQARsYD0NNc9wN+Bm0olDbOurE8f\nuOyy1Gi+444wfnzeEZkt1lyJYwBwG/AZixPFFsASwAERMbsiETbDJQ7rDh57LFVdHX10avvo2TPv\niKzWlfVx3Oyx2Z2BjUmN1s9HxP3tvVhnc+Kw7mLOnNTfY6mlUn+PFT0jjnVAOR/HPTCS+4E/RsRv\nqylpmHUn/fun9o4NN0wN6NNKjrlgVhnNtXGcVbA8udyBmFnzevWCiy9Ow7PvsUealtYsD62ZjwPS\nUCNmVgUOPRQefBAuuABOOgn+/e+WjzHrTM0ljqUkfUXSFoXL2b9N9cEwswrYcEN44gl46y0YNgxm\nNdnTyazzNfdUVQNFc3AUbm+qH0YluXHcursI+OUv05AlN9yQhiwxa4nn46jh+M06y+TJcMQRcNpp\ncOqpIFcuWzPKljgkbQXMauyvIembwEHADKA+It5p70U7ixOH2WKvvw4HHQSDBsGYMbDccnlHZNWq\nnD3HryCbY1zSTsDPgWtIgxxe0d4Lmll5fOlLMGUKLL88bL01vPhi3hFZV9Vc4uhRUKoYAfwhIm6J\niJ8C65U/NDNrqyWXhCuvTONc7bQT3Hxz3hFZV9Rc4ugpqXe2vBvwQMG25kbVNbOcHXcc/PnPMGoU\nfP/7adBEs87SXOK4GXhQ0p3Ax8AUSBM8Ae9WIDYz64AttoAnn4RXX02P7L7+et4RWVfRXOL4GvBD\n4Gpgh4honFZGwClljsvMOsEKK6TpaA88EIYOhbvvzjsi6wqae6pqWkQMqXA8beKnqsxa76GH4LDD\n4JhjPMpud1fOx3FnAb+i9HAjERG/au9FO4sTh1nbzJkDI0dCjx6pw+Cqq+YdkeWhnI/j9gSWA5Yt\n8fIT4mY1qH9/uPde2Gab1Aby8MN5R2S1qGxVVZLGAHsD/yyaOvYU4CRgITAhIkaVOHYGqb/IQmB+\nRAxt4houcZi104QJqdrqRz9Kj++6t3n3Uc6qqo4mjh2BD4FrGxOHpJ2BM4G9ImK+pFUiYm6JY18F\ntmipd7oTh1nHvPZaml1w9dVh7NjUedC6vnJWVe3W3pMCRMQUYF7R6m8DF0TE/GyfLySNAv77x6zM\n1lor9TZfYw1PEGWt12TiiIi3y3C99YCdJD0mqUHSlk1dHpgsaaqk48sQh5llllgCfvtbOP/8NEHU\nFVekUXfNmlLpHuC9gBUiYptsEMU/AWuX2G/7iJgtaRXgXkkvZiWYL6ivr/98ua6ujrq6us6P2qwb\nGDECNtssTRT1wAPwhz9A3755R2WdoaGhgYaGhk47X1mHVZc0CBhf0MbxZ+DnEfFg9n46sHVzpRtJ\nZwMfRsTFJba5jcOsk33ySRqm5P774U9/giFV3ZvL2qOcbRzlcDuwC4CkwUCf4qQhaWlJy2XLywB7\nAM9VOE6zbmuppVJp49xzU9XV737nqiv7T2VLHJLGAY8CgyXNlHQ0MAZYW9JzwDjgyGzfgZImZIcO\nAKZIehp4HLgrIiaVK04zK+3rX4dHH02j7R56KLz3Xt4RWbXwDIBm1qxPP02zCt59N9x0U3r6ympb\nrVVVmVmNWXJJuPxy+MUvYK+94De/cdVVd+cSh5m12iuvpKev1lwzTU+7wgp5R2Tt4RKHmVXMOuvA\nI4+kaWqHDIHHH887IsuDE4eZtckSS6Tqqksugf32g4svdtVVd+OqKjNrtxkz0tNXq6ySxrpaeeW8\nI7LWcFWVmeVm0KA0QdQGG6SqqwceyDsiqwSXOMysU0yaBEcdtXiGwd69847ImlK2YdVrgROHWXWZ\nMwe++U14//00w+CgQXlHZKW4qsrMqkb//jBxIhx0EAwdCjffnHdEVg4ucZhZWUydmuY3r6tLT2At\ns0zeEVkjlzjMrCptuSU89RT8+99p+Zln8o7IOosTh5mVzXLLwbXXwplnwm67wWWXuc9HV+CqKjOr\niOnTU9UBZxAJAAAONUlEQVTVwIFpuJKVVso7ou7LVVVmVhPWXTcNVzJ4MGy+OXTihHRWYS5xmFnF\n3XMPHH10etXXu89HpbnEYWY1Z889Ydo0ePpp2H57ePnlvCOytijnDIBjJM3JZvsrXH+KpBck/U3S\nL5o4drikFyW9LGlUuWI0s/z07w933ZU6DG63HYwe7YbzWlG2qipJOwIfAtdGxCbZup2BM4G9ImK+\npFUiYm7RcT2Bl4DdgDeAvwIjI+KFEtdwVZVZF/D883DYYakd5Ior3HBeblVbVRURU4B5Rau/DVwQ\nEfOzfeZ+4UAYCkyPiBnZfjcC+5crTjPL30YbwRNPpCFKNt8c7rsv74isOZVu41gP2EnSY5IaJJWa\nvXh1YGbB+1nZOjPrwpZYIs3tMWZMqr467bTUedCqT68crrdCRGwjaSvgT8DaRfu0qe6pvr7+8+W6\nujrq6uo6GKKZ5Wn33VOj+fHHw9Zbp8ESN9ww76hqW0NDAw2d+PxzWR/HlTQIGF/QxvFn4OcR8WD2\nfjqwdUS8XXDMNkB9RAzP3p8BLIqILzSku43DrOuKgKuugjPOgHPOgW9/G9TuWnkrVLVtHE24HdgF\nQNJgoE9h0shMBdaTNEhSH2AEcGdlwzSzvEmp1PHII6n6at994Z//zDsqg/I+jjsOeBQYLGmmpKOB\nMcDa2SO644Ajs30HSpoAEBELgJOBe4C/AzeVeqLKzLqH9deHRx+FTTZJDecTJ+YdkbnnuJnVjIaG\n1HC+115w0UUeqr29aq2qysys3erq4Nln4eOPU+njscfyjqh7conDzGrS//4vfOc7cOKJcNZZHu+q\nLVziMLNu6eCD02O7U6fCttvCC24JrRgnDjOrWautBhMmwHHHwY47wqWXwqJFeUfV9bmqysy6hJdf\nhm98I806OHYsrLFG3hFVL1dVmZkB660HDz8Mw4bBV74C48blHVHX5RKHmXU5Tz4JRxyRnry6/HJY\nccW8I6ouLnGYmRXZYgt46ilYdVXYbDOYNCnviLoWlzjMrEubPBmOOSYNWfKLX8Cyy+YdUf5c4jAz\na8Zuu6VOgx9+mEofU6bkHVHtc4nDzLqNO+5Io+yOHAnnnQdLLZV3RPlwicPMrJX23z+VPmbNSk9e\nPfFE3hHVJpc4zKxbuukm+O5309DtP/sZ9OmTd0SV4xKHmVk7jBgBzzyTSiBbbZWWrXWcOMys2xow\nILV7/PCHqRH9vPNgwYK8o6p+rqoyMwNmzoRjj4V58+Caa7r2POdVW1UlaYykOdlsf43r6iXNkjQt\new1v4tgZkp7N9nHzlZmV3Zprwj33pAETd9opTRS1cGHeUVWnspU4JO0IfAhcGxGbZOvOBj6IiF+1\ncOyrwBYR8U4L+7nEYWad7h//gKOPhvnz04CJ66+fd0Sdq2pLHBExBZhXYlNrg233hzIz64i114YH\nHkj9PbbfHi680KWPQnk0jp8i6RlJoyUt38Q+AUyWNFXS8ZUMzswMoEcPOOWU1Ndj4sSUQDxZVNKr\nwtf7PXBOtnwucDFwbIn9to+I2ZJWAe6V9GJWgvmC+vr6z5fr6uqoq6vr1IDNrHtbe2247z74wx/S\nZFGnnZZevSr969kBDQ0NNDQ0dNr5yvpUlaRBwPjGNo7Wbiva72zgw4i4uMQ2t3GYWcXMmJEaz997\nL7V9bLxx3hG1T9W2cZQiabWCtwcAz5XYZ2lJy2XLywB7lNrPzKzSBg2Ce+9Nvc133hnOPz81oHc3\n5XyqahwwDFgZmAOcDdQBm5PaMF4FToyIOZIGAldGxN6S1gZuzU7TC7g+Ii5o4houcZhZLl5/PSWQ\nuXPh6qth003zjqj1OlricAdAM7N2ikhVVqNGwcknwxln1MaYVzVVVWVm1pVIaZKoadPS01dDh6bl\nrs4lDjOzThAB110HN98Md96Zkkq1clVVDcdvZl1PRHUnDXBVlZlZVan2pNEZnDjMzKxNnDjMzKxN\nnDjMzKxNnDjMzKxNnDjMzKxNnDjMzKxNnDjMzKxNnDjMzKxNnDjMzKxNnDjMzKxNnDjMzKxNypY4\nJI2RNEfScwXr6iXNkjQtew1v4tjhkl6U9LKkUeWK0czM2q6cJY6xQHFiCOBXETEke91dfJCknsBl\n2bEbAiMlbVDGOLuEzpyIvtb5Xizme7GY70XnKVviiIgpwLwSm1oaO3IoMD0iZkTEfOBGYP/Ojq+r\n8f8Ui/leLOZ7sZjvRefJo43jFEnPSBotafkS21cHZha8n5WtMzOzKlDpxPF74MvA5sBs4OIS+3hm\nJjOzKlbWGQAlDQLGR8Qmrd0maRugPiKGZ+/PABZFxC9KnMNJxsysHToyA2CvzgykJZJWi4jZ2dsD\ngOdK7DYVWC9LLG8CI4CRpc7XkQ9uZmbtU7bEIWkcMAxYWdJM4GygTtLmpOqoV4ETs30HAldGxN4R\nsUDSycA9QE9gdES8UK44zcysbcpaVWVmZl1PVfYcb6kDoKTTCjoRPidpQeMTWpJmSHo22/ZE5aPv\nXK24FytLulvS05L+Jumo1h5bazp4L7rb92IFSbdlTzA+Lmmj1h5bazp4L7rM96JUp+sS+1ya3adn\nJA0pWN+270REVNWLVD01HRgE9AaeBjZoZv99gMkF718FVsz7c1TqXgD1wAXZ8srA26QqyDbdx2p/\ndeRedNPvxYXAWdny+o3/j3TT70XJe9EFvxc7AkOA55rYvhcwMVveGnisvd+JaixxtLUD4GHAuKJ1\nXaXRvDX3YjbQN1vuC7wdEQtaeWwt6ci9aNSdvhcbAA8ARMRLwCBJq7by2FrS3nuxSsH2LvG9iKY7\nXTfaD7gm2/dxYHlJA2jHd6IaE0erOwBKWhrYE7ilYHUAkyVNlXR82aKsjNbciyuBjSS9CTwDfK8N\nx9aSjtwL6H7fi2eAAwEkDQXWAtZo5bG1pCP3ArrW96IlTd2rgU2sb1JFH8dtpba01u8LPBwR7xas\n2z4iZmd/Udwr6cUsE9ei1tyLM4GnI6JO0jqkz7xZmePKQ7vvRUR8QPf7Xvwc+I2kaaTH3qcBC1t5\nbC3pyL0A2CEi3uwi34vW6JTSVTWWON4A1ix4vyYpA5bydYqqqSLrJxIRc4HbSMWwWtWae7EdcDNA\nRLxCqrNdP9uvtfexFnTkXnS770VEfBARx0QaTPRIYBXgldYcW2Paey/+kW17M/u3K3wvWlJ8r9Yg\n3as2fyeqMXF83gFQUh9SB8A7i3eS1A/YCbijYN3SkpbLlpcB9qB0J8Na0Zp78SKwG4Ck/qQfyn+0\n8tha0u570R2/F5L6ZdvIqmAejIgPW3NsjWn3veiC34uW3AkcCZ+P0PFuRMyhHd+JqquqiiY6AEo6\nMdv+h2zXrwH3RMQnBYf3B26TBOmzXR8RkyoXfedq5b34b2CspGdIfwj8KCLeASh1bB6fozN05F5I\nWhu4tZt9LzYErlYaludvwLHNHZvH5+gMHbkXdLHfC5XudN0b0n2IiImS9pI0HfgIODrb1ubvhDsA\nmplZm1RjVZWZmVUxJw4zM2sTJw4zM2sTJw4zM2sTJw4zM2sTJw4zM2sTJw6rWZIWavHQ+n+StFTe\nMTVH0lGSfpstnyjpGwXrVyva96Zs2JTmztcgaYts+RBJf5d0XzP795c0seOfxLo7Jw6rZR9nw0hs\nAnwGfKucF5PU0Q6zn3eayjpkXZe9/SZpoLnG66wLLJMNm9LS+RrPeSxwXETs2uTOqZfwPElfaU/w\nZo2cOKyreBhYV2nSntuziWr+ImkTgGyynr5K3i74a/9aSbtK6iHpQklPZMeekG2vkzRF0h3A84UX\nlNRT0tVZiedZSd/L1jdIuqSgNLRV4yEFx9ZLOlXSQcCWwPWSnpK0JGkMtjsL9v2dpL8qTU5VX/S5\nJelnwPbAGEm/lLSWpIckPZm9ti3Y/05gZAfvtXVzThxW87KSwHDgWeAc4MmI2Iw0Wu612W6PADsA\nG5EG+9shW78N8ChwHGnsnqGkge6OlzQo22cI8N2IWL/o0psDAyNik4jYFBibrQ9gqYgYApwEjCkR\ndgAREbeQxgo6LCK+EhGfkpLA1IJ9fxIRWwGbAcMak2HjeSLinIJz/Aj4J7B7RGxBSkKXFuz/BGmM\nN7N2q7qxqszaYCmlobIBHiL9QD9ONvdCRDwgaaVsILsppB/M14DfAydIGgjMi4hPJO0BbCLp4Ox8\nfYF1gQXAExHxWonrvwKsLelSYAJQOM7RuCyGKVlJp18Ln6VwuOu1SJNSNRqRDc7XC1iNNDFRqcH4\nGs/RB7hMaXj9hcDggn1mk2Z6M2s3Jw6rZZ9kf9V/LhuwrnjOgSAllpOBGcBPgAOAg7P1jU6OiHuL\nzldHGhDuCyLiXUmbkko73wIOZfEAesUWtfBZigeNU3b9LwOnAltGxHuSxgJLtnCOHwCzI+IbknoC\nnxad1wPUWYe4qsq6minA4fD5j/7ciPgwImaR5iFfNyJeJbWJnMbixHEPcFJjA7ikwUozTDZJ0kqk\nOc1vBc4iVWlB+nEeke2zA6kK7IPiw1mc4D5g8ZS3kEpFjU9Z9SUlrveVhor/aivuQV/grWz5SNKI\np41Wy85v1m4ucVgtK/WXcz2pkfgZ0g/uNwu2PcbiP5YeJg3D/nD2/ipSFc5TSsWWf5JKJYVPLhVb\nnTSMe+M5f1wQ16eSniL9P3ZMwfoosXw18D+SPiZNRvUwqcH8yYh4JquOe5E0vWdjvM35HXCLpCOB\nu4EPC7YN5T9LWWZt5mHVzTqZpAeAUyPiqXYevzbw24jYu3MjA0nXAxdFxLQWdzZrgquqzKpMRPwD\n+KClDoBtJWlVYHknDesolzjMzKxNXOIwM7M2ceIwM7M2ceIwM7M2ceIwM7M2ceIwM7M2ceIwM7M2\n+f/STtLiZhT/9gAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x5abb490>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAasAAAEZCAYAAAApEwoTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XucHXV9//HXO5t7SDYJgQgIrBZEQEAqjfQn1i3ww4go\nrZdSbMVbEfVHta1WitZm1bZesWqxSv1xUYvQi1LwVwUvZRVvYBQCKKFECAQSQsj9trnt5/fH93t2\nZ0/O2T2b5OzO2X0/H4/zODPz/c7Md77nO/OZ75w5cxQRmJmZldmE0S6AmZnZUByszMys9ByszMys\n9ByszMys9ByszMys9ByszMys9FomWEnqlfTsPHydpA83cV2dklY0a/mtQNJxku6RtEnSpSOwvvsl\n/U6z11NjvddKWifpp6Ow7o7crhvaD4ebv5VI6pL0ldEux1AkvVHSHQ3k+7ykvx6JMlWtd0B7lvR2\nSavzfjxX0mZJHUMs46icTyNR5kYN2eglLZe0LRf+yVwZM0aicIOI/DogioFwpDU78O6H9wLfi4hZ\nEXHlgVxwrW2OiOdFxA8O5HoaKMeLgbOBwyPi9BrpDR2Y9jV/A8tbLunMA7W8Ya57pIPHmPrBZ0S8\nPSL+diTXWWjPR0TE6ZImAVcAZ+X9eF1EzIyI5YMtJyIey/lK9Zk0coYWwHkRMRP4TeA0YK8zBkkT\nD3DZhjLsqD9EGff5LGIsnukCRwO/qpc4Rrb5aGB5RPSMdkHqCPavXbYdwLK0jDHSNvdFpT1vz+PP\nAKYCD4xekQ6giBj0BTwCnFkY/wRwSx7uBd4BPAT8Ok+7OI+vBW4GDivM2wv8KfBrYA3wcUCF9DeT\nDpDrgFuBo6rmfXYevhb4cCHtPOAeYD3wI+CkQtpyUi/hXmA7MKFq+36Ql70F2Ay8FugEVgB/AawG\nVgJvLMxzHfB54Jt5vrOAbuAthTxvBO7IwwL+IS9rYy7LicBbgZ3Ajrzum+t8Bp8BHsvzLgbOKKQt\nyNM2Ak8CV9RZxmzg/wFP5fr9BukMrFbe/wZ25/raBBxbY5vPBF4O3J3X/RiwqGo5ZwA/zp/LY8Ab\ncvvYa5vz53RWHp4CfBp4Ir/+AZic0zqBx+t9NjW25XDgFlJ7fAj4kzz9LXn7dudyVJf9+Kr0dXl6\nO/DlXI/Lgffnz7de/rp1BHSQ2t6EGuX+CrAH2JaX955C/ouAR0n70PsK83QB/5Hn3Zi38ToG7iud\nwIrC+GW5PjcBS/PnujB/Pjvzuu+uU7d/BSzL8/4S+L2q9v9D0vFiHfAwsLCQ/izg+3nebwP/CHyl\nznoqn/nleZsfAV43yP54Zv48uklt737gFVX5v5DXuynnO6reZ0Jh32bo/fqEwjo+nIcfAF5eWN7E\nvB3Pz+On07+f3AO8ZIj2/DVS+3sY+NM67fmruS568/h3axxHp5F6XsuBDcAdpH1vQB2Q2vzVpH3t\nceDDhbShPue5pOP1Ezn963n6/aROUCXfJOBp4JS6214vobCQR+g/iByZV/LBwobfRjoQTsmNZA3w\nfGAy8Fng+1UB53s5/5HAg4VGcD7pYHIcqcf3fuBHVfPuFayAU3Nj+a3ceC7KZZ5UOAj+AjgCmFJn\nG/uWXdg5dpF2/jbgZcBWoL3QEDcAv104uN4OvLlqZ6006peSAsqsPH4c8IzCtnxoiM/gj4A5uV7+\nAlhF/8H7J8Af5eHpwAvrLGMu8PukM62DgH8DbhpkndXbU2ubXwKcmMdPIgXL8/P40aQDwQW5DueS\nG2KtbaZwUgR8iLTzzsuvH1XyD/XZ1NiOHwBXktrjKaSd/Hdz2hsqn1GdefdKJwWqm4AZeRsfrNRT\nnfyD1VEHdYJVdZ1U5b8q1//JQA9wXE7vIgWYV+bxqdV1TSFYkdrhY/S3xaPo38cWAV8eol2+pjDv\nH5AOjvML7X8n6SAq4G3AE4V5fwJ8knSQenFuKzXXV/jMK/l/J6/rOXXa5kxSEP0rUmD43bz8Yv5N\npJOpyaQTozuq6rgYrPr2BfZhvwY+APxLYXkvB36Zh48gHaQX5vGz8/i8GvUwAfg56crWRFLA/zVw\nTq32R2qf1dtSPI5+jnRielhe9um5PgbUAam9f54U3A4B7gTe2uDn/F/ADaSANxF4cZ7+l8CNhXzn\nA0sGbW+DJeaFLCdF5vV5+EryQT9vUGch79XARwvjM/KGHFXIf04h/e30R/xvMfDgOIF0EDqyRiUX\nG8Ln2fvAt7RQKY8wyJl39bILO8e2qg95NbCg0Nivq1rGYMHqTNJB7YXs3bPrC7yNvkhnKCfl4e+T\nDlJ7Ne4hlvF88tl/nfTbGdhTvLZ6m2vM82ngU3n4cuBrdfLttc0MDFbLGHh2dg7wSCOfTdUyjySd\nac4oTPt74Nrqz6hOOQekk4LjDuC5hWlvBW5vZHk16qiDfQtWhxem3Qn8QR7uAroHq2sGBqtjct2d\nRT65K+Trok5PZ5Btu5v+QPlG4KFC2vRc9kNJQXEXMK2Qfn299dEfrIr5/xX46zx8XbFtkoLfqqpl\nfJXcq835v1pIm5HbyRG1PhPqB6uG9utcz5uAqYVtrZT9MqqCNOmq0kU16uGFwKNV0y4HrqnTXmtt\nSy/wbNLxdRuFq1C15gPmk06IphbSLwT+u4HP+TDS1YG9TiRJPcTNwEF5/D+A9wzWvhr9zur8iJgT\nER0RcWlE7CikF++aO4x0eSLNGLGVdPnliDr5H8uFhnQW8BlJ6yWtz/NRNW8tRwPvrsyX531mYbnV\n62zU2ojoLYxvI/VIINVJw8uMiP8mBfnPAaslXSVpZqPzS3qPpF9J2pC3r53U44B0RvMc4AFJd0l6\neZ1lTM/rXS5pIynItQ9xx09UjQ/YZkkvlHS7pKckbQAuAQ7OyUeSLgnsi8MptCMGthMY/LOpXs66\n3A6LyxqqTdUzj3RmX122ussboo721ZOF4eptf7zRhUTEMuDPSIFptaQbJB3W6PySLpJ0d2G/ex4D\nt62vnBGxLQ8eRPpc1kf/dyswsE5rqZW/UtZg4HYfzt7756P0t6EB+XP7WMfANjakRvfrXM8PAK+U\nNB14BSl4Qjp+vbbq+PUi0vdN1Y4GDq/KezkpMAzXPFLP+9dD5Dua1OZXFdb5BVIPq6Le53wkaf/b\nWL3QiFhJumLyGkmzSZeerx+sIAfii8jiAW0lKSoDkO8aPJh0vbLiqKrhStpjpK7lnMJrRkQMdUvx\nY8DfVc13UET8a50yNstW0hlaxYDGFhH/GBGnASeQgstfNlK2fIfPXwKvjYjZETGHdH1cebnLIuJ1\nEXEI8DHgPyRNq7God+f1LoiIdtLlKbEfX+CTdrj/BJ4ZEbNJjbiyvMeA36gz31Cfx4B2RGonK/eh\nfCuBuZKKB/OjaPyAXl3Op0ln+NVle7xOfqhdR43ud/vSbqvn2Uo6262obpc3RMSLSQelILWhIdct\n6Wjgn4H/A8zN7fJ+GmtPq4A5+cBdUVl/PbXyF9tE9XHoyKoTsaPpP9aIdCCtbMtBpMvUK0n1BYPU\nWdEg+3W1G0g9kvOBX0VE5UTuMVKPsnj8mhkRH6+xjMdIVxiKeWdFxHn1yjeIp0k9pmOGyLeCdDXh\n4MI62yPipAbWsYK0/7XXSf8S8Mek+wR+HBGrBlvYgb5r5gbgTZJOkTSFdMnlpxHxWCHPeyTNlnQk\n8E5Sdx7STvw+SScASGqX9No66ykeZL8IvE3SAiUzJL286gA1lNXUP7DWW3+1e4BXSZom6RhSjycA\nJJ2Wz7Ankc6Ee0jd48q6B7ttfibpEsXTkiZL+htgVl9BpD+WVDnL2ZjX2bv3YjiI9AXsRklzSd9J\nDGc7a23zQaQz3p2SFgCvK6R9FThb0mslTZR0sKRTctpQ23wD8NeS5kmaB/wN6aaBYYmIFaTvvj4i\naYqkk0k38fxLg4t4Enhm/tyIiD2k7/r+TtJB+YD954XlrS7mz2rVUaNBaLjtspZ7gHMlzZH0DFJP\nCgBJz5F0Zt5XdzCwXT4JdAzS855B2o6ngQmS3kTqWQ0pIh4lfdfzQUmTJJ1BuklqKJX8LyZ97/Pv\nlU2pyvdT0n723py/My//xkKecyW9SNJk0g0DP4mIJyJiDSmovV5Sm6Q3U+czGGK/ri7TjaTvuN7G\nwB7EvwCvkHROXt9Upd951uqt3wVslvTefJxpk/Q8SafVrK1B5CsT1wCfknRYXtZv5/oo5ltFuhHl\nU5JmSpog6TfUwG8i87zfAv4pH/MnVc13E+kO83eSvgse1P4GqwE7XUR8j/Rl4tdIZynPAv6wap6b\nSV8S3k26O+2aPO9/ks7qbsyXqe4jfbi11hWV8Yj4OekOsytJXfmHSDdZDOestAv4Uu7mvqa4/Dpq\npf8D6fu51aTr1cUD4izSWeg60vd+T5PunoH0Pd8Jed1fr7GuW/Prf/K820lnWBUvBe6XtDmX4Q+r\nLtNWfJr0BenTpAP4t4bYRqhT5wXvAD4kaRPpc+/rzeYTlHNJPbq1pM/75Jw81Db/Lelgdm9+Lc7T\napVrKBeSekIrga8Df5Mv39TbpqL/Jt3l9qSkp/K0PyWdfT9MunvqetLnDenmoer8deuogW35CClo\nr5f0Fw3kr7U9XwGWkNrOraSDZiXPlLyONaTezjzSZSXoDwRrJS3ea0URvyLdSfYTUmB7HumusMHK\nUhx/Hek7mHWkk5EvDbJd5HWsJ32OXwEuiYj/qbWuiNhFutT2srxtVwKvr8r/VdIJ21rSTVp/XFjX\nxaQe0tOkHtOP6mzXYPt1dZmeJO13v83A/eRxUm/rfaSbfx4j7TN7HZtzgDmP9H3zw3nb/pn+k9eh\n6rx6/D2k4+zPcj18hP4gW8x3EenGi8qd2v9Of29zqHW+nnQ1Yinp2PjOwvb0kPbJjvw+KOUvtw44\nSdeQzn6eqnQZJfUCx0TEw5KeS9rJTwXeHxFXFOZdTvpCcg+wKyIWNKWQZlZ6uWf0lYg4cqi8DS7v\nWuDxiPjAgVie7TtJHwCOjYiLhsrbzB/yXkv67US97t1a0lnq79VIC9JdhuuaVDYzG79K9Rih8Sp/\nHfFmUu9rSE37pXdE3EHqtg+YXEhfExGLSV3EWtygzKziQF4CGuryrzWZpItJlzy/FRE/HCo/NPEy\nYC5QB/CNwe4ckbQI2FJ1GfBh0s0Ce4CrIuKLTSukmZmV3kg/z69RL4qIVfkut+9IWpp7amZmNg6V\nMlhV7rePiDWSbiI9/25AsJLkbryZ2T6IiJb7mqUMTyceUGlKT1qYmYdnkB61c1+tGQd7NEcrvrbt\n3EbPrp5hz7do0aJRL3tZXq4L14XrYvBXq2paz0rSDaSnJMxT+iPDRaTHdhARV+UfKP6M9BuBXknv\nIv2m4VDg6/m3iBOB6yPi280qZ5lce8+1/Nmtf8b0SdOZf9B8Dp1xaHpNP7R/uPCaf9B8Zk+dPdrF\nNjNruqYFq4i4cIj0Jyk88qRgC+lHb+POO37rHbz9tLezoWcDq7eu5qmtTw14/XLNL7l9+e0Dpm3Z\nuYXJP57MTV+4iUOmH8KhMw7tf5/RP14Zbp/STv2HEpiZlVMpv7MazyQxZ9oc5kybw3PnPXfI/Dv3\n7OTmW2/m2Bccy1Nbn2LN1jXpfdsalj+xnKe2DZy2fdd2DplxSF8AO2R6//C86fMGTp9xCHOnzWVC\nC/2XXWdn52gXoTRcF/1cF62vqbeuN5OkaNWyj6Ydu3ewZtuavsC2Ztsant72dN/wmm1r+oe3rmHT\njk3MmTZnQBCrBLV50+el4Rn9w/Omz2P6pOlDF8TMRoUkogVvsBhzwWqsXOIqy+eya88u1m5fOyCA\nrd2+ti/APb39aZ7e9nT/+LanaZvQ1h/Ipg8MZPOmz+PgaQcPHJ9+MJPbJg9dGDPbbw5WI2ywYNWq\n21TRytsQEWzdtXVAAFuzbQ1rt63tC3LF19rta1m7bS1TJ07dK4DNm5beD552cN97Je3gaQczbVKt\nf0Ixs8E4WI0wB6uxIyLYtGPTgABW6cEVg1xxfO22tUzQhBTUcm+tL7Dl4bnT5u41PHvqbNomtI32\nJpuNGgerEeZgNb5FBNt2basZxIrv67avGzC+ecdmZk2Z1RfY5k6bm4LZ1Ll9Qa3yqqTPnTaX9qnt\nLXWjiVk9DlYjzMHK9sWe3j2s71nP2m0DA9m67ev6XpUgVxzfunMrs6fOHhDMql9zps4ZOD5tDnOm\nzmFS26ShC2Y2QhysRpiDlY2kXXt2saFnw14Bbf329f2BrWfdXtM29Gxg+qTpewWwSnCbM63+8Kwp\ns9ybswPOwWqEtWKw6ujo4KmnnqKtrY0ZM2bwspe9jCuvvJIZM2YMyFfmbbDh6Y1eNu/YPKCntr5n\nfV9QW9+zPgW3nqrx7evYtmsb7VPb+4LXgPeqabOnzh4wzZctrR4HqxHWisHqWc96FldffTVnnnkm\nK1eu5KUvfSnnnXceH/nIRwbkK/M22MjZ3bubDT0bBgSxWu8bejbsNX3Lzi3MnDKzL4jNnjo7BbUp\ns/uC215phWlTJ04dMz8DsYFaNVj5CRaj5PDDD2fhwoXcf//9o10UK6mJEyb23co/XHt697Bxx0Y2\n9GzoC3iVoFYZX7V5VQpuhWmVYaAveA14Tdl72pxpc2if0t433j61nWkTpznY2QHlYDXCKj2mFStW\n8K1vfYtXv/rVo1wiG4vaJrT1fU+2L3p29/QFunqvRzc+2hfcNvRsYGNPf3Dsjd6+wNUXxKbsPVxJ\nb5/SPmB41pRZ/omBDTAuLwMeqBO+4VZdR0cHa9euZeLEibS3t3PeeedxxRVXMGXKlKry+TKgtbae\n3T19wavYwysGtEpaJX1jT//wlp1bmDFpBu1T2/cKZJXxwd5nTZnlgFeHLwO2kNGKA5K4+eabOfPM\nM0enAGYjZOrEqUw9aCrzD5q/T/NXbkypFciKw49ufDQFvJ6NbNqxqW94446NbNm5hemTpu8VxAYM\nF6ZVxovDs6bM8vd3JTEug5WZldsETUhBZWo7R7UftU/LKAa8YjDbtGNTX0DbtGMTy9YtY9POTTXz\nbNqxid7o7eupVV6VgDZrctX4lFmccMgJnHb4aQe4RszByszGpGLAo33fl7Nj9w427dg04FUJaH3j\nPRt5YtMTbNq5ia07tzpYNYGDlZnZIKZMnMIhE9Nf5NjoGZc3WJTdWNgGMyunVr3Bwj9xNzOz0mta\nsJJ0jaTVku6rk/5cST+R1CPp3VVpCyUtlfSQpMuaVUYzM2sNzexZXQssHCR9LfCnwCeLEyW1AVfm\neU8ALpR0fLMKaWZm5de0YBURdwDrB0lfExGLgV1VSQuAZRGxPCJ2ATcC5zernGZmVn5l/M7qCGBF\nYfzxPM3MzMapMt663vBtcF1dXX3DnZ2ddHZ2NqE4Zmatq7u7m+7u7tEuxn5r6q3rkjqAb0TESYPk\nWQRsiYgr8vjpQFdELMzjlwO9EfGxqvl867qZ2TD51vV9V11pi4FjJXVImgxcANwy8sUyM7OyaNpl\nQEk3AC8B5klaASwCJgFExFWSngH8DJgF9Ep6F3BCRGyRdClwG9AGXB0RDzSrnGZmVn5+gkUJjYVt\nMLNy8mVAG5bdu3ePdhHMzFqGg9UI6ujo4OMf/zgnn3wyM2fOpLe3d7SLZGbWEnwZcAR1dHQwd+5c\nvvGNbzBv3ry9/iG4oszbYGatrVUvA5bxd1ZNpw8emM8pFg0voEjine98J0cc4d84m5kNx7gMVsMN\nMgfSkUceOWrrNjNrVf7OaoRJLdf7NjMbdQ5WZmZWeg5WZmZWer4bsITGwjaYWTm16t2A7lmZmVnp\nOViZmVnpOViZmVnpOViZmVnpOViZmVnpOViZmVnpjcnHLfkpEWZmY8uYC1b+fZKZ2djjy4BmZlZ6\nDlZmZlZ6TQtWkq6RtFrSfYPk+aykhyQtkXRqYfpySfdKulvSXc0qo5mZtYZm9qyuBRbWS5R0LnBM\nRBwLvBX4fCE5gM6IODUiFjSxjGZm1gKaFqwi4g5g/SBZXgl8Kee9E5gtaX4h3bf0mZkZMLrfWR0B\nrCiMP56nQepZfVfSYkkXj3jJzMysVEb71vV6vaczImKlpEOA70hamntqA3R1dfUNd3Z20tnZ2ZRC\nmpm1qu7ubrq7u0e7GPutqf9nJakD+EZEnFQj7QtAd0TcmMeXAi+JiNVV+RYBWyLiiqrpNf/PyszM\n6vP/WQ3fLcBFAJJOBzZExGpJ0yXNzNNnAOcAde8oNDOzsa9plwEl3QC8BJgnaQWwCJgEEBFXRcQ3\nJZ0raRmwFXhTnvUZwNfzI5MmAtdHxLebVU4zMyu/Mfe39mZmVp8vA5qZmTWJg5WZmZWeg5WZmZWe\ng5WZmZWeg5WZmZWeg5WZmZWeg5WZmZWeg5WZmZWeg5WZmZWeg5WZmZWeg5WZmZWeg5WZmZWeg5WZ\nmZWeg5WZmZWeg5WZmZWeg5WZmZWeg5WZmZWeg5WZmZWeg5WZmZVe04KVpGskrZZ03yB5PivpIUlL\nJJ1amL5Q0tKcdlmzymhmZq2hmT2ra4GF9RIlnQscExHHAm8FPp+ntwFX5nlPAC6UdHwTy2lmZiXX\ntGAVEXcA6wfJ8krgSznvncBsSc8AFgDLImJ5ROwCbgTOb1Y5zcys/EbzO6sjgBWF8cfztMPrTDcz\ns3Fq4iivX/szc1dXV99wZ2cnnZ2d+1kcM7Oxpbu7m+7u7tEuxn5TRDRv4VIH8I2IOKlG2heA7oi4\nMY8vBV4CPAvoioiFefrlQG9EfKxq/mhm2c3MxiJJRMR+dRRGw2heBrwFuAhA0unAhohYDSwGjpXU\nIWkycEHOa2Zm41TTLgNKuoHUU5onaQWwCJgEEBFXRcQ3JZ0raRmwFXhTTtst6VLgNqANuDoiHmhW\nOc3MrPyaehmwmXwZ0Mxs+HwZ0MzMrEkcrMzMrPQcrMzMrPQcrMzMrPQcrMzMrPQcrMzMrPQcrMzM\nrPQcrMzMrPQcrMzMrPQcrMzMrPQcrMzMrPQcrMzMrPQcrMzMrPQcrMzMrPQcrMzMrPQcrMzMrPQc\nrMzMrPQcrMzMrPQcrMzMrPSaGqwkLZS0VNJDki6rkT5H0k2Slki6U9KJhbTlku6VdLeku5pZTjMz\nKzdFRHMWLLUBDwJnA08APwMujIgHCnk+AWyKiA9LOg74XEScndMeAV4QEevqLD+aVXYzs7FKEhGh\n0S7HcDWzZ7UAWBYRyyNiF3AjcH5VnuOB2wEi4kGgQ9IhhfSWq1AzMzvwmhmsjgBWFMYfz9OKlgCv\nApC0ADgaeGZOC+C7khZLuriJ5TQzs5Kb2MRlN3KN7qPAZyTdDdwH3A3syWlnRMTK3NP6jqSlEXFH\nceaurq6+4c7OTjo7Ow9Euc3Mxozu7m66u7tHuxj7rZnfWZ0OdEXEwjx+OdAbER8bZJ5HgJMiYkvV\n9EXAloi4ojDN31mZmQ2Tv7Pa22LgWEkdkiYDFwC3FDNIas9p5Et934+ILZKmS5qZp88AziH1vMzM\nbBxq2mXAiNgt6VLgNqANuDoiHpB0SU6/CjgBuE5SAPcDb8mzzwduklQp4/UR8e1mldXMzMqtaZcB\nm82XAc3Mhs+XAc3MzJrEwcrMzErPwcrMzErPwcrMzErPwcrMzErPwcrMzErPwcrMzErPwcrMzErP\nwcrMzErPwcrMzErPwcrMzErPwcrMzErPwcrMzErPwcrMzErPwcrMzErPwcrMzErPwcrMzErPwcrM\nzErPwcrMzEqvacFK0kJJSyU9JOmyGulzJN0kaYmkOyWd2Oi8ZmY2vjQlWElqA64EFgInABdKOr4q\n2/uAX0TEKcBFwGeGMa+ZmY0jzepZLQCWRcTyiNgF3AicX5XneOB2gIh4EOiQdGiD85qZ2TjSrGB1\nBLCiMP54nla0BHgVgKQFwNHAMxuc18zMxpGJTVpuNJDno8BnJN0N3AfcDexpcF4Aurq6+oY7Ozvp\n7OwcViHNzMa67u5uuru7R7sY+00RDceGxhcqnQ50RcTCPH450BsRHxtknkeAk4DnNTKvpGhG2c3M\nxjJJRIRGuxzD1azLgIuBYyV1SJoMXADcUswgqT2nIeli4PsRsaWRec3MbHxpymXAiNgt6VLgNqAN\nuDoiHpB0SU6/inSn33WSArgfeMtg8zajnGZm1hqachlwJPgyoJnZ8PkyoJmZWZM4WJmZWek5WJmZ\nWek5WJmZWek5WJmZWek5WJmZWek5WJmZWek5WJmZWek5WJmZWek5WJmZWek5WJmZWek5WJmZWek5\nWJmZWek5WJmZWek5WJmZWek5WJmZWek5WJmZWek5WJmZWek5WJmZWek1NVhJWihpqaSHJF1WI32e\npFsl3SPpfklvLKQtl3SvpLsl3dXMcpqZWbkpIpqzYKkNeBA4G3gC+BlwYUQ8UMjTBUyJiMslzcv5\n50fEbkmPAC+IiHV1lh/NKruZ2VgliYjQaJdjuJrZs1oALIuI5RGxC7gROL8qzypgVh6eBayNiN2F\n9JarUDMzO/CaGayOAFYUxh/P04q+CJwoaSWwBHhXIS2A70paLOniJpbTzMxKbmITl93INbr3AfdE\nRKek3wC+I+mUiNgMvCgiVkk6JE9fGhF3FGfu6urqG+7s7KSzs/PAld7MbAzo7u6mu7t7tIux35r5\nndXpQFdELMzjlwO9EfGxQp5vAn8XET/K498DLouIxVXLWgRsiYgrCtP8nZWZ2TD5O6u9LQaOldQh\naTJwAXBLVZ6lpBswkDQfOA54WNJ0STPz9BnAOcB9TSyrmZmVWNMuA+Y7+i4FbgPagKsj4gFJl+T0\nq4C/B66VtIQUON8bEeskPRv4uqRKGa+PiG83q6xmZlZuTbsM2Gy+DGhmNny+DGhmZtYkDlZmZlZ6\nDlZmZlZ6DlZmZlZ6DlZmZlZ6DlZmZlZ6DlZmZlZ6DlZmZlZ6DlZmZlZ6DlZmZlZ6DlZmZlZ6DlZm\nZlZ6DlZmZlZ6DlZmZlZ6DlZmZlZ6DlZmZlZ6DlZmZlZ6DlZmZlZ6TQtWkhZKWirpIUmX1UifJ+lW\nSfdIul/SGxud18zMxhdFxIFfqNQGPAicDTwB/Ay4MCIeKOTpAqZExOWS5uX884EYat48fzSj7GZm\nY5kkIkKZkfh5AAAK9klEQVSjXY7halbPagGwLCKWR8Qu4Ebg/Ko8q4BZeXgWsDYidjc4r5mZjSPN\nClZHACsK44/naUVfBE6UtBJYArxrGPOamdk40qxg1cj1ufcB90TE4cDzgc9Jmtmk8piZWQub2KTl\nPgEcWRg/ktRDKvpfwN8BRMSvJT0CHJfzDTUvAF1dXX3DnZ2ddHZ27mexzczqi4Ddu2HnzoGvHTv6\nh2fPho6O0S5pv+7ubrq7u0e7GPutWTdYTCTdJHEWsBK4i71vsPgUsDEiPihpPvBz4GRg01Dz5vl9\ng4XZGNPb23/gr/c+3LTq4Vrj9abVmj5xIkyeXPs1ZQq8+tXwgQ+Mdk3W16o3WDSlZxURuyVdCtwG\ntAFXR8QDki7J6VcBfw9cK2kJ6XLkeyNiHUCteZtRTrPxqjooDPUaKl8xvTpgNDJced+zJx3wKwf+\n4bzXm3bQQXtPq85fa9qkSXsve9IkmOBfp46KpvSsRoJ7VtZKImDXrnRA7unpP1g3MlxrvN60waYX\nX7t39x+AG3kNlbeY3uhwrbSJE0Etd87fWlq1Z+VgZWNeb286c+/pqf2qHOBrjQ9nuPq9elpbWzog\nT53af5BuZLjWeL1pg00vviZNclAYrxysRpiDVeuofCnd0wPbt/cfyIvDjYzv62vXroEH8urXlCkw\nbdrAA30xbV+m1XpvaxvtT8LMwWrEOVjtmz17UhCovCpBYbjDw32XUkCoHNzrDRfHK0GkXr7BAk9x\n2uTJ7kWYVThYjbCxGKw2bICVKwcGk8pr27ba06vTioGl1mv37nQwr7wqQWC447Xe66VVLjuZ2ehr\n1WDVrN9Z2T742tfgk58cGByKr+nTB47PnVs/b72Xexlm1orcszIzG0datWflXwyYmVnpOViZmVnp\nOViZmVnpOViZmVnpOViZmVnpOViZmVnpOViZmVnpOViZmVnpOViZmVnpOViZmVnpOViZmVnpOViZ\nmVnpOViZmVnpNTVYSVooaamkhyRdViP9PZLuzq/7JO2WNDunLZd0b067q5nlNDOzcmtasJLUBlwJ\nLAROAC6UdHwxT0R8MiJOjYhTgcuB7ojYUEkGOnP6gmaVcyzo7u4e7SKUhuuin+uin+ui9TWzZ7UA\nWBYRyyNiF3AjcP4g+V8H3FA1reX+c2U0eEfs57ro57ro57pofc0MVkcAKwrjj+dpe5E0HXgp8LXC\n5AC+K2mxpIubVkozMyu9Zv6t/XD+xvcVwA8LlwABXhQRqyQdAnxH0tKIuOPAFtHMzFpB0/7WXtLp\nQFdELMzjlwO9EfGxGnlvAv41Im6ss6xFwJaIuKIwzf9pb2a2D1rxb+2bGawmAg8CZwErgbuACyPi\ngap87cDDwDMjYnueNh1oi4jNkmYA3wY+GBHfbkphzcys1Jp2GTAidku6FLgNaAOujogHJF2S06/K\nWX8PuK0SqLL5wE2SKmW83oHKzGz8alrPyszM7EAp5RMs/GPifg3UxTxJt0q6R9L9kt7Y6LytZj/r\nYry1izmSbpK0RNKdkk5sdN5Ws591MWbahaRrJK2WdN8geT6b62mJpFML08vfJiKiVC/SJcNlQAcw\nCbgHOH6Q/OcB3y2MPwLMHe3tGKm6ALqAj+ThecBa0qXTYdVj2V/7UxfjtF18AvhAHj6uso+M03ZR\nsy7GYLt4MXAqcF+d9HOBb+bhFwI/baU2UcaelX9M3K+RulgFzMrDs4C1EbG7wXlbyf7URcV4ahfH\nA7cDRMSDQIekQxuct5Xsa10cUkgfE+0i0k971g+S5ZXAl3LeO4HZkp5Bi7SJMgYr/5i4XyN18UXg\nREkrgSXAu4YxbyvZn7qA8dculgCvApC0ADgaeGaD87aS/akLGFvtYij16urwOtNLpZk/Ct5X/jFx\nv0bq4n3APRHRKek3SNt8SpPLNRr2uS4iYjPjr118FPiMpLuB+4C7gT0NzttK9qcuAM6IiJVjpF00\nomV7kWXsWT0BHFkYP5IU6Wv5Q6ouAUbEqvy+BriJ1MVtVY3Uxf8C/h0gIn5NugZ/XM7XaD22gv2p\ni3HXLiJic0S8OdKDoC8CDgF+3ci8LWZf6+LhnLYyv4+FdjGU6rp6JqmuWqJNlDFYLQaOldQhaTJw\nAXBLdab8Y+LfAW4uTJsuaWYengGcQzqTalWN1MVS4GwASfNJB+eHG5y3lexzXYzHdiGpPaeRL299\nPyK2NDJvi9nnuhiD7WIotwAXQd8ThjZExGpapE2U7jJg+MfEfRqsi78HrpW0hHTy8d6IWAdQa97R\n2I4DYX/qQtKzga+Ps3ZxAnCd0mPJ7gfeMti8o7EdB8L+1AVj7Hgh6QbgJcA8SSuARaS7+4iIqyLi\nm5LOlbQM2Aq8Kae1RJvwj4LNzKz0yngZ0MzMbAAHKzMzKz0HKzMzKz0HKzMzKz0HKzMzKz0HKzMz\nKz0HK2tZkvao/29i/k3StNEu02AkvVHSP+bhSyS9vjD9sKq8/5ofGTXY8rolvSAPv1bSryR9b5D8\n8yV9c/+3xGzkOVhZK9uWH6FzErATeFszVyZpf39E3/ejxvwjza/k0TeQHiZaWc8xwIz8yKihlldZ\n5luAP4mIs+pmTk8rWC/pN/el8GajycHKxoofAsco/dHef+Y/l/uJpJMA8h/szVKyttCr+bKksyRN\nkPQJSXfled+a0zsl3SHpZuCXxRVKapN0Xe7Z3SvpXXl6t6RPF3p9v1WZpTBvl6R3S3o1cBpwvaRf\nSJpKeublLYW8/yTpZ0p/KNlVtd2S9DfAi4BrJH1c0tGSfiDp5/n124X8twAX7mddm404BytrebnH\nsxC4F/gQ8POIOIX0FPYv52w/As4ATiQ90PWMPP104MfAn5CelbaA9DDTiyV15DynAu+MiOOqVv18\n4PCIOCkiTgauzdMDmBYRpwLvAK6pUewAIiK+Rno22+si4jcjoocUeBYX8r4/In4LOAV4SSUAV5YT\nER8qLOO9wFPA/46IF5AC32cL+e8iPVPTrKWU7tmAZsMwTelvHwB+QAoKd5L/uygibpd0cH5Y6R2k\ng/SjwOeBt0o6HFgfEdslnQOcJOk1eXmzgGOA3cBdEfFojfX/Gni2pM8C/wUUnyt3Qy7DHblH1z7E\nthT/uuFo0h9JVlyQH8A6ETiM9GeCtR64WlnGZOBKpb+K2QM8p5BnFekfYc1aioOVtbLtuffSJz+U\ntPo/e4IUzC4FlgPvB34feE2eXnFpRHynanmdpId+7iUiNkg6mdSrexvwB/Q/JLVa7xDbUv2QTuX1\nPwt4N3BaRGyUdC0wdYhl/DmwKiJeL6kN6Klarh8Iai3HlwFtrLkD+CPoCzRrImJLRDwOzAOOiYhH\nSN9xvYf+YHUb8I7KTRSSnqP0T9R1SToYmBgRXwc+QLpcCCkgXJDznEG6vLi5enb6g+pmUk+u4lFS\nD4o8fSuwSelvT17WQB3MAp7MwxeRnqRdcVhevllLcc/KWlmtHkIX6UaDJaSD/BsKaT+l/wTth6S/\nFPlhHv+/pMtjv1Dqnj1F6n0V77irdgTpL0kqy/yrQrl6JP2CtI+9uTA9agxfB3xB0jbSH0j+kHTT\nxc8jYkm+1LmU9NfjlfIO5p+Ar0m6CLgV2FJIW8DA3qRZS/BfhJgdYJJuB94dEb/Yx/mfDfxjRLz8\nwJYMJF0PfDIi7h4ys1mJ+DKgWclExMPA5qF+FDxckg4FZjtQWStyz8rMzErPPSszMys9ByszMys9\nByszMys9ByszMys9ByszMys9ByszMyu9/w9P6hXBX++EYwAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x5b24e50>"
       ]
      }
     ],
     "prompt_number": 3
    }
   ],
   "metadata": {}
  }
 ]
}