chap15_solved.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "romantic-speech",
      "metadata": {
        "id": "romantic-speech"
      },
      "source": [
        "# Cooling Coffee"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "higher-button",
      "metadata": {
        "tags": [],
        "id": "higher-button"
      },
      "source": [
        "*Modeling and Simulation in Python*\n",
        "\n",
        "Copyright 2021 Allen Downey\n",
        "\n",
        "License: [Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International](https://creativecommons.org/licenses/by-nc-sa/4.0/)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "id": "dependent-monitoring",
      "metadata": {
        "tags": [],
        "id": "dependent-monitoring",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "390c304f-b9ed-442f-c1ef-5021b1479fba"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Downloaded modsim.py\n"
          ]
        }
      ],
      "source": [
        "# download modsim.py if necessary\n",
        "\n",
        "from os.path import basename, exists\n",
        "\n",
        "def download(url):\n",
        "    filename = basename(url)\n",
        "    if not exists(filename):\n",
        "        from urllib.request import urlretrieve\n",
        "        local, _ = urlretrieve(url, filename)\n",
        "        print('Downloaded ' + local)\n",
        "\n",
        "download('https://github.com/AllenDowney/ModSimPy/raw/master/' +\n",
        "         'modsim.py')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "id": "continental-decrease",
      "metadata": {
        "tags": [],
        "id": "continental-decrease"
      },
      "outputs": [],
      "source": [
        "# import functions from modsim\n",
        "\n",
        "from modsim import *"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "moving-trademark",
      "metadata": {
        "tags": [],
        "id": "moving-trademark"
      },
      "source": [
        "This chapter is available as a Jupyter notebook where you can read the text, run the code, and work on the exercises.\n",
        "Click here to access the notebooks: <https://allendowney.github.io/ModSimPy/>."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "capital-needle",
      "metadata": {
        "id": "capital-needle"
      },
      "source": [
        "So far the systems we have studied have been physical in the sense that they exist in the world, but they have not been physics in the sense of what physics classes are usually about. In the next few chapters, we'll do some physics, starting with *thermal systems*, that is, systems where the temperature of objects changes as heat transfers from one to another."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "working-senator",
      "metadata": {
        "id": "working-senator"
      },
      "source": [
        "## The Coffee Cooling Problem\n",
        "\n",
        "The coffee cooling problem was discussed by Jearl Walker in\n",
        "\"The Amateur Scientist\", *Scientific American*, Volume 237, Issue 5, November 1977. Since then it has become a standard example of modeling and simulation.\n",
        "\n",
        "Here is my version of the problem:\n",
        "\n",
        "> Suppose I stop on the way to work to pick up a cup of coffee and a small container of milk. Assuming that I want the coffee to be as hot as possible when I arrive at work, should I add the milk at the coffee shop, wait until I get to work, or add the milk at some point in between?"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "backed-steel",
      "metadata": {
        "id": "backed-steel"
      },
      "source": [
        "To help answer this question, I made a trial run with the milk and\n",
        "coffee in separate containers and took some measurements (not really):\n",
        "\n",
        "-   When served, the temperature of the coffee is 90 °C. The volume is\n",
        "    300 mL.\n",
        "\n",
        "-   The milk is at an initial temperature of 5 °C, and I take about\n",
        "    50 mL.\n",
        "\n",
        "-   The ambient temperature in my car is 22 °C.\n",
        "\n",
        "-   The coffee is served in a well insulated cup. When I arrive at work after 30 minutes, the temperature of the coffee has fallen to 70 °C.\n",
        "\n",
        "-   The milk container is not as well insulated. After 15 minutes, it\n",
        "    warms up to 20 °C, nearly the ambient temperature.\n",
        "\n",
        "To use this data and answer the question, we have to know something\n",
        "about temperature and heat, and we have to make some modeling decisions."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "retired-noise",
      "metadata": {
        "id": "retired-noise"
      },
      "source": [
        "## Temperature and Heat\n",
        "\n",
        "To understand how coffee cools (and milk warms), we need a model of\n",
        "temperature and heat. *Temperature* is a property of an object or a\n",
        "system; in SI units it is measured in degrees Celsius (°C). Temperature quantifies how hot or cold the object is, which is related to the average velocity of the particles that make it up.\n",
        "\n",
        "When particles in a hot object contact particles in a cold object, the\n",
        "hot object gets cooler and the cold object gets warmer as energy is\n",
        "transferred from one to the other. The transferred energy is called\n",
        "*heat*; in SI units it is measured in joules (J).\n",
        "\n",
        "Heat is related to temperature by the following equation (see\n",
        "<https://modsimpy.com/thermass>):\n",
        "\n",
        "$$Q = C~\\Delta T$$\n",
        "\n",
        "where $Q$ is the amount of heat transferred to an object, $\\Delta T$ is the resulting change in temperature, and $C$ is the object's *thermal mass*, which is a property of the object that determines how much energy it takes to heat or cool it. In SI units, thermal mass is measured in joules per degree Celsius (J/°C)."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "upper-johnson",
      "metadata": {
        "id": "upper-johnson"
      },
      "source": [
        "For objects made primarily from one material, thermal mass can be\n",
        "computed like this:\n",
        "\n",
        "$$C = m c_p$$\n",
        "\n",
        "where $m$ is the mass of the object and $c_p$ is the *specific heat capacity* of the material, which is the amount of thermal mass per gram (see <https://modsimpy.com/specheat>).\n",
        "\n",
        "We can use these equations to estimate the thermal mass of a cup of\n",
        "coffee. The specific heat capacity of coffee is probably close to that\n",
        "of water, which is 4.2 J/g/°C. Assuming that the density of coffee is\n",
        "close to that of water, which is 1 g/mL, the mass of 300 mL of coffee is 300 g, and the thermal mass is 1260 J/°C.\n",
        "\n",
        "So when a cup of coffee cools from 90 °C to 70 °C, the change in\n",
        "temperature, $\\Delta T$ is 20 °C, which means that 25 200 J of heat\n",
        "energy was transferred from the cup and the coffee to the surrounding environment\n",
        "(the cup holder and air in my car).\n",
        "\n",
        "To give you a sense of how much energy that is, if you were able to\n",
        "harness all of that heat to do work (which you cannot), you could\n",
        "use it to lift a cup of coffee from sea level to 8571 m, just shy of the height of Mount Everest, 8848 m."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "closed-vegetable",
      "metadata": {
        "id": "closed-vegetable"
      },
      "source": [
        "## Heat Transfer\n",
        "\n",
        "In a situation like the coffee cooling problem, there are three ways\n",
        "heat transfers from one object to another (see <https://modsimpy.com/transfer>):\n",
        "\n",
        "-   Conduction: When objects at different temperatures come into\n",
        "    contact, the faster-moving particles of the higher-temperature\n",
        "    object transfer kinetic energy to the slower-moving particles of the lower-temperature object.\n",
        "\n",
        "-   Convection: When particles in a gas or liquid flow from place to\n",
        "    place, they carry heat energy with them. Fluid flows can be caused\n",
        "    by external action, like stirring, or by internal differences in\n",
        "    temperature. For example, you might have heard that hot air rises,\n",
        "    which is a form of \"natural convection\".\n",
        "\n",
        "-   Radiation: As the particles in an object move due to thermal energy,\n",
        "    they emit electromagnetic radiation. The energy carried by this\n",
        "    radiation depends on the object's temperature and surface properties\n",
        "    (see <https://modsimpy.com/thermrad>)."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "valuable-shark",
      "metadata": {
        "id": "valuable-shark"
      },
      "source": [
        "For objects like coffee in a car, the effect of radiation is much\n",
        "smaller than the effects of conduction and convection, so we will ignore it.\n",
        "\n",
        "Convection can be a complex topic, since it often depends on details of fluid flow in three dimensions. But for this problem we will be able to get away with a simple model called \"Newton's law of cooling\"."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "freelance-target",
      "metadata": {
        "id": "freelance-target"
      },
      "source": [
        "## Newton's Law of Cooling\n",
        "\n",
        "*Newton's law of cooling* asserts that the temperature rate of change for an object is proportional to the difference in temperature between the object and the surrounding environment:\n",
        "\n",
        "$$\\frac{dT}{dt} = -r (T - T_{env})$$\n",
        "\n",
        "where $t$ is time, $T$ is the temperature of the object, $T_{env}$ is the temperature of the environment, and $r$ is a constant that characterizes how quickly heat is transferred between the object and the environment."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "corrected-attachment",
      "metadata": {
        "id": "corrected-attachment"
      },
      "source": [
        "Newton's so-called \"law \" is really a model: it is a good approximation in some conditions and less good in others.\n",
        "\n",
        "For example, if the primary mechanism of heat transfer is conduction,\n",
        "Newton's law is \"true\", which is to say that $r$ is constant over a\n",
        "wide range of temperatures. And sometimes we can estimate $r$ based on\n",
        "the material properties and shape of the object.\n",
        "\n",
        "When convection contributes a non-negligible fraction of heat transfer, $r$ depends on temperature, but Newton's law is often accurate enough, at least over a narrow range of temperatures. In this case $r$ usually has to be estimated experimentally, since it depends on details of surface shape, air flow, evaporation, etc.\n",
        "\n",
        "When radiation makes up a substantial part of heat transfer, Newton's\n",
        "law is not a good model at all. This is the case for objects in space or in a vacuum, and for objects at high temperatures (more than a few\n",
        "hundred degrees Celsius, say).\n",
        "\n",
        "However, for a situation like the coffee cooling problem, we expect Newton's model to be quite good.\n",
        "\n",
        "With that, we have just one more modeling decision to make: whether to treat the coffee and the cup as separate objects or a single object. If the cup is made of paper, it has less mass than the coffee, and the specific heat capacity of paper is lower, too. In that case, it would be reasonable to treat the cup and coffee as a single object. For a cup with substantial thermal mass, like a ceramic mug, we might consider a model that computes the temperature of coffee and cup separately."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "selective-video",
      "metadata": {
        "id": "selective-video"
      },
      "source": [
        "## Implementing Newtonian Cooling\n",
        "\n",
        "To get started, we'll focus on the coffee. Then, as an exercise, you can simulate the milk. In the next chapter, we'll put them together, literally.\n",
        "\n",
        "Here's a function that takes the parameters of the system and makes a `System` object:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "id": "incredible-greeting",
      "metadata": {
        "tags": [],
        "id": "incredible-greeting"
      },
      "outputs": [],
      "source": [
        "def make_system(T_init, volume, r, t_end):\n",
        "    return System(T_init=T_init,\n",
        "                  T_final=T_init,\n",
        "                  volume=volume,\n",
        "                  r=r,\n",
        "                  t_end=t_end,\n",
        "                  T_env=22,\n",
        "                  t_0=0,\n",
        "                  dt=1)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "gross-appeal",
      "metadata": {
        "id": "gross-appeal"
      },
      "source": [
        "In addition to the parameters, `make_system` sets the temperature of the environment, `T_env`, the initial time stamp, `t_0`, and the time step, `dt`, which we will use to simulate the cooling process.\n",
        "Here's a `System` object that represents the coffee."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "id": "fuzzy-support",
      "metadata": {
        "id": "fuzzy-support"
      },
      "outputs": [],
      "source": [
        "coffee = make_system(T_init=90, volume=300, r=0.01, t_end=30)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "rental-pearl",
      "metadata": {
        "id": "rental-pearl"
      },
      "source": [
        "The values of `T_init`, `volume`, and `t_end` come from the statement of the problem.\n",
        "I chose the value of `r` arbitrarily for now; we will see how to estimate it soon.\n",
        "\n",
        "Strictly speaking, Newton's law is a differential equation, but over a short period of time we can approximate it with a difference equation:\n",
        "\n",
        "$$\\Delta T = -r (T - T_{env}) dt$$\n",
        "\n",
        "where $dt$ is the time step and $\\Delta T$ is the change in temperature during that time step.\n",
        "\n",
        "Note: I use $\\Delta T$ to denote a change in temperature over time, but in the context of heat transfer, you might also see $\\Delta T$ used to denote the difference in temperature between an object and its\n",
        "environment, $T - T_{env}$. To minimize confusion, I avoid this second\n",
        "use.\n",
        "\n",
        "The following function takes the current time `t`, the current temperature, `T`, and a `System` object, and computes the change in temperature during a time step:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "id": "tight-baptist",
      "metadata": {
        "tags": [],
        "id": "tight-baptist"
      },
      "outputs": [],
      "source": [
        "def change_func(t, T, system):\n",
        "    r, T_env, dt = system.r, system.T_env, system.dt\n",
        "    return -r * (T - T_env) * dt"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "flexible-crest",
      "metadata": {
        "id": "flexible-crest"
      },
      "source": [
        "We can test it with the initial temperature of the coffee, like this:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "id": "fiscal-artwork",
      "metadata": {
        "id": "fiscal-artwork",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "a7524022-d65e-4ce6-f302-b3eeb95a700d"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "-0.68"
            ]
          },
          "metadata": {},
          "execution_count": 6
        }
      ],
      "source": [
        "change_func(0, coffee.T_init, coffee)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "sporting-morocco",
      "metadata": {
        "id": "sporting-morocco"
      },
      "source": [
        "With `dt=1` minute, the temperature drops by about 0.7 °C, at least for this value of `r`.\n",
        "\n",
        "Now here's a version of `run_simulation` that simulates a series of time steps from `t_0` to `t_end`:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "id": "cardiac-independence",
      "metadata": {
        "tags": [],
        "id": "cardiac-independence"
      },
      "outputs": [],
      "source": [
        "def run_simulation(system, change_func):\n",
        "    t_array = linrange(system.t_0, system.t_end, system.dt)\n",
        "    n = len(t_array)\n",
        "\n",
        "    series = TimeSeries(index=t_array)\n",
        "    series.iloc[0] = system.T_init\n",
        "\n",
        "    for i in range(n-1):\n",
        "        t = t_array[i]\n",
        "        T = series.iloc[i]\n",
        "        series.iloc[i+1] = T + change_func(t, T, system)\n",
        "\n",
        "    system.T_final = series.iloc[-1]\n",
        "    return series"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "together-adapter",
      "metadata": {
        "id": "together-adapter"
      },
      "source": [
        "There are two things here that are different from previous versions of `run_simulation`.\n",
        "\n",
        "First, we use  `linrange` to make an array of values from `t_0` to `t_end` with time step `dt`.\n",
        "`linrange` is similar to `linspace`; they both take a start value and an end value and return an array of equally spaced values.\n",
        "The difference is the third argument: `linspace` takes an integer that indicates the number of points in the range; `linrange` takes a step size that indicates the interval between values.\n",
        "When we make the `TimeSeries`, we use the keyword argument `index` to indicate that the index of the `TimeSeries` is the array of time stamps, `t_array`."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "fantastic-object",
      "metadata": {
        "id": "fantastic-object"
      },
      "source": [
        "Second, this version of `run_simulation` uses `iloc` rather than `loc` to specify the rows in the `TimeSeries`.\n",
        "Here's the difference:\n",
        "\n",
        "* With `loc`, the label in brackets can be any kind of value, with any start, end, and time step.  For example, in the world population model, the labels are years starting in 1960 and ending in 2016.\n",
        "\n",
        "* With `iloc`, the label in brackets is always an integer starting at 0.  So we can always get the first element with `iloc[0]` and the last element with `iloc[-1]`, regardless of what the labels are.\n",
        "\n",
        "In this version of `run_simulation`, the loop variable is an integer, `i`, that goes from `0` to `n-1`, including `0` but not including `n-1`.\n",
        "So the first time through the loop, `i` is `0` and the value we add to the `TimeSeries` has index 1.\n",
        "The last time through the loop, `i` is `n-2` and the value we add has index `n-1`.\n",
        "\n",
        "We can run the simulation like this:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "id": "numerous-metabolism",
      "metadata": {
        "id": "numerous-metabolism"
      },
      "outputs": [],
      "source": [
        "results = run_simulation(coffee, change_func)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "earned-primary",
      "metadata": {
        "id": "earned-primary"
      },
      "source": [
        "The result is a `TimeSeries` with one row per time step.\n",
        "Here are the first few rows:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "id": "infectious-carolina",
      "metadata": {
        "id": "infectious-carolina",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 237
        },
        "outputId": "ae21100a-9276-4935-d094-888ceec439bf"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "       Quantity\n",
              "Time           \n",
              "0.0   90.000000\n",
              "1.0   89.320000\n",
              "2.0   88.646800\n",
              "3.0   87.980332\n",
              "4.0   87.320529"
            ],
            "text/html": [
              "\n",
              "  <div id=\"df-604f6f63-b228-4928-9926-894c5918ada8\" class=\"colab-df-container\">\n",
              "    <div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>Quantity</th>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>Time</th>\n",
              "      <th></th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0.0</th>\n",
              "      <td>90.000000</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1.0</th>\n",
              "      <td>89.320000</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2.0</th>\n",
              "      <td>88.646800</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3.0</th>\n",
              "      <td>87.980332</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>4.0</th>\n",
              "      <td>87.320529</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>\n",
              "    <div class=\"colab-df-buttons\">\n",
              "\n",
              "  <div class=\"colab-df-container\">\n",
              "    <button class=\"colab-df-convert\" onclick=\"convertToInteractive('df-604f6f63-b228-4928-9926-894c5918ada8')\"\n",
              "            title=\"Convert this dataframe to an interactive table.\"\n",
              "            style=\"display:none;\">\n",
              "\n",
              "  <svg xmlns=\"https://www.w3.org/2000/svg\" height=\"24px\" viewBox=\"0 -960 960 960\">\n",
              "    <path d=\"M120-120v-720h720v720H120Zm60-500h600v-160H180v160Zm220 220h160v-160H400v160Zm0 220h160v-160H400v160ZM180-400h160v-160H180v160Zm440 0h160v-160H620v160ZM180-180h160v-160H180v160Zm440 0h160v-160H620v160Z\"/>\n",
              "  </svg>\n",
              "    </button>\n",
              "\n",
              "  <style>\n",
              "    .colab-df-container {\n",
              "      display:flex;\n",
              "      gap: 12px;\n",
              "    }\n",
              "\n",
              "    .colab-df-convert {\n",
              "      background-color: #E8F0FE;\n",
              "      border: none;\n",
              "      border-radius: 50%;\n",
              "      cursor: pointer;\n",
              "      display: none;\n",
              "      fill: #1967D2;\n",
              "      height: 32px;\n",
              "      padding: 0 0 0 0;\n",
              "      width: 32px;\n",
              "    }\n",
              "\n",
              "    .colab-df-convert:hover {\n",
              "      background-color: #E2EBFA;\n",
              "      box-shadow: 0px 1px 2px rgba(60, 64, 67, 0.3), 0px 1px 3px 1px rgba(60, 64, 67, 0.15);\n",
              "      fill: #174EA6;\n",
              "    }\n",
              "\n",
              "    .colab-df-buttons div {\n",
              "      margin-bottom: 4px;\n",
              "    }\n",
              "\n",
              "    [theme=dark] .colab-df-convert {\n",
              "      background-color: #3B4455;\n",
              "      fill: #D2E3FC;\n",
              "    }\n",
              "\n",
              "    [theme=dark] .colab-df-convert:hover {\n",
              "      background-color: #434B5C;\n",
              "      box-shadow: 0px 1px 3px 1px rgba(0, 0, 0, 0.15);\n",
              "      filter: drop-shadow(0px 1px 2px rgba(0, 0, 0, 0.3));\n",
              "      fill: #FFFFFF;\n",
              "    }\n",
              "  </style>\n",
              "\n",
              "    <script>\n",
              "      const buttonEl =\n",
              "        document.querySelector('#df-604f6f63-b228-4928-9926-894c5918ada8 button.colab-df-convert');\n",
              "      buttonEl.style.display =\n",
              "        google.colab.kernel.accessAllowed ? 'block' : 'none';\n",
              "\n",
              "      async function convertToInteractive(key) {\n",
              "        const element = document.querySelector('#df-604f6f63-b228-4928-9926-894c5918ada8');\n",
              "        const dataTable =\n",
              "          await google.colab.kernel.invokeFunction('convertToInteractive',\n",
              "                                                    [key], {});\n",
              "        if (!dataTable) return;\n",
              "\n",
              "        const docLinkHtml = 'Like what you see? Visit the ' +\n",
              "          '<a target=\"_blank\" href=https://colab.research.google.com/notebooks/data_table.ipynb>data table notebook</a>'\n",
              "          + ' to learn more about interactive tables.';\n",
              "        element.innerHTML = '';\n",
              "        dataTable['output_type'] = 'display_data';\n",
              "        await google.colab.output.renderOutput(dataTable, element);\n",
              "        const docLink = document.createElement('div');\n",
              "        docLink.innerHTML = docLinkHtml;\n",
              "        element.appendChild(docLink);\n",
              "      }\n",
              "    </script>\n",
              "  </div>\n",
              "\n",
              "\n",
              "<div id=\"df-38cadfa8-edae-45c8-9563-4218da13c85b\">\n",
              "  <button class=\"colab-df-quickchart\" onclick=\"quickchart('df-38cadfa8-edae-45c8-9563-4218da13c85b')\"\n",
              "            title=\"Suggest charts\"\n",
              "            style=\"display:none;\">\n",
              "\n",
              "<svg xmlns=\"https://www.w3.org/2000/svg\" height=\"24px\"viewBox=\"0 0 24 24\"\n",
              "     width=\"24px\">\n",
              "    <g>\n",
              "        <path d=\"M19 3H5c-1.1 0-2 .9-2 2v14c0 1.1.9 2 2 2h14c1.1 0 2-.9 2-2V5c0-1.1-.9-2-2-2zM9 17H7v-7h2v7zm4 0h-2V7h2v10zm4 0h-2v-4h2v4z\"/>\n",
              "    </g>\n",
              "</svg>\n",
              "  </button>\n",
              "\n",
              "<style>\n",
              "  .colab-df-quickchart {\n",
              "      --bg-color: #E8F0FE;\n",
              "      --fill-color: #1967D2;\n",
              "      --hover-bg-color: #E2EBFA;\n",
              "      --hover-fill-color: #174EA6;\n",
              "      --disabled-fill-color: #AAA;\n",
              "      --disabled-bg-color: #DDD;\n",
              "  }\n",
              "\n",
              "  [theme=dark] .colab-df-quickchart {\n",
              "      --bg-color: #3B4455;\n",
              "      --fill-color: #D2E3FC;\n",
              "      --hover-bg-color: #434B5C;\n",
              "      --hover-fill-color: #FFFFFF;\n",
              "      --disabled-bg-color: #3B4455;\n",
              "      --disabled-fill-color: #666;\n",
              "  }\n",
              "\n",
              "  .colab-df-quickchart {\n",
              "    background-color: var(--bg-color);\n",
              "    border: none;\n",
              "    border-radius: 50%;\n",
              "    cursor: pointer;\n",
              "    display: none;\n",
              "    fill: var(--fill-color);\n",
              "    height: 32px;\n",
              "    padding: 0;\n",
              "    width: 32px;\n",
              "  }\n",
              "\n",
              "  .colab-df-quickchart:hover {\n",
              "    background-color: var(--hover-bg-color);\n",
              "    box-shadow: 0 1px 2px rgba(60, 64, 67, 0.3), 0 1px 3px 1px rgba(60, 64, 67, 0.15);\n",
              "    fill: var(--button-hover-fill-color);\n",
              "  }\n",
              "\n",
              "  .colab-df-quickchart-complete:disabled,\n",
              "  .colab-df-quickchart-complete:disabled:hover {\n",
              "    background-color: var(--disabled-bg-color);\n",
              "    fill: var(--disabled-fill-color);\n",
              "    box-shadow: none;\n",
              "  }\n",
              "\n",
              "  .colab-df-spinner {\n",
              "    border: 2px solid var(--fill-color);\n",
              "    border-color: transparent;\n",
              "    border-bottom-color: var(--fill-color);\n",
              "    animation:\n",
              "      spin 1s steps(1) infinite;\n",
              "  }\n",
              "\n",
              "  @keyframes spin {\n",
              "    0% {\n",
              "      border-color: transparent;\n",
              "      border-bottom-color: var(--fill-color);\n",
              "      border-left-color: var(--fill-color);\n",
              "    }\n",
              "    20% {\n",
              "      border-color: transparent;\n",
              "      border-left-color: var(--fill-color);\n",
              "      border-top-color: var(--fill-color);\n",
              "    }\n",
              "    30% {\n",
              "      border-color: transparent;\n",
              "      border-left-color: var(--fill-color);\n",
              "      border-top-color: var(--fill-color);\n",
              "      border-right-color: var(--fill-color);\n",
              "    }\n",
              "    40% {\n",
              "      border-color: transparent;\n",
              "      border-right-color: var(--fill-color);\n",
              "      border-top-color: var(--fill-color);\n",
              "    }\n",
              "    60% {\n",
              "      border-color: transparent;\n",
              "      border-right-color: var(--fill-color);\n",
              "    }\n",
              "    80% {\n",
              "      border-color: transparent;\n",
              "      border-right-color: var(--fill-color);\n",
              "      border-bottom-color: var(--fill-color);\n",
              "    }\n",
              "    90% {\n",
              "      border-color: transparent;\n",
              "      border-bottom-color: var(--fill-color);\n",
              "    }\n",
              "  }\n",
              "</style>\n",
              "\n",
              "  <script>\n",
              "    async function quickchart(key) {\n",
              "      const quickchartButtonEl =\n",
              "        document.querySelector('#' + key + ' button');\n",
              "      quickchartButtonEl.disabled = true;  // To prevent multiple clicks.\n",
              "      quickchartButtonEl.classList.add('colab-df-spinner');\n",
              "      try {\n",
              "        const charts = await google.colab.kernel.invokeFunction(\n",
              "            'suggestCharts', [key], {});\n",
              "      } catch (error) {\n",
              "        console.error('Error during call to suggestCharts:', error);\n",
              "      }\n",
              "      quickchartButtonEl.classList.remove('colab-df-spinner');\n",
              "      quickchartButtonEl.classList.add('colab-df-quickchart-complete');\n",
              "    }\n",
              "    (() => {\n",
              "      let quickchartButtonEl =\n",
              "        document.querySelector('#df-38cadfa8-edae-45c8-9563-4218da13c85b button');\n",
              "      quickchartButtonEl.style.display =\n",
              "        google.colab.kernel.accessAllowed ? 'block' : 'none';\n",
              "    })();\n",
              "  </script>\n",
              "</div>\n",
              "\n",
              "    </div>\n",
              "  </div>\n"
            ],
            "application/vnd.google.colaboratory.intrinsic+json": {
              "type": "dataframe",
              "summary": "{\n  \"name\": \"show(results\",\n  \"rows\": 5,\n  \"fields\": [\n    {\n      \"column\": \"Time\",\n      \"properties\": {\n        \"dtype\": \"number\",\n        \"std\": 1.5811388300841898,\n        \"min\": 0.0,\n        \"max\": 4.0,\n        \"num_unique_values\": 5,\n        \"samples\": [\n          1.0,\n          4.0,\n          2.0\n        ],\n        \"semantic_type\": \"\",\n        \"description\": \"\"\n      }\n    },\n    {\n      \"column\": \"Quantity\",\n      \"properties\": {\n        \"dtype\": \"number\",\n        \"std\": 1.0591620604283003,\n        \"min\": 87.32052868000001,\n        \"max\": 90.0,\n        \"num_unique_values\": 5,\n        \"samples\": [\n          89.32,\n          87.32052868000001,\n          88.6468\n        ],\n        \"semantic_type\": \"\",\n        \"description\": \"\"\n      }\n    }\n  ]\n}"
            }
          },
          "metadata": {},
          "execution_count": 9
        }
      ],
      "source": [
        "show(results.head())"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "color-olive",
      "metadata": {
        "id": "color-olive"
      },
      "source": [
        "And the last few rows:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "id": "accompanied-melbourne",
      "metadata": {
        "id": "accompanied-melbourne",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 237
        },
        "outputId": "f9582bcc-5554-46cb-f42e-c80c782a255a"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "       Quantity\n",
              "Time           \n",
              "26.0  74.362934\n",
              "27.0  73.839305\n",
              "28.0  73.320912\n",
              "29.0  72.807702\n",
              "30.0  72.299625"
            ],
            "text/html": [
              "\n",
              "  <div id=\"df-1a5005cb-18cd-4b07-85f2-e3680f0e9d28\" class=\"colab-df-container\">\n",
              "    <div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>Quantity</th>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>Time</th>\n",
              "      <th></th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>26.0</th>\n",
              "      <td>74.362934</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>27.0</th>\n",
              "      <td>73.839305</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>28.0</th>\n",
              "      <td>73.320912</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>29.0</th>\n",
              "      <td>72.807702</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>30.0</th>\n",
              "      <td>72.299625</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>\n",
              "    <div class=\"colab-df-buttons\">\n",
              "\n",
              "  <div class=\"colab-df-container\">\n",
              "    <button class=\"colab-df-convert\" onclick=\"convertToInteractive('df-1a5005cb-18cd-4b07-85f2-e3680f0e9d28')\"\n",
              "            title=\"Convert this dataframe to an interactive table.\"\n",
              "            style=\"display:none;\">\n",
              "\n",
              "  <svg xmlns=\"https://www.w3.org/2000/svg\" height=\"24px\" viewBox=\"0 -960 960 960\">\n",
              "    <path d=\"M120-120v-720h720v720H120Zm60-500h600v-160H180v160Zm220 220h160v-160H400v160Zm0 220h160v-160H400v160ZM180-400h160v-160H180v160Zm440 0h160v-160H620v160ZM180-180h160v-160H180v160Zm440 0h160v-160H620v160Z\"/>\n",
              "  </svg>\n",
              "    </button>\n",
              "\n",
              "  <style>\n",
              "    .colab-df-container {\n",
              "      display:flex;\n",
              "      gap: 12px;\n",
              "    }\n",
              "\n",
              "    .colab-df-convert {\n",
              "      background-color: #E8F0FE;\n",
              "      border: none;\n",
              "      border-radius: 50%;\n",
              "      cursor: pointer;\n",
              "      display: none;\n",
              "      fill: #1967D2;\n",
              "      height: 32px;\n",
              "      padding: 0 0 0 0;\n",
              "      width: 32px;\n",
              "    }\n",
              "\n",
              "    .colab-df-convert:hover {\n",
              "      background-color: #E2EBFA;\n",
              "      box-shadow: 0px 1px 2px rgba(60, 64, 67, 0.3), 0px 1px 3px 1px rgba(60, 64, 67, 0.15);\n",
              "      fill: #174EA6;\n",
              "    }\n",
              "\n",
              "    .colab-df-buttons div {\n",
              "      margin-bottom: 4px;\n",
              "    }\n",
              "\n",
              "    [theme=dark] .colab-df-convert {\n",
              "      background-color: #3B4455;\n",
              "      fill: #D2E3FC;\n",
              "    }\n",
              "\n",
              "    [theme=dark] .colab-df-convert:hover {\n",
              "      background-color: #434B5C;\n",
              "      box-shadow: 0px 1px 3px 1px rgba(0, 0, 0, 0.15);\n",
              "      filter: drop-shadow(0px 1px 2px rgba(0, 0, 0, 0.3));\n",
              "      fill: #FFFFFF;\n",
              "    }\n",
              "  </style>\n",
              "\n",
              "    <script>\n",
              "      const buttonEl =\n",
              "        document.querySelector('#df-1a5005cb-18cd-4b07-85f2-e3680f0e9d28 button.colab-df-convert');\n",
              "      buttonEl.style.display =\n",
              "        google.colab.kernel.accessAllowed ? 'block' : 'none';\n",
              "\n",
              "      async function convertToInteractive(key) {\n",
              "        const element = document.querySelector('#df-1a5005cb-18cd-4b07-85f2-e3680f0e9d28');\n",
              "        const dataTable =\n",
              "          await google.colab.kernel.invokeFunction('convertToInteractive',\n",
              "                                                    [key], {});\n",
              "        if (!dataTable) return;\n",
              "\n",
              "        const docLinkHtml = 'Like what you see? Visit the ' +\n",
              "          '<a target=\"_blank\" href=https://colab.research.google.com/notebooks/data_table.ipynb>data table notebook</a>'\n",
              "          + ' to learn more about interactive tables.';\n",
              "        element.innerHTML = '';\n",
              "        dataTable['output_type'] = 'display_data';\n",
              "        await google.colab.output.renderOutput(dataTable, element);\n",
              "        const docLink = document.createElement('div');\n",
              "        docLink.innerHTML = docLinkHtml;\n",
              "        element.appendChild(docLink);\n",
              "      }\n",
              "    </script>\n",
              "  </div>\n",
              "\n",
              "\n",
              "<div id=\"df-24e0c126-9c66-4ff4-9cdd-ad50ca56c0b0\">\n",
              "  <button class=\"colab-df-quickchart\" onclick=\"quickchart('df-24e0c126-9c66-4ff4-9cdd-ad50ca56c0b0')\"\n",
              "            title=\"Suggest charts\"\n",
              "            style=\"display:none;\">\n",
              "\n",
              "<svg xmlns=\"https://www.w3.org/2000/svg\" height=\"24px\"viewBox=\"0 0 24 24\"\n",
              "     width=\"24px\">\n",
              "    <g>\n",
              "        <path d=\"M19 3H5c-1.1 0-2 .9-2 2v14c0 1.1.9 2 2 2h14c1.1 0 2-.9 2-2V5c0-1.1-.9-2-2-2zM9 17H7v-7h2v7zm4 0h-2V7h2v10zm4 0h-2v-4h2v4z\"/>\n",
              "    </g>\n",
              "</svg>\n",
              "  </button>\n",
              "\n",
              "<style>\n",
              "  .colab-df-quickchart {\n",
              "      --bg-color: #E8F0FE;\n",
              "      --fill-color: #1967D2;\n",
              "      --hover-bg-color: #E2EBFA;\n",
              "      --hover-fill-color: #174EA6;\n",
              "      --disabled-fill-color: #AAA;\n",
              "      --disabled-bg-color: #DDD;\n",
              "  }\n",
              "\n",
              "  [theme=dark] .colab-df-quickchart {\n",
              "      --bg-color: #3B4455;\n",
              "      --fill-color: #D2E3FC;\n",
              "      --hover-bg-color: #434B5C;\n",
              "      --hover-fill-color: #FFFFFF;\n",
              "      --disabled-bg-color: #3B4455;\n",
              "      --disabled-fill-color: #666;\n",
              "  }\n",
              "\n",
              "  .colab-df-quickchart {\n",
              "    background-color: var(--bg-color);\n",
              "    border: none;\n",
              "    border-radius: 50%;\n",
              "    cursor: pointer;\n",
              "    display: none;\n",
              "    fill: var(--fill-color);\n",
              "    height: 32px;\n",
              "    padding: 0;\n",
              "    width: 32px;\n",
              "  }\n",
              "\n",
              "  .colab-df-quickchart:hover {\n",
              "    background-color: var(--hover-bg-color);\n",
              "    box-shadow: 0 1px 2px rgba(60, 64, 67, 0.3), 0 1px 3px 1px rgba(60, 64, 67, 0.15);\n",
              "    fill: var(--button-hover-fill-color);\n",
              "  }\n",
              "\n",
              "  .colab-df-quickchart-complete:disabled,\n",
              "  .colab-df-quickchart-complete:disabled:hover {\n",
              "    background-color: var(--disabled-bg-color);\n",
              "    fill: var(--disabled-fill-color);\n",
              "    box-shadow: none;\n",
              "  }\n",
              "\n",
              "  .colab-df-spinner {\n",
              "    border: 2px solid var(--fill-color);\n",
              "    border-color: transparent;\n",
              "    border-bottom-color: var(--fill-color);\n",
              "    animation:\n",
              "      spin 1s steps(1) infinite;\n",
              "  }\n",
              "\n",
              "  @keyframes spin {\n",
              "    0% {\n",
              "      border-color: transparent;\n",
              "      border-bottom-color: var(--fill-color);\n",
              "      border-left-color: var(--fill-color);\n",
              "    }\n",
              "    20% {\n",
              "      border-color: transparent;\n",
              "      border-left-color: var(--fill-color);\n",
              "      border-top-color: var(--fill-color);\n",
              "    }\n",
              "    30% {\n",
              "      border-color: transparent;\n",
              "      border-left-color: var(--fill-color);\n",
              "      border-top-color: var(--fill-color);\n",
              "      border-right-color: var(--fill-color);\n",
              "    }\n",
              "    40% {\n",
              "      border-color: transparent;\n",
              "      border-right-color: var(--fill-color);\n",
              "      border-top-color: var(--fill-color);\n",
              "    }\n",
              "    60% {\n",
              "      border-color: transparent;\n",
              "      border-right-color: var(--fill-color);\n",
              "    }\n",
              "    80% {\n",
              "      border-color: transparent;\n",
              "      border-right-color: var(--fill-color);\n",
              "      border-bottom-color: var(--fill-color);\n",
              "    }\n",
              "    90% {\n",
              "      border-color: transparent;\n",
              "      border-bottom-color: var(--fill-color);\n",
              "    }\n",
              "  }\n",
              "</style>\n",
              "\n",
              "  <script>\n",
              "    async function quickchart(key) {\n",
              "      const quickchartButtonEl =\n",
              "        document.querySelector('#' + key + ' button');\n",
              "      quickchartButtonEl.disabled = true;  // To prevent multiple clicks.\n",
              "      quickchartButtonEl.classList.add('colab-df-spinner');\n",
              "      try {\n",
              "        const charts = await google.colab.kernel.invokeFunction(\n",
              "            'suggestCharts', [key], {});\n",
              "      } catch (error) {\n",
              "        console.error('Error during call to suggestCharts:', error);\n",
              "      }\n",
              "      quickchartButtonEl.classList.remove('colab-df-spinner');\n",
              "      quickchartButtonEl.classList.add('colab-df-quickchart-complete');\n",
              "    }\n",
              "    (() => {\n",
              "      let quickchartButtonEl =\n",
              "        document.querySelector('#df-24e0c126-9c66-4ff4-9cdd-ad50ca56c0b0 button');\n",
              "      quickchartButtonEl.style.display =\n",
              "        google.colab.kernel.accessAllowed ? 'block' : 'none';\n",
              "    })();\n",
              "  </script>\n",
              "</div>\n",
              "\n",
              "    </div>\n",
              "  </div>\n"
            ],
            "application/vnd.google.colaboratory.intrinsic+json": {
              "type": "dataframe",
              "summary": "{\n  \"name\": \"show(results\",\n  \"rows\": 5,\n  \"fields\": [\n    {\n      \"column\": \"Time\",\n      \"properties\": {\n        \"dtype\": \"number\",\n        \"std\": 1.5811388300841898,\n        \"min\": 26.0,\n        \"max\": 30.0,\n        \"num_unique_values\": 5,\n        \"samples\": [\n          27.0,\n          30.0,\n          28.0\n        ],\n        \"semantic_type\": \"\",\n        \"description\": \"\"\n      }\n    },\n    {\n      \"column\": \"Quantity\",\n      \"properties\": {\n        \"dtype\": \"number\",\n        \"std\": 0.8156004849296761,\n        \"min\": 72.2996253904031,\n        \"max\": 74.36293391475057,\n        \"num_unique_values\": 5,\n        \"samples\": [\n          73.83930457560307,\n          72.2996253904031,\n          73.32091152984704\n        ],\n        \"semantic_type\": \"\",\n        \"description\": \"\"\n      }\n    }\n  ]\n}"
            }
          },
          "metadata": {},
          "execution_count": 10
        }
      ],
      "source": [
        "show(results.tail())"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "formal-headset",
      "metadata": {
        "id": "formal-headset"
      },
      "source": [
        "With `t_0=0`, `t_end=30`, and `dt=1`, the time stamps go from `0.0` to `30.0`.\n",
        "\n",
        "Here's what the `TimeSeries` looks like."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "id": "important-constitution",
      "metadata": {
        "id": "important-constitution",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 309
        },
        "outputId": "4f1e5c99-c6eb-43bf-9070-a10edea9671a"
      },
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 450x300 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAboAAAEkCAYAAABZm/S2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAuJAAALiQE3ycutAABCx0lEQVR4nO3dd1xWdf/H8dfFFhkqiCLDhcgSwY0IuHBkqBCWad45qJ9lZTm6LU1vw5aa3qaIlZbmznCboiioKKWiuEgcKbLcOFFknN8f3l5lgqBwceHF5/l4XI+HHM54H07x4XvO+X6/KkVRFIQQQggdpaftAEIIIYQmSaETQgih06TQCSGE0GlS6IQQQug0KXRCCCF0mhQ6IYQQOk0KnRAaNHnyZGrXro2ZmRlXr17l1KlTtG7dGnNzc0aPHq3teE9NpVKRlJQEwOeff86rr76q3UBClIJK+tEJUbz4+Hg+++wzfvvtNxRFoX79+gwcOJD3338fIyOjJ26blpZGw4YNOXPmDPXr1wcgLCwMRVFYsGCBxjJv2LCBadOmcejQIYyMjGjcuDFDhw5l+PDhZd63SqXi0KFDeHl5lT2oEBVEWnRCFGPjxo307NmT7t27c+rUKa5fv87KlStJTk4mKyurxO3PnTuHmZmZusgBnD17lmbNmmksc2RkJK+//jpvvPEG6enpXLlyhcjISDZu3KixYwpR6SlCiMcUFhYqDRs2VMLDw5+43v79+5X27dsrlpaWiqurq7Js2TJFURRlzZo1iomJiQIo1atXVzp16qS0bt1a0dPTU4yMjJTq1asr27ZtUxRFUZYvX640a9ZMsbS0VFq1aqXs2bNHvf/79+8rn3zyidKoUSOlVq1aSlBQkJKRkVFklps3byrm5ubK4sWLn5g5Ojpa8fLyUiwsLBRvb291jofHGzdunOLg4KBYW1srL7/8snLp0iX19wHl0KFDiqIoyqRJk5Q+ffo88r3IyEjF3d1dMTc3V4KCgpTr16+rv79z507Fw8NDMTMzU4KDg5WhQ4cqr7/++hOzClEepNAJUYSUlBQFUE6fPl3sOtnZ2YqVlZXyzTffKPfv31fi4uKU6tWrK/Hx8YqiKEpsbKxiaWn5yDYBAQHKzJkz1V9v2rRJsbOzUxITE5WCggIlKipKqVWrlnLlyhVFURRl7NixSufOnZXMzEwlNzdXGT16tOLn51dknujoaEVfX1/Jzc0tNvOpU6cUExMTJSoqSsnLy1NWrVqlVKtWTfnzzz8VRVGUyZMnKx4eHkpqaqpy69Yt5ZVXXlECAwPV25dU6Dp16qRcvHhRyc7OVry9vZVJkyYpiqIo165dUywtLZUFCxYoeXl5yqZNmxRjY2MpdKJCyK1LIYpw+fJlAOzs7IpdZ9OmTdSuXZt3330XQ0NDAgICGDBgAIsWLSr1cSIiIhg7diwtWrRAT0+PkJAQXFxc+PXXX1EUhblz5zJjxgxsbW0xMjJiypQp7Nmzh7S0tCIzW1tbP/HZ4cqVK+nYsSMhISEYGBgQGhpKhw4dWL58OQCLFy9mwoQJODo6YmZmxowZM9i2bRuZmZmlOp8PP/wQGxsbatSowUsvvURiYiLw4Dawvb09Q4cOxcDAgBdeeIEuXbqU+uckRFlIoROiCNbW1gBkZGQUu056ejoNGjR4ZFmjRo1IT08v9XHOnTvHxx9/TI0aNdSfpKQkMjIyuHLlCnfu3MHf31/9vbp162JkZFRkobO2tubKlSvcv3//mTP/8/v16tXD2Ni41OdUt25d9b+rV6/OrVu3AMjMzMTBweGRdR0dHUu1TyHKSgqdEEVwdnamQYMGrFixoth17O3tOXfu3CPLzp07h729famP4+DgwNdff83169fVnzt37jBu3DisrKwwNTXl999/f+T7d+/epX379o/tq3379piamrJq1apnzvzP71+4cIHc3NynOqei1KtX77HifP78+TLtU4jSkkInRBFUKhWzZ8/myy+/ZPbs2Vy9ehWAkydPMmzYMFJTU3nhhRe4dOkSc+fOJT8/n927d7N06VL+9a9/lfo4I0aMYNq0aSQmJqIoCjk5OcTExJCeno6enh7Dhw9n9OjR6iJx9epVVq5cWeS+zM3N+eqrr3jvvfdYunQpN2/eRFEUkpKS6N27NwCvvPIKcXFxrFu3jvz8fFavXs2uXbvo378/AK+99hqff/45aWlp3L59m1GjRtG1a1fq1atXlh8nvXr1Ii0tjYULF5Kfn8+WLVvYsWNHmfYpRGlJoROiGC+++CKbN29m06ZNNG7cmBo1ahAaGoqLiwu2trbUrFmTzZs3s2TJEqysrHjzzTeJjIykQ4cOpT5GUFAQX375JW+88QY1a9akYcOGzJo1i8LCQgC++OILfHx86Ny5M+bm5rRs2ZKtW7cWu7+33nqLH3/8kXnz5lGvXj2sra158803efHFFwFwcnJi9erVTJo0iVq1avHpp5+yZs0aGjVqBMBHH31E9+7d8fHxoUGDBuTl5bFkyZIy/BQfqFWrFmvXrmX69OnUqFGD7777jn79+mFsbFzmfQtREukwLoTQiu7du+Pv78/48eO1HUXoOGnRCSEqxNatW7ly5Qr5+fmsWLGCHTt2EBISou1Yogow0HYAIUTVkJiYyMCBA8nJyaFhw4YsX74cV1dXbccSVYDcuhRCCKHT5NalEEIInSaFTgghhE7TuWd0FhYWZe7cKoQQ4vmSnp7OzZs3i/yezhU6e3t7kpOTtR1DCCFEBXJzcyv2e3LrUgghhE6TQieEEEKn6dytSyGE0DXSC+wvKpXqqbeRQieEEJVUYWEhaWlp5OTkaDtKpWFqaoqDgwN6eqW/ISmFTgghKqnLly+jUqlwdnZ+ql/suqqwsJCMjAwuX75MnTp1Sr2dxn9y8+bNw8PDA3d3d95++20KCgoAGDduHE5OTjg7OxMVFVXktvn5+QwZMgQnJyfc3d3ZvXu3puMKIUSlcePGDerUqYO+vj4qlarKf/T19alTpw43btx4qp+jRlt0x48fZ/r06SQmJmJpacnIkSNZunQp9erVY+/evaSkpHDhwgV8fHzo3r07ZmZmj2y/cOFCcnNzOX36NEeOHCE0NJSUlJRnukcrhBDPE0VRKCgowNDQUNtRKhVDQ0MKCgpQFKXUtUCjLbrk5GTatm2LpaUl8GBajpUrVxIVFcXgwYPR19fHzs4OX1/fIufYioqKIiwsDABPT0+srKw4dOiQJiMDcOriLY0fQwghSkP+sH/Us/w8NFroPD09iY+PJzMzk4KCAlatWkVaWhrp6ek4ODio13N0dFTPoPx3pV2vPCWlXaf7f3fx/opD3Libp9FjCSGELlq9ejWurq54eXlRUFDARx99hJubG0OGDNFKHo3eumzatCmfffYZvXv3xtDQkMDAQJKSkkr9qmxp1ouIiCAiIkL9dXZ29jPnBWhax5xB7eqzKCGV389eY3q/5vg6WZdpn0IIUZV89913zJ49m65duwIwa9Ysrl+/jpGRkVbyaPxllNdee40DBw6QkJCAt7c3Li4uODg4PNIyO3/+fJHjU5ZmvREjRpCcnKz+1KxZs0x5qxnpM7mPB4uHtUFRYOD835m84Tj38grKtF8hhHie7dixgzZt2tC8eXNat27NtWvXmDhxIh4eHnh4ePDpp58CMH78eOLj43nrrbcYMmQI3bt35+7du7Rp04aIiAgKCwsZP348bdq0wdPTk48//lh9jLi4ODp06EDLli3p1q1bud3B03j3gosXL1KnTh1u3rzJV199xYwZM7hz5w7h4eG8/vrrXLhwgfj4eL777rvHtg0JCWH+/Pl07tyZI0eOcPnyZVq0aKHpyAD4NalN9Pv+fLLuGD/uOcfuU1eY+bIXzewtK+T4Qgjxdx/+cpiTF2+X+36d65gxNbT5E9e5cuUKgwYNYvv27bi4uHDr1i127NhBTEwMBw4cAMDf35/WrVvz2WefsWfPHqZMmUKHDh0AMDAwICkpCYAffvgBgH379lFYWEjfvn3ZvHkzbdu25eOPP2bz5s1YWlqyatUqxowZw8qVK8t8jhovdP379+fixYsoisLYsWNp3749ANu2bVP3DZkxYwbm5uYAhIWF0bt3b3r37s3gwYOJj4/HyckJIyMjFixYUKEPZi1NDfnmVW8C3eowYe0xgufu4b0uTXi7Y2MM9KVPixCiavjtt99o164dLi4uAJibmxMXF8fAgQMxMTEBYMCAAcTGxtKzZ88n7uvXX3/l8OHDbNq0CYA7d+5w6tQpCgoKSElJISAgAHjQZ+5hXSgrjRe62NjYIpdPnTqVqVOnPrZ8/vz56n8bGhqyaNEijWUrraDm9WjdoBZjfznMjG0n2XHiEjNebk6j2mYlbyyEEOWgpFaXJhX1vsQ/Gx2lbYQoisK0adPo27fvI8s3bNiAn58fa9eufdaYxZJmSSnVtTThp6FtCO/jzokLN+n1TTyLE87JGHRCCJ3n4+PDb7/9xokTJwC4desWHTt2ZNmyZeTm5nLv3j2WL19O586dS9xXz549iYyM5N69ewBkZmaq+1Pv27ePY8eOAZCXl8fRo0fLJb8MAfYUVCoVg3wa4OtkzQc/H+aTdcfZ9sclpoV6UsfCRNvxhBBCI6ytrVm8eDGvvfYaeXl5mJiYsHnzZg4cOEDLli0BePnll+nRo0eJ+xo2bBjp6em0atUKlUpF9erVWbhwIS4uLixfvpywsDDu3r1Lfn4+b7/9Ns2aNStzfpWiY00SNze3Cpl4Nb+gkLlxZ/hm+ymqGxsQ3teD3s3rafy4QoiqQVEUTpw4gYuLi3Qa/5vifi5P+t0vty6fkYG+Hu91acLqt9tT29yY95YfYsSyg2Tfua/taEIIIf5GCl0ZedrXYOO7HRjWoSG/Hs2i2393EXvikrZjCSGE+B8pdOXAxFCfT150Y/kb7TA20GPIwv18tPoIt3PztR1NCCGqPCl05ahdIyu2vO9P/9YOLN+XRs9Zu/j9z6vajiWEeI7p2GsUZfYsPw8pdOXMzNiAL1/yZMHrrbh7v5D+3//GlI3JMoSYEOKpqFQqjI2NuXr1KoWFhSiKUuU/hYWFXL16FWNj46d6QUe6F2hIF9c6bP2gJp+sPcb8+LPsPHmZGTKEmBDiKTwc7/fKlSvajlJpGBsbPzKrTWlI9wINUxSF9Ycz+WTtMXLuF/Bu5ya83akxhjKEmBCilHTs13SZFNeSe9LvfmnRaZhKpaKPlx1tG1rxYdQRZsacJOaPi8x4uTlN6pTPOG5CCN0m/ejKRpoVFaSupQmLhrTms2APzly+Ta/Z8Xy36wwFhfKXmhBCaJIUugqkUqkY2LY+W0b64+VQg89/PcEr3yZw7sodbUcTQgidJYVOCxytTFnxRjs+edGNoxk36DlrNz8lnKNQWndCCFHupNBpiZ6e6sFoKiP9aFrXnInrjvPagt9Jz87RdjQhhNApUui0rHFtM34Z7sOHPZqy/9w1evx3Nyv3n5e3rIQQopxIoasEDPT1eLujE+vf6YBjLVP+HXWUoQv3c/HmPW1HE0KI554UukrE1daCtSN8ea9LE3adukK3mbtYl5QhrTshhCgDKXSVjJGBHqMCnVnzv+l/Rq5IYviSRC7fytV2NCGEeC5JoaukHk7/838BjdiWfJFuM3ey8UimtmMJIcRzRwpdJWZiqM9HPV1ZNbw9NU2NeGfZIUYsPcjV29K6E0KI0tJooVuyZAmenp54eXnh5+dHSkoK8+fPx8vLS/0xMjJi3bp1j20bFxeHubm5er3g4GBNRq3UWtavya8j/Qjr0JBfj2XRbeYuNh/N0nYsIYR4LmhsUOecnBwcHBxISUnB2tqaefPmsX37dlatWqVeJy0tDU9PT7KysjAxMXlk+7i4OKZMmUJMTMxTHbeyDepc3vafu8bYVYc5dzWH3s3rMbm3OzWrG2k7lhBCaNWTfvdrrEX3cP6k27dvA3Djxg1sbW0fWWfZsmWEhIQ8VuRE8Vo3qMXmkf4M8W3A+sOZBM7cxdbjF7QdSwghKi2NFTozMzPmzJmDh4cHdnZ2LFq0iPDw8EfWWbp0KYMGDSp2H4mJiXh7e+Pv7090dHSR60RERODm5qb+ZGdnl+t5VEbVjPSZFOTOijfbYWqkz5uLE/lgZRLXc+5rO5oQQlQ6Grt1mZeXR6dOnfj+++9xdXVl4sSJXLp0iXnz5gFw+PBhgoKCSE1NLXIKips3bwJgYWHB8ePH6dGjB7t27aJhw4ZPPK6u37r8p5z7+Xy1+QSLElKxMTfmi5BmdHGto+1YQghRobRy6zIpKQlFUXB1dQWgf//+7N27V/39JUuWMHDgwGLnWbKwsMDCwgIAd3d3fH19OXjwoKbiPrdMjQyY3MeDZW+0xchAj2GLDjDq5yRu5ORpO5oQQlQKGit09vb2pKSkkJGRAcC2bdtwc3MDHjy/W758+RNvW2ZlZalHBMnIyCAhIQF3d3dNxX3utW9szZb3/RnUrj6rD2YQOHMnMckXtR1LCCG0TmMzjNva2vLll18SGBiIoaEhtWvX5ocffgAevFFpY2OjLnwPTZw4kXr16jF8+HCioqKIjIzE0NAQgClTpuDi4qKpuDrBzNiA8L4e9GxWlw9/OULYTwcI9rZjUpAbNUzlzUwhRNWksWd02lLVntEV505uPl9tOcFPCanUNjfm8+BmBLrJszshhG7SyjM6oV3VjQ34tI8Hy99oh4mhHm/8dID3Vxwi+468mSmEqFqk0Ok4n8ZWbBnpz+s+9Vmb9KDfXbT0uxNCVCFS6KqA6sYP3sx82O/u/xYn8t7yQ1yT1p0QogqQQleFtGtkxZb3/Rjc/sGoKt1m7pQxM4UQOk8KXRVjamTAf3q7s/LNdlQ3NuCtpQd5e2kiV2RGBCGEjpJCV0W1bfTg2V1Yh4ZsPnaBwBk7ZTZzIYROkkJXhVUz0mfCi278Mrw9taobMXJFEm8uTuTSzXvajiaEEOVGCp2gZf2abHrPj7c6NmbHiUt0nbGTVQfSpHUnhNAJUugE8GA283/3cGHN2+2pV6MaY385wuAf95N5/a62owkhRJlIoROP8LSvwfp3OvB+1ybsOX2FbjN3sez389K6E0I8t6TQiccYGejxfldnNrzbgQbWpny85igD5/9O2rUcbUcTQoinJoVOFMvV1oK1b/sytntTDpzLptvMXfy45ywFhdK6E0I8P6TQiScy0NdjRCcnfh3ZARdbcyZvSOblbxM4femWtqMJIUSpSKETpeJkY84vw9sz8UU3kjNv8sKseCJiT5NXUKjtaEII8URS6ESp6eupGNqhIdHv+9O6YU2mRafQZ84ejmXc0HY0IYQolhQ68dQcrUxZMqwtX4Y0I+1aDn0i9jB1ywnu5RVoO5oQQjxGCp14JiqViv5tHNk2KoBOTW2YG3eGF77ZzYFz17QdTQghHiGFTpRJXUsTvv9XS7551ZvrOXn0+zaB/6w/zp3cfG1HE0IIQAqdKAcqlYrezeux7QN/ejevx8K95+g2cxe7Tl7WdjQhhJBCJ8qPlZkxs/p7M/9frSgoVPjXD/sY/fNhrufIBK9CCO3ReKFbsmQJnp6eeHl54efnR0pKCnFxcZibm+Pl5YWXlxfBwcHFbj9u3DicnJxwdnYmKipK03FFOejqVoeto/wZ2NaRqIPpdJ2xk01HsmQYMSGEVmi00OXk5DBy5Eh27NhBUlISAwcOZMKECQC0bduWpKQkkpKSWLNmTZHbx8TEsHfvXlJSUoiNjeWDDz7g9u3bmowsyomFiSGfBTdjxZvtMDcxZMSyg7y5OJELN2QKICFExdJooSssLERRFHVxunHjBra2tqXePioqisGDB6Ovr4+dnR2+vr5s3bpVU3GFBrRrZMXmkX9NARQ4YyfLfj9PoQwjJoSoIBotdGZmZsyZMwcPDw/s7OxYtGgR4eHhACQmJuLt7Y2/vz/R0dFFbp+eno6Dg4P6a0dHR9LS0jQZWWjAwymA1o3wpf7/Bol+9fvfOHvljrajCSGqAI0Wury8PObOncv+/fvJyMggNDSUf//737Ro0YLU1FQOHTpEZGQkYWFhnD179rHtS/NMJyIiAjc3N/UnOztbE6ciyoGHnSVr3/ZlXE8XktKu0+O/u4iMO0O+DCMmhNAgjRa6pKQkFEXB1dUVgP79+7N3714sLCywsLAAwN3dHV9fXw4ePPjY9g4ODo+04M6fP4+9vf0j64wYMYLk5GT1p2bNmho8I1FWBvp6DA9ozJb3/fF2rMFXW07QJ0KGERNCaI5GC529vT0pKSlkZGQAsG3bNtzc3MjK+usNvIyMDBISEnB3d39s+5CQEBYuXEhBQQEZGRnEx8fTrVs3TUYWFaShdXWWhbXji5BmnL/6YBixLzb/wd37MoyYEKJ8GWhy57a2tnz55ZcEBgZiaGhI7dq1+eGHH4iKiiIyMhJDQ0MApkyZgouLCwATJ06kXr16DB8+nMDAQLZt24azszN6enrMmDEDc3NzTUYWFUhPT8WrbRzp1NSGieuO8e3OP9ly7AKfBzfD18la2/GEEDpCpehY5yY3NzeSk5O1HUM8gy3Hsvhk3XEu38qlX0t7xvdypYapkbZjCSGeA0/63f/EFl1OTg6rV69m9+7dZGZmUq1aNTw8PAgKCsLb21sjYUXV1cPDFp/G1ny5+QTL950nNuUSE4PcCfK0RaVSaTueEOI5VWyLLjw8nJ9//pnAwEBatWqFjY0N9+7d4+TJk2zfvp3c3FxmzZpV5LM1bZIWnW74/c+rfLT6KH9euUNnFxvC+3pgV6OatmMJISqpZ2rROTs7c+TIkSL/kh41ahTp6emkpqaWX0oh/qZtIyt+HelHROxpIuPO0G3GTsZ2b8ognwbo60nrTghRevKMTlR6Jy7cZFzUUZLSruPlUIOvXvKkaV15KUkI8Zcn/e4vtnvBRx99xNy5cx9bHhkZyccff1x+6YQogUtdC6Leas+kIDdOXrxFr292Mz06RWY0F0KUSrEtulatWrFv3z709B6thQUFBXh5eXH06NEKCfi0pEWn29Kzc/hk7TFiUy7T0Lo6nwV70L6xdEUQoqp7phZdbm7uY0UOQF9fX6ZbEVpjX9OUHwa3Zvar3ty6l8eA739n7KrDZN+ROe+EEEUrttAVFBSQmZn52PLMzEwKC2VsQqE9KpWKoOb12D6qI6+2cWBV4oM579YlZcgfYUKIxxRb6N555x169epFbGwsN2/e5ObNm8TGxhIUFMSIESMqMqMQRbI0NeSLEE9WvtkOS1NDRq5I4l8/7OP81RxtRxNCVCJPfOtywYIFhIeHqwdWdnR0ZPz48YSFhVVYwKclz+iqptz8AiLjzjA39gx6evBBV2eGdWiIgb5Gh3MVQlQST/rdX6ruBQ8nTjUzMyvfZBogha5qO33pFh+vPsa+c9dws7Xgi5BmNHeooe1YQggNe6aXUf4+bY6ZmdljRe7evXukpKSUU0QhyoeTjTkr3nwwK0Jadg7Bc/cwecNxbufmazuaEEJLih0ZZdasWWRmZhISEvLIEGCnTp0iJiaGvXv38vXXX9O0adOKzCtEiR7OitDF1YZPNyTz455zRB+7wOQ+HgS61dF2PCFEBXvircvTp0+zePHiIgd1Dg4OVk+zU5nIrUvxT7EnLjFh7TEyrt+lu3sd/tPbHVtLGTdTCF1S5md0zxMpdKIoOffzmRVzivnxZzEx0GNM96b8S8bNFEJnPNMzOiF0iamRAR+94MqGdzrgVMecyRuSCZ67h2MZN7QdTQihYVLoRJXiVs+C1W+159M+7vx5+Q6958QzZWMyd+RlFSF0lhQ6UeXo66n4l08Dto8OoIdHXebHn6XbzF1s/+OitqMJITSgVIUuOjqaOXPmAHDx4kVOnjyp0VBCVIQ6FibMHdiSBa+3AmDYogO8tSSRCzfuaTmZEKI8lVjoJkyYwNy5c5k9ezYAhYWFDB48WNO5hKgwXVzrsPUDf97wa8jW5It0nbGTH/ecpaBQp97TEqLKKrHQbdy4kbVr12JqagqAra2teqQUIXRFdWMDxvdyY90IXxrbmDF5QzJ9I/ZwJP26tqMJIcqoxEJnYmKCSvXXK9i5ubmlHiF+yZIleHp64uXlhZ+fHykpKSxdupTmzZvj6elJq1at2LFjR5HbxsXFYW5ujpeXF15eXgQHB5fylIR4dh52lqx+qz3hfdw5d+UOfSP28J/1x7l1L0/b0YQQz6jYkVEe6tq1K5988gk5OTmsX7+eyMhIQkNDS9xxTk4OI0eOJCUlBWtra+bNm8eECRP44IMP2LFjB1ZWVhw9epSuXbuSlZVV5Nx3bdu2JSYm5tnOTIhnpK+nYpBPA7q71yV80x8s3HuOX49mMSnInRea1X3kDz8hROVXYosuPDycRo0a4enpyU8//URoaCiTJk0qcceFhYUoiqK+zXnjxg1sbW1p3749VlZWAHh4eJCbm8udO3fKeBpClD8bCxNmv+rNoqFtMDHUZ8SygwxZuF+mARLiOfPEFl1BQQHdunVj+/btDBky5Kl2bGZmxpw5c/Dw8MDS0hJLS0sSEhIeWWfFihU0a9YMc3PzIveRmJiIt7c35ubmjB8/nu7duz+2TkREBBEREeqvs7OznyqnECUJcK7N1g/8iYg9zbydZwicuZP3ujThDb9GGBlIDx0hKrsShwALCAhgy5YtVKv2dGMD5uXl0alTJ77//ntcXV2ZOHEily5dYt68eQAcOnSIvn37EhMTQ5MmTR7b/ubNmwBYWFhw/PhxevTowa5du2jYsOETjytDgAlNOn3pNuPXHOX3s9doYmPGZ8HNaNOwlrZjCVHlPel3f4nP6BwcHGjTpg29e/emevXq6uUff/zxE7dLSkpCURRcXV0B6N+/P/379wfg5MmTvPTSS6xYsaLIIgcPCtxD7u7u+Pr6cvDgwRILnRCa5GRjxoo32xF1MIPPf/2Dl79NoF9Lez56wZVa1Y20HU8IUYQS77s0adKE0NBQjIyMyMvLU39KYm9vT0pKChkZGQBs27YNNzc30tPT6dWrF5GRkfj4+BS7fVZWlvrtzoyMDBISEnB3dy/teQmhMSqVitCW9mwfFcArrRxYlZhO56/jWLHvPIXS906ISkejsxfMnz+fGTNmYGhoSO3atfnhhx8IDw/n559/fqRltnbtWho0aMDEiROpV68ew4cPZ86cOURGRqqnAho9ejSDBg0q8Zhy61JUtAPnrjFh7TFOXLhFy/o1mdLXA1dbi5I3FEKUmzJN0/Ppp58WuXzixIllT6YBUuiENuQVFLJwzzlmxpwkN7+QIe0b8H6gM2bGJT4dEEKUgzJN06Ovr6/+5OXlsW7dOs6ePVvuIYV4nhnq6/GGfyNiRgUQ6FqH+fFn6fr1TjYfzSr1AAtCCM146luXOTk5BAUFsX37dk1lKhNp0YnKIPbEJSauP0batbt0alqbyb09cLQy1XYsIXRWuU68qq+vz/nz58scSghd1snFhq3vB/BOJyfiT18hcOZOZm8/RW5+gbajCVHllPgAITAwUD3kUWFhISdOnOD111/XeDAhnnfVjPQZ070pfb3t+GTtMb7edpI1SRmE9/HA18la2/GEqDJKvHW5c+dO9b8NDAxo0KABdnZ2Gg/2rOTWpaiMFEVhXVImUzb9wZXbuQQ1r8eEXq7UsTDRdjQhdEKZbl1GR0cTEBBAQEAAvr6+2NnZldhZXAjxKJVKRV9vO3aMCWBw+wZsOpJJl693siD+LPkFhdqOJ4ROK7HQbdmy5bFlv/76q0bCCKHrLEwM+U9vd9a/0wEnGzPCNybz4ux4Dpy7pu1oQuisYp/RLV68mJ9++okzZ87QrVs39fKbN29iY2NTIeGE0FUP5737+UAaX245Qei8BEJb2vNRTxeszIy1HU8InVJsofP398fBwYFRo0Yxfvx49XJzc3OaN29eIeGE0GV6eir6t3Gkm3tdpm45wYr9aWxLvsiHPZrSv7Uj+noy750Q5UGjQ4Bpg7yMIp5XB89n88naYxzPvElze0um9G1GM3tLbccS4rlQppdRUlNTGTBgAJ6enjg7O6s/Qojy1cKxJuvf6cDk3u78efkOvSPimbD2KNdz7ms7mhDPtRIL3ZAhQ+jXrx8FBQWsXr2a7t27M2DAgIrIJkSVo6+n4vX2Ddg+JoC+XnYs+e08nb/eyc/702RmBCGeUYmF7saNGwQHB6Ovr4+HhwezZ89m06ZNFZFNiCrLxtyEma94sfLNdlibGfFh1BFemreXYxk3tB1NiOdOiYXOwODB+yrW1tbs27eP9PR0rl69qvFgQgho28iKTe/5MaGXK6cu3qb3nHgmrjvGjZyS54QUQjxQ4hBgb731FtnZ2Xz22Wf079+fW7duMXXq1IrIJoTgwcwIYX6NCGpej89//YOfElLZdCSLcT1deKmFPXrydqYQT/TEQldQUICZmRk1a9bEx8eHP//8s6JyCSH+oY6FCbP6e9O/tSMT1x1j7C9HWLE/jU/7uONeT97OFKI4T7x1qa+vzxdffFFRWYQQpeDT2IpfR/ox/gVXTmTdJGh2PJPWHePGXbmdKURRSnxG98ILL7Bo0SLu35dXnIWoLB5O9Lp9dEd6edZjUUIqXb6OY9UBeTtTiH8qscO4oaEhBQUP5tAyMDBAURRUKlWlLXzSYVxURXvPXGHiuuOcvnSbFo41+LSPBx52cjtTVB1l6jCel5dHYWEhhYWF3L9/n7y8vEpb5ISoqto3tmbz/25nply4RdCceMavkc7mQkApZxiPjo5mzpw5AFy6dImTJ0+W+gBLlizB09MTLy8v/Pz8SElJAWDcuHE4OTnh7OxMVFRUkdvm5+czZMgQnJyccHd3Z/fu3aU+rhBVzcPbmbFjOtLXy46lv5+n0/Q4lv1+ngK5nSmqsBIL3YQJE5g7dy6zZ88GHryJOXjw4FLtPCcnh5EjR7Jjxw6SkpIYOHAgEyZMICYmhr1795KSkkJsbCwffPABt2/ffmz7hQsXkpuby+nTp1m+fDnDhg1Dx4bmFKLc2Vg86Gz+8//5UMfChI/XHCV47h4Onc/WdjQhtKLEQrdx40bWrl2LqakpALa2tkUWpaIUFhaiKIp6/Rs3bmBra0tUVBSDBw9GX18fOzs7fH192bp162PbR0VFERYWBoCnpydWVlYcOnSo1CcnRFXWpmEtNr77YOzMs1fuEDx3Lx/+cpgrt3O1HU2IClVioTMxMUGl+qtDam5ubqlbVWZmZsyZMwcPDw/s7OxYtGgR4eHhpKen4+DgoF7P0dGRtLS0x7Yv7XpCiKIZ6OvxevsGxI7pyMut7Pn5QDqdp8exaO85mdlcVBklFrquXbvyySefkJOTw/r16+nbty+hoaGl2nleXh5z585l//79ZGRkEBoayr///e9SF8rSrBcREYGbm5v6k50tt2eE+CdrM2OmhjZnzdvtqW9VnUnrj/Pi7Hh+/1OG8xO6r8RCFx4eTqNGjfD09OSnn34iNDSUSZMmlWrnSUlJKIqCq6srAP3792fv3r04ODg80jI7f/489vb2j21fmvVGjBhBcnKy+lOzZs1SZROiKvJ2rMnaEb58HtyMizfv8cp3v/He8kNk3bir7WhCaEypJ169e/fB/wjVqlUr9c6zsrJo1qwZhw8fxs7OjlmzZpGQkMCwYcMIDw8nNjaWCxcu0K5dO5KTkzE3N39k+++//57Y2FiWLVvGkSNHCAkJ4dSpU4/cSv0n6UcnROlcz7nPjG0nWfJbKiaG+ozo5ESYX0OMDfS1HU2Ip1amfnTHjh2jdevWODg44ODgQJs2bTh27FipDmxra8uXX35JYGAgzZs3Z8OGDUydOpXAwEDatWuHs7MzHTt2ZMaMGeoiFxYWxvr16wEYPHgwhoaGODk50b9/fxYsWPDEIieEKL0apkZ82seDje/64WFnybToFLrP3MWOExe1HU2IclVii65Vq1ZMnDiR3r17Aw/ewvzPf/7DgQMHKiTg05IWnRBPT1EUNhzJ4vNNf3Dh5j06u9jwyYtuNLSuru1oQpRKmVp09+/fVxc5gBdffJG8PBk8VghdolKp6N28HttHB/B2x8bEn7pC95m7+GrLCe7k5ms7nhBlUmKhCwoK4ttvv+Xu3bvcvXuX77777pHCJ4TQHdWNDfiwhwtbP/DHr4k1kXFn6PL1TtYlZchgDeK59VSDOj+kr//gYXVlHNxZbl0KUX5iT1xi8objnLuaQ5sGtZjU203mvhOVUrkN6vzwk5eXJ4M7C1EFdHKxIfoDfz7s0ZRjmTcImh3Px2uOcu2O/L8vnh9PnGH8ofv373PhwgXy8/+6V9+oUSONhRJCVB7GBvq83dGJl1rY89XmEyz7/TwbD2fyQaAzr7Wrj6F+qcaGF0JrSrx1+d///pfw8HBq166tfrVfpVJV2tuDcutSCM1KTM1m8objHEm/QRMbMyYFudOhibW2Y4kq7km/+0ssdE5OTuzbt49atWppJFx5k0InhOYVFir8kpjO1OgTXLl9n+7udZjQyw2HWqbajiaqqDI9o2vUqBFmZmblHkoI8fzS01PxcmsHdozpyBt+Ddn+xyW6zNjJ11tTyLkv3RFE5VJiiy4pKYl3330Xf39/jI2N1csnTpyo8XDPQlp0QlS805duE74xmZ0nL2NracK4ni70bl5PRjISFaZMLboxY8ZgY2ODiYkJ+vr66o8QQjzkZGPGwiGtWfB6K4wM9Bi5Iol+8xI4mn5D29GEKPmty6ysLGJiYioiixDiOaZSqejiWocOTaz5If4cc3acondEPP1a2jOme1NszE20HVFUUSW26Pz9/UlMTKyILEIIHWBsoM9bHRsTO7YjoS3sWZWYTufpO/l25xly8wtK3oEQ5azEZ3SNGzcmNTUVOzs7jI2NURQFlUrFyZMnKyrjU5FndEJULkfSrzN5QzKJqdk0sDJlQi83urjayPM7Ua7K1L0gNTW1yOX169cvezINkEInROWjKArrD2fy5eYTZN24h18Taya+6EaTOuYlbyxEKZTpZZT69etz9+5dkpKSqF+/PlZWVo9NkCqEEE+iUqno42XH9tEBvNelCfvOXqPHrN38Z/1xrufIcGJCs0osdHPnzmXQoEGMHj0agAsXLhAaGqrxYEII3WNqZMCoQGe2jw6gp0ddFu49R8fpcfyUcI78gkJtxxM6qsRC9+2337Jnzx51K87JyYnLly9rPJgQQnfZ1zRlzoAW/Px/PtjVqMbEdcfpOWs3u07K7xZR/ootdHv27AGgWrVqGBkZqZcXFspfXUKI8tGmYS3Wv9OBr15qRnbOff71wz6GLdzPmcu3tR1N6JBiC927774LgLe3N/Pnzyc/P58jR44QFhZGp06dKiygEEK36eupeKW1I7FjOjI8oDG7/ze7+acbkrmRk6fteEIHFFvoHr6MOXPmTNLS0jAxMWHYsGHY29szffr0CgsohKgazE0MGdfThZhRAQS61eGHPWfpOD2WxfL8TpRRsd0LrK2tCQkJAf4qeuqNVCq+++47zad7BtK9QAjdkHDmKuEbk0nOuolzHTM+edENvya1tR1LVFJP+t1f7BBg1atXx9fXt0wHvnTpEt26dVN/feHCBXx8fOjVqxdz5sxRL09OTmbVqlX06dPnke3j4uIICgqicePGADRs2JA1a9aUKZMQ4vng09iKDe924JfENKZFpzBowT66uNjwcS9XGteWGVVE6RXbomvRogUHDx4s14N17dqVoUOHMmDAAPWytLQ0PD09ycrKwsTk0bHw4uLimDJlylONtSktOiF0z617eUTEnuGH+LMUKgqDfOozsksTapgalbyxqBKeqcN4CQOmPLXMzEwOHDhA3759H1m+bNkyQkJCHityQgjx0N+f33V1rcOPex70v/txz1ny5PmdKEGxhW737t3leqDly5fTt29fTE0fnYF46dKlDBo0qNjtEhMT8fb2xt/fn+jo6Me+HxERgZubm/qTnZ1drrmFEJWHo5Up8wa1ZOWb7bCvWY3JG5Lp/t9dbP/jYrn/cS50R4ljXZYXb29vpk2bRteuXdXLDh8+TFBQEKmpqUUO8Hrz5k0ALCwsOH78OD169GDXrl00bNiw2OPIrUshqobCQoXVhzKYFn2Cizdz6eBkzYQXXXGpa6HtaEILyjTWZXlITk7m8uXLdO7c+ZHlS5YsYeDAgcWOYm5hYYGFxYP/aN3d3fH19S3354ZCiOeTnp6K0Jb2xI7pyHtdmnAg9RovzNrNR6uPcvlWrrbjiUqkQgrd4sWLGTBgAHp6fx2usLCQ5cuXP/G2ZVZWlvp2REZGBgkJCbi7u2s8rxDi+fFw/MzYMR3p42XH8n3n6TQ9jrlxp7mXJ/PfiQoodIqisGzZMl577bVHlsfFxWFjY4Obm9sjyydOnMi8efMAiIqKwsPDAy8vL3r16sWUKVNwcXHRdGQhxHPI1rIaM1/xYt0IX1zqmjN1SwpdZ+xk45FMeX5XxVXYM7qKIs/ohBCKovDr0Qt8sfkP0rPv0sKxBhNedKOFY01tRxMaovVndEIIUZFUKhW9PG2JGRXARz1dOHXxNiFz9/LOsoOkXcvRdjxRwaTQCSF0lomhPv8X0Ji4sR35l099Nh+7QJcZO/li8x/cvCcDRlcVUuiEEDrPysyYT/t4EP2+H35O1ny78086TouTAaOrCCl0Qogqw8nGnAWDW7M0rC11LEz4ZN1xuv93FztOSIdzXSaFTghR5fg6WbPx3Q5MfcmTW/fyGbrwAK8t+J3kzJvajiY0QAqdEKJK0tdT8XJrB3WH88TUbHrN3s3YVYe5cOOetuOJciSFTghRpVU3ftDhPG5MJ15qYc8vB9PpND2OGdtOcic3X9vxRDmQQieEEEBdSxOm92vOxnc70KJ+Db7ZfoqAaXEs33deXlh5zkmhE0KIv3GvZ8mSYW35cUhralU35KPVR3nhm93EplySF1aeU1LohBDiH1QqFZ2a2vDre358EdKMa3fyGPLjfgYt2MfxzBvajieekhQ6IYQohoG+Hq+2cWTn2L9eWHlxdjxjVh0m68ZdbccTpSSFTgghSvDwhZXYMR3p19KeqP+9sDI9OoXb8sJKpSeFTgghSqmupQlTQ5vz63t+tG5QizmxpwmYGsvihHPkyQsrlZYUOiGEeEquthYsHtaWn4a2oba58YMRVmbuIvr4BXlhpRKSQieEEM/I37k2m97zY3q/5tzNK+D/Fify8rcJHDyfre1o4m+k0AkhRBno66kIbWlP7JiOfNijKSeybhEydy9vL03k3JU72o4nkEInhBDlwsRQn7c7OhE3tiOD2zdg6/GLdJ2xk/+sP861O/e1Ha9Kk0InhBDlyMrMmP/0didmVADd3euycO85AqbGEhF7mnt5BdqOVyVJoRNCCA1oYF2diIEtWP12e5rWNWdadAqdpsfx84E0CgrlhZWKJIVOCCE0qIVjTVYN9+HbQS2pZqTPh78c4YVZu2UOvAqk0UJ36dIlvLy81J+6desSHBxMXFwc5ubm6uXBwcHF7mPcuHE4OTnh7OxMVFSUJuMKIYRGqFQqurvXZev7/nwW7MG1nPsMXXiA/t/9RlLadW3H03kGmty5jY0NSUlJ6q+7du1Kv379AGjbti0xMTFP3D4mJoa9e/eSkpLChQsX8PHxoXv37piZmWkythBCaISBvh4D29Yn2NuO+bvP8u3OM/SN2EMvT1vGdmtKA+vq2o6okyrs1mVmZiYHDhygb9++pd4mKiqKwYMHo6+vj52dHb6+vmzdulVzIYUQogKYGhnwXpcm7PywE6/71Cf62AW6ztjJpHXHuHI7V9vxdE6FFbrly5fTt29fTE1NAUhMTMTb2xt/f3+io6OL3CY9PR0HBwf1146OjqSlpVVIXiGE0DRrM2Mm9/Fg2//e0FyUkErA1Fhmbz9Fzn0ZQ7O8VFihW7JkCa+99hoALVq0IDU1lUOHDhEZGUlYWBhnz559bJvSPKiNiIjAzc1N/cnOlhEJhBDPl4b/e0Nz7QhfPOws+XrbSQKmxbH091QZQ7McVEihS05O5vLly3Tu3BkACwsLLCwsAHB3d8fX15eDBw8+tp2Dg8MjLbjz589jb2//yDojRowgOTlZ/alZs6YGz0QIITTHy6EGK95sxw+DW1HT1JDxa47RfeYuNh/Nkjc0y6BCCt3ixYsZMGAAenoPDpeV9ddFy8jIICEhAXd398e2CwkJYeHChRQUFJCRkUF8fDzdunWriMhCCKEVKpWKzi512DzSn2mhntzLK+CtpQfpO3cvCWeuajvec0mjb13Cg9uPy5YtY8OGDeplUVFRREZGYmhoCMCUKVNwcXEBYOLEidSrV4/hw4cTGBjItm3bcHZ2Rk9PjxkzZmBubq7pyEIIoXX6eir6tXIgqHk9FiekMif2NK9+/xsdm9bmw+4uuNWz0HbE54ZK0bH2sJubG8nJydqOIYQQ5erG3Tzm7TzDj3vOkptfSF8vO0YFOuNQy1Tb0SqFJ/3ul5FRhBDiOWBZzZB/93Ahbkwn+rd2YF1SBl2+3snkDce5Kl0SnkgKnRBCPEfqWprwRYgnWz8IoJNLbX7cc46AaXF8s/0Ud3KlS0JRpNAJIcRzyMnGjG8HtSLqrfa41bNgxv+6JPyUcI77+dIl4e+k0AkhxHOsZf2arHyzHT8Obo21mRET1x2n64ydrEvKoFBmSQCk0AkhxHNPpVLRycWGX9/z47+veKGgMHJFEr1mxxN74lKV74MnhU4IIXSEnp6Kvt52bB/Vkcm93bl86x5DFu7nle9+IzG16o4aJYVOCCF0jJGBHq+3b8DOsZ0YFehMcuZNXorcyxs/HeDkxVvajlfhpNAJIYSOqm78YJaEXR92YliHhuxMuUz3/+5i9M+HSc/O0Xa8CiMdxoUQoorIuH6XmdtOsvpgOgZ6egxo68g7nZ2wNjPWdrQykw7jQgghsKtRjen9mhP9vj+dXGqzcO85/KfG8vXWFG7ey9N2PI2RQieEEFVMkzrmfDuoFWtH+OLtWIPZO07jPzWWb3ee4V5egbbjlTspdEIIUUV5OdRgaVg7loa1pb5Vdb7YfIKAabE6Nw+eFDohhKjifJ2sWft2e+a91hILkwfz4AXqUKdzKXRCCCFQqVT08KjLlvf9md6vOXkFf3U633Hi4nPd6VwKnRBCCDV9PRWhLe3ZMSbgf53Ocxm68AD95iXw+5/P58SvUuiEEEI8xthAn9fbN2DXhx0Z270pJy/e4pXvfmPQgt85kn5d2/GeihQ6IYQQxTI1MmBEJyd2f9iZEZ0ac+BcNr3n7OH/Fj8/o6xIoRNCCFEiS1NDxnZ3YdeHnRji24DYEw9GWRm1MonzVyv3KCsyMooQQoinlnH9Lt/EnOKXg+mogFdaO/Bu5ybUtTTRSp4n/e6XQieEEOKZ/Xn5NjNjTrHhcCbG/xtMenhAY2pVN6rQHDIEmBBCCI1oVNuM2a968+t7fvg1sea7XX/iPzWWGdtOVpphxTRa6C5duoSXl5f6U7duXYKDg1m6dCnNmzfH09OTVq1asWPHjiK3j4uLw9zcXL19cHCwJuMKIYR4Rm71LJj/emui3mpPcwdLvtl+Cr+vYomIPU3O/XytZjPQ5M5tbGxISkpSf921a1f69etHgwYN2LFjB1ZWVhw9epSuXbuSlZWFnt7jdbdt27bExMRoMqYQQohy0rJ+TZaGtWPvmStMj05hWnQKP+45y1sdnRjY1hETQ/0Kz1Rhty4zMzM5cOAAffv2pX379lhZWQHg4eFBbm4ud+7cqagoQgghNKx9Y2ui3mrPj4NbU8fChPCNyXScFsfS31O5n1+x42hWWKFbvnw5ffv2xdTU9JHlK1asoFmzZpibmxe5XWJiIt7e3vj7+xMdHf3Y9yMiInBzc1N/srOr7nTxQghRmahUKjq52LDx3Q7Me60F5iYGjF9zjC4z4vglMZ2CChpHs8LeuvT29mbatGl07dpVvezQoUP07duXmJgYmjRp8tg2N2/eBMDCwoLjx4/To0cPdu3aRcOGDYs9jrx1KYQQlVNBocKGw5nMjDlJ6tUcGteuzgeBzrzgYYuenqpM+9b6W5fJyclcvnyZzp07q5edPHmSl156iRUrVhRZ5OBBgbOwsADA3d0dX19fDh48WBGRhRBClDN9PRV9ve2IGRXAlyHNuHu/gHeWHeKHPWc1etwKKXSLFy9mwIAB6pdN0tPT6dWrF5GRkfj4+BS7XVZWlnrE7IyMDBISEnB3d6+IyEIIITTEUF+P/m0ciR3bkU/7uPNSC3uNHk/jhU5RFJYtW8Zrr72mXjZ58mQuXbrEv//9b3XXgXPnzgEwceJE5s2bB0BUVBQeHh54eXnRq1cvpkyZgouLi6YjCyGEqADGBvr8y6cBNTXcuVxGRhFCCPHc0/ozOiGEEEJbpNAJIYTQaVLohBBC6DQpdEIIIXSaFDohhBA6TQqdEEIInSaFTgghhE7TuX50FhYW2NuXvZd9dnY2NWvWLIdElZucp26R89Qtcp6ll56erh4f+Z90rtCVl6rS8VzOU7fIeeoWOc/yIbcuhRBC6DQpdEIIIXSaFLpijBgxQtsRKoScp26R89Qtcp7lQ57RCSGE0GnSohNCCKHTpNAJIYTQaVLo/iEuLg53d3ecnJwICwujoKBA25E0okGDBri7u6snvj169Ki2I5WLkSNHYm9vj4GBwSPLx40bh5OTE87OzkRFRWkpXfkp6jzj4uIwNzdXX9Pg4GAtJiwfaWlpdOnSBVdXV9zd3fnoo4/U39Ola1rceeriNe3WrRteXl40a9aM0NBQdd83jV5PRagVFBQojRs3Vo4fP64oiqL069dPWbhwoZZTaUb9+vWVtLQ0bccod7t371aysrIUfX199bJt27Ypfn5+Sn5+vpKenq44ODgot27d0mLKsivqPGNjY5UuXbpoMVX5y8zMVPbv368oiqLk5uYqHTp0UNauXatz17S489TFa3r9+nX1v0eOHKlMmjRJ49dTWnR/s3//furVq4ebmxsAw4YNe+7/UqxqOnToQN26dR9ZFhUVxeDBg9HX18fOzg5fX1+2bt2qpYTlo6jz1EW2tra0atUKACMjI7y9vTl//rzOXdPizlMXWVpaAlBYWMi9e/dQqVQav55S6P4mPT0dBwcH9deOjo6kpaVpMZFmBQUF4eXlxfjx48nLy9N2HI2pStc1MTERb29v/P39iY6O1naccnXt2jXWrl1LYGCgTl/Tv58n6OY1DQ4OxsbGhpSUFEaPHq3x6ymF7m+UKtTTYvfu3Rw6dIg9e/aQkpLC9OnTtR1JY6rKdW3RogWpqakcOnSIyMhIwsLCOHv2rLZjlYv79+8TGhrKyJEjcXFx0dlr+s/z1NVrumbNGjIzM7G3t+eXX37R+PWUQvc3Dg4Oj/wVcf78+XIZILoyevjXU/Xq1QkLC2Pv3r1aTqQ5VeW6WlhYYGFhAYC7uzu+vr4cPHhQy6nKrqCggAEDBuDl5cXo0aMB3bymRZ2nrl5TeHCLtn///qxZs0bj11MK3d+0atWK9PR09eCiCxYsICQkRMupyt+dO3fUbzoVFBQQFRWFp6enllNpTkhICAsXLqSgoICMjAzi4+Pp1q2btmOVu6ysLPVfxhkZGSQkJODu7q7lVGX35ptvYm5uztdff61epovXtKjz1LVreuvWLbKysoAHz+jWr1+Pu7u7xq+nQcmrVB36+vrMnz+f0NBQcnNzCQgIYNCgQdqOVe4uXrxISEgIhYWFFBQU4OPjw/jx47Udq1z83//9H5s2baKgoAB7e3v69OlDREQE27Ztw9nZGT09PWbMmIG5ubm2o5ZJUefp6upKZGQkhoaGAEyZMgUXFxctJy2bPXv28MMPP+Dh4YG3tzcAQ4cO5b333tOpa1rceerp6enUNb116xZ9+vQhNzeXwsJC2rZty4QJEzA1NdXo9ZQhwIQQQug0uXUphBBCp0mhE0IIodOk0AkhhNBpUuiEEELoNCl0QgghdJoUOiE05NatW+pR5x0dHbGyssLLywtnZ2c6deqksePu2LGDd99996m2OXDgAMOHDy9xvbfeeotdu3Y9azQhtEK6FwhRARYuXEh8fDzz58/X+LH8/f356aefaNCgQbnv+/Tp07z55pvs2LGj3PcthKZIi06ICnbu3DmcnJzU/27UqBHDhw/Hzc2Nzp07c+jQIQIDA2nUqBEzZ85UbxcXF0eHDh1o2bIl3bp1K3LQ2zNnznD37l11kVu4cCFBQUH07NmTRo0aMWrUKNavX0/79u1p3LixunUWFxdH165d1f/29fVlwIABuLm50bNnT+7duweAk5MT2dnZpKamavJHJES5kkInhJadO3eOoUOHkpycTI0aNXjnnXdYv349+/fvJzw8nHv37nHt2jU+/vhjNm3aRGJiIm+88QZjxox5bF/x8fHq6V4eSkpKYtmyZSQnJ7Ny5Uri4uLYu3cvs2fPZtKkSUVmSkpKYsqUKSQnJ6Ovr//IdFWtWrUiPj6+fH8IQmiQDAEmhJbZ2dnRpk0bALy8vLh79y7VqlWjWrVqWFtbc+HCBY4dO0ZKSgoBAQHAg3ECixoiKTMzExsbm0eWdezYkZo1awLQtGlTunfvrj5WcSPht2jRgkaNGgHQunXrR9azsbEhMzOzjGctRMWRQieElhkbG6v/raen99jX+fn5KIqCn58fa9eufeK+qlWrxvXr10u1/4f7LimTvr7+I+vdu3ePatWqlXheQlQWcutSiOeAj48P+/bt49ixYwDk5eVx9OjRx9Zzd3fn1KlTGs1y8uRJPDw8NHoMIcqTFDohngPW1tYsX76csLAwmjdvjpeXV5Gv+fv5+XH06FGNzRh///59Tpw4gY+Pj0b2L4QmSPcCIXTMuHHjaNWqFaGhoeW+7+XLl5OcnEx4eHi571sITZEWnRA65qOPPtJYi66wsJAPP/xQI/sWQlOkRSeEEEKnSYtOCCGETpNCJ4QQQqdJoRNCCKHTpNAJIYTQaVLohBBC6DQpdEIIIXTa/wO98G8VKrsrcAAAAABJRU5ErkJggg==\n"
          },
          "metadata": {}
        }
      ],
      "source": [
        "results.plot(label='coffee')\n",
        "\n",
        "decorate(xlabel='Time (min)',\n",
        "         ylabel='Temperature (C)',\n",
        "         title='Coffee Cooling')"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "absent-arkansas",
      "metadata": {
        "id": "absent-arkansas"
      },
      "source": [
        "The temperature after 30 minutes is 72.3 °C, which is a little higher than the measurement we're trying to match, which is 70 °C."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "id": "absolute-desire",
      "metadata": {
        "id": "absolute-desire",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "c4dd65fe-21ca-4cf9-f1ea-993b19bbd535"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "72.2996253904031"
            ]
          },
          "metadata": {},
          "execution_count": 12
        }
      ],
      "source": [
        "coffee.T_final"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "light-carpet",
      "metadata": {
        "id": "light-carpet"
      },
      "source": [
        "To find the value of `r` where the final temperature is precisely 70 °C, we could proceed by trial and error, but it is more efficient to use a root-finding algorithm."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "funny-weekly",
      "metadata": {
        "id": "funny-weekly"
      },
      "source": [
        "## Finding Roots\n",
        "\n",
        "The ModSim library provides a function called `root_scalar` that finds the roots of non-linear equations. As an example, suppose you want to find the roots of the polynomial\n",
        "\n",
        "$$f(x) = (x - 1)(x - 2)(x - 3)$$\n",
        "\n",
        "A *root* is a value of $x$ that makes $f(x)=0$. Because of the way I wrote this polynomial, we can see that if $x=1$, the first factor is 0; if $x=2$, the second factor is 0; and if $x=3$, the third factor is 0, so those are the roots.\n",
        "\n",
        "I'll use this example to demonstrate `root_scalar`. First, we have to\n",
        "write a function that evaluates $f$:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "id": "small-shark",
      "metadata": {
        "id": "small-shark"
      },
      "outputs": [],
      "source": [
        "def func(x):\n",
        "    return (x-1) * (x-2) * (x-3)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "cleared-cylinder",
      "metadata": {
        "id": "cleared-cylinder"
      },
      "source": [
        "Now we call `root_scalar` like this:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "id": "defensive-content",
      "metadata": {
        "id": "defensive-content",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "ecdb62dd-4974-4156-bac4-81a067a00228"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "      converged: True\n",
              "           flag: converged\n",
              " function_calls: 3\n",
              "     iterations: 2\n",
              "           root: 2.0"
            ]
          },
          "metadata": {},
          "execution_count": 14
        }
      ],
      "source": [
        "res = root_scalar(func, bracket=[1.5, 2.5])\n",
        "res"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "adult-management",
      "metadata": {
        "id": "adult-management"
      },
      "source": [
        "The first argument is the function whose roots we want. The second\n",
        "argument is an interval that contains or *brackets* a root. The result is an object that contains several variables, including the Boolean value `converged`, which is `True` if the  search converged successfully on a root, and `root`, which is the root that was found."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "id": "behind-perth",
      "metadata": {
        "id": "behind-perth",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "e2713f1e-10f7-45b5-d112-57658ce95fe2"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "2.0"
            ]
          },
          "metadata": {},
          "execution_count": 15
        }
      ],
      "source": [
        "res.root"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "disturbed-shade",
      "metadata": {
        "id": "disturbed-shade"
      },
      "source": [
        "If we provide a different interval, we find a different root."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 16,
      "id": "following-sentence",
      "metadata": {
        "id": "following-sentence",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "fec835f0-3659-454b-cec4-a8f3523f57b0"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "2.9999771663211003"
            ]
          },
          "metadata": {},
          "execution_count": 16
        }
      ],
      "source": [
        "res = root_scalar(func, bracket=[2.5, 3.5])\n",
        "res.root"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "eligible-updating",
      "metadata": {
        "id": "eligible-updating"
      },
      "source": [
        "If the interval doesn't contain a root, you'll get a `ValueError` and a message like `f(a) and f(b) must have different signs`.\n",
        "\n",
        "Now we can use `root_scalar` to estimate `r`."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "commercial-correlation",
      "metadata": {
        "id": "commercial-correlation"
      },
      "source": [
        "## Estimating r\n",
        "\n",
        "What we want is the value of `r` that yields a final temperature of\n",
        "70 °C. To use `root_scalar`, we need a function that takes `r` as a parameter and returns the difference between the final temperature and the goal:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 17,
      "id": "genetic-compound",
      "metadata": {
        "id": "genetic-compound"
      },
      "outputs": [],
      "source": [
        "def error_func(r, system):\n",
        "    system.r = r\n",
        "    results = run_simulation(system, change_func)\n",
        "    return system.T_final - 70"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "happy-bridal",
      "metadata": {
        "id": "happy-bridal"
      },
      "source": [
        "This is called an *error function* because it returns the\n",
        "difference between what we got and what we wanted, that is, the error.\n",
        "With the right value of `r`, the error is 0.\n",
        "\n",
        "We can test `error_func` like this, using the initial guess `r=0.01`:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 18,
      "id": "french-financing",
      "metadata": {
        "id": "french-financing",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "3597f06d-cb95-4460-b3f5-b86e6fadd46f"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "2.2996253904030937"
            ]
          },
          "metadata": {},
          "execution_count": 18
        }
      ],
      "source": [
        "coffee = make_system(T_init=90, volume=300, r=0.01, t_end=30)\n",
        "error_func(0.01, coffee)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "statutory-activation",
      "metadata": {
        "id": "statutory-activation"
      },
      "source": [
        "The result is an error of 2.3 °C, which means the final temperature with `r=0.01` is too high."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 19,
      "id": "excellent-fellow",
      "metadata": {
        "id": "excellent-fellow",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "4f36ae0c-7886-4d0a-b35a-aa9af43e4c84"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "-10.907066281994297"
            ]
          },
          "metadata": {},
          "execution_count": 19
        }
      ],
      "source": [
        "error_func(0.02, coffee)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "unlimited-living",
      "metadata": {
        "id": "unlimited-living"
      },
      "source": [
        "With `r=0.02`, the error is  about -11°C, which means that the final temperature is too low. So we know that the correct value must be in between.\n",
        "\n",
        "Now we can call `root_scalar` like this:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 20,
      "id": "french-decline",
      "metadata": {
        "id": "french-decline",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "8f17315d-4c62-4f8b-e9f2-019d850b2af5"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "      converged: True\n",
              "           flag: converged\n",
              " function_calls: 6\n",
              "     iterations: 5\n",
              "           root: 0.011543084190599507"
            ]
          },
          "metadata": {},
          "execution_count": 20
        }
      ],
      "source": [
        "res = root_scalar(error_func, coffee, bracket=[0.01, 0.02])\n",
        "res"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "tamil-absorption",
      "metadata": {
        "id": "tamil-absorption"
      },
      "source": [
        "The first argument is the error function.\n",
        "The second argument is the `System` object, which `root_scalar` passes as an argument to `error_func`.\n",
        "The third argument is an interval that brackets the root.\n",
        "\n",
        "Here's the root we found."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 21,
      "id": "afraid-ordering",
      "metadata": {
        "id": "afraid-ordering",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "25cb32f3-9216-4764-dabf-53fb39dde750"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "0.011543084190599507"
            ]
          },
          "metadata": {},
          "execution_count": 21
        }
      ],
      "source": [
        "r_coffee = res.root\n",
        "r_coffee"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "partial-definition",
      "metadata": {
        "id": "partial-definition"
      },
      "source": [
        "In this example, `r_coffee` turns out to be about `0.0115`, in units of min$^{-1}$ (inverse minutes).\n",
        "We can confirm that this value is correct by setting `r` to the root we found and running the simulation."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 22,
      "id": "derived-annual",
      "metadata": {
        "id": "derived-annual",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "f69aadd1-903e-4df9-c620-9a51025bf18d"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "70.00000057308064"
            ]
          },
          "metadata": {},
          "execution_count": 22
        }
      ],
      "source": [
        "coffee.r = res.root\n",
        "run_simulation(coffee, change_func)\n",
        "coffee.T_final"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "prerequisite-conditioning",
      "metadata": {
        "id": "prerequisite-conditioning"
      },
      "source": [
        "The final temperature is very close to 70 °C."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "outer-powder",
      "metadata": {
        "id": "outer-powder"
      },
      "source": [
        "## Summary\n",
        "\n",
        "This chapter presents the basics of heat, temperature, and Newton's law of cooling, which is a model that is accurate when most heat transfer is by conduction and convection, not radiation.\n",
        "\n",
        "To simulate a hot cup of coffee, we wrote Newton's law as a difference equation, then wrote a version of `run_simulation` that implements it. Then we used `root_scalar` to find the value of `r` that matches the measurement from my hypothetical experiment.\n",
        "\n",
        "All that is the first step toward solving the coffee cooling problem I posed at the beginning of the chapter. As an exercise, you'll do the next step, which is simulating the milk. In the next chapter, we'll model the mixing process and finish off the problem."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "smart-yeast",
      "metadata": {
        "id": "smart-yeast"
      },
      "source": [
        "## Exercises\n",
        "\n",
        "This chapter is available as a Jupyter notebook where you can read the text, run the code, and work on the exercises.\n",
        "You can access the notebooks at <https://allendowney.github.io/ModSimPy/>."
      ]
    },
    {
      "cell_type": "markdown",
      "id": "scientific-attachment",
      "metadata": {
        "id": "scientific-attachment"
      },
      "source": [
        "### Exercise 1\n",
        "\n",
        "Simulate the temperature of 50 mL of milk with a starting temperature of 5 °C, in a vessel with `r=0.1`, for 15 minutes, and plot the results. Use `make_system` to make a `System` object that represents the milk, and use `run_simulation` to simulate it.\n",
        "By trial and error, find a value for `r` that makes the final temperature close to 20 °C."
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "from numpy import linspace\n",
        "import matplotlib.pyplot as plt"
      ],
      "metadata": {
        "id": "jJJaC9W6Qene"
      },
      "id": "jJJaC9W6Qene",
      "execution_count": 23,
      "outputs": []
    },
    {
      "cell_type": "code",
      "execution_count": 27,
      "id": "committed-angel",
      "metadata": {
        "id": "committed-angel"
      },
      "outputs": [],
      "source": [
        "# Solution goes here\n",
        "# Define the System class and necessary imports\n",
        "class System:\n",
        "    def __init__(self, T_init, r, volume, T_final=None):\n",
        "        self.T_init = T_init\n",
        "        self.r = r\n",
        "        self.volume = volume\n",
        "        self.T_final = T_final\n",
        "\n",
        "# Define the simulation and change functions\n",
        "def run_simulation(system, change_func):\n",
        "    \"\"\"Runs a simulation of the system.\n",
        "\n",
        "    Adds TimeSeries to `system` with the results.\n",
        "\n",
        "    system: System object\n",
        "    change_func: function that computes the change in temperature\n",
        "    \"\"\"\n",
        "    timesteps = linspace(0, 15, 100)  # simulate for 15 minutes\n",
        "    temperatures = []  # to store temperatures at each timestep\n",
        "    T = system.T_init  # initial temperature\n",
        "    for t in timesteps:\n",
        "        dT = change_func(T, system)  # compute the temperature change\n",
        "        T += dT  # update the temperature\n",
        "        temperatures.append(T)  # record the temperature\n",
        "    system.T_final = T  # store the final temperature for later use\n",
        "    return timesteps, temperatures  # return the time and temperature arrays for plotting\n",
        "\n",
        "def change_func(T, system):\n",
        "    \"\"\"Computes the change in temperature using Newton's law of cooling.\n",
        "\n",
        "    T: current temperature\n",
        "    system: System object containing parameters of the simulation\n",
        "    \"\"\"\n",
        "    return -system.r * (T - 20)  # environment temperature is assumed to be 20°C\n"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Create the system object for the milk\n",
        "milk = System(T_init=5, r=0.1, volume=50)"
      ],
      "metadata": {
        "id": "OTEFlog0Q2O3"
      },
      "id": "OTEFlog0Q2O3",
      "execution_count": 28,
      "outputs": []
    },
    {
      "cell_type": "code",
      "execution_count": 29,
      "id": "military-military",
      "metadata": {
        "id": "military-military"
      },
      "outputs": [],
      "source": [
        "# Solution goes here\n",
        "# Run the simulation and capture the time and temperatures for plotting\n",
        "timesteps, temperatures = run_simulation(milk, change_func)"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Print the final temperature\n",
        "print(\"Final temperature:\", milk.T_final)\n",
        "\n",
        "# Plotting results\n",
        "plt.figure(figsize=(8, 4))\n",
        "plt.plot(timesteps, temperatures, label='Temperature of Milk')\n",
        "plt.xlabel('Time (minutes)')\n",
        "plt.ylabel('Temperature (°C)')\n",
        "plt.title('Cooling of Milk Over Time')\n",
        "plt.legend()\n",
        "plt.grid(True)\n",
        "plt.show()"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 328
        },
        "id": "NArvRoVjQknz",
        "outputId": "a8e96557-5523-48cf-f3e7-21488fc8a1c9"
      },
      "id": "NArvRoVjQknz",
      "execution_count": 30,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Final temperature: 19.999601579016687\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 600x300 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgIAAAElCAYAAACWBIrfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAuJAAALiQE3ycutAABCyElEQVR4nO3deVhU9f4H8Pcwww4C4oYMgoiooCxK7vualWamN6/ibmqheVtuadnu/bl0s+4tNM0Ud7OrYWqLilkupeGuJLghIIKCgizDwMx8f3+gkyirznCYM+/Xow9zljnn/QEezmfOqhBCCBAREZFVspE6ABEREUmHjQAREZEVYyNARERkxdgIEBERWTE2AkRERFaMjQAREZEVYyNAVEv8/PwQGxsLAFi/fj26du0qWZYvv/wSXl5ecHFxwfHjxx96OSkpKXBxcUFubi4AYMKECfjHP/4BAEhOToZCoUBOTo4JEtcNUv/ciMyBjQDRPQ4cOIDBgwfDw8MD7u7uCA0NxaJFi1BcXGzS9YwZMwaHDh0y6TKrq6SkBC+99BI2b96M/Px8hIeHPzBP7969oVAosGfPnjLjP/roIygUCuPGvlmzZsjPz4ebm5tJsu3Zswc9evSAi4sL3NzcMHjwYBw7dswky66O4OBguLi4wMXFBba2trCzszMOBwcHS/pzIzIXNgJEd+zYsQODBw/GoEGDcP78eeTk5ODrr79GQkICrl27JnU8k8nIyEBRURHatWtX6XytWrXCqlWryoxbtWoVWrdubZZc3333HYYNG4Zx48bh2rVrSE5ORs+ePdGzZ0/Ex8ebfH1CCOj1+jLjzp49i/z8fOTn52PMmDF48cUXjcNnz541eQaiuoCNABFKNwovvfQS3njjDfzjH/9AgwYNAACtW7dGTEwMfH19AQDx8fHo1q0b3N3dERQUhI0bN5ZZxscff4wWLVqgfv36ePzxx3Hp0qVy1xcTE4OwsDDjsJ+fHxYtWoTOnTvD1dUVvXr1QmpqqnH62bNnjdP69OmD119/Hb17966wngsXLmDQoEGoX78+WrRogU8//RQAcPz4ceOGXK1Wo0WLFhUuY9SoUfjhhx+Mu/0PHz4MAOjUqZNxnprs/v/999/h7e2NrVu3PjBNCIFZs2Zh9uzZeP755+Hq6goPDw/MmTMHzz33HF577TUAwKxZszBp0qQy7124cCEGDRpkXM5///tftG7dGu7u7ujduzf+/PNP47x+fn6YP38+OnfuDCcnJyQkJFSZ+17l/dzmz5+Pxx57DM7Ozhg8eDBu3ryJF198Ee7u7mjZsmWZPQglJSV455130KJFC3h6emLo0KFIT0+vUQYiU2MjQATg/PnzuHz5Mv7+979XOE9OTg4ef/xxjBo1Cjdu3MDSpUvx/PPP4+DBgwCAtWvXYvHixYiNjUV6ejqCg4MxZMgQ6HS6amVYt24dNm7ciBs3bsDZ2Rlvv/02gNKNx9ChQzF48GBkZ2djwYIFWLlyZYXL0el0eOqppxAaGor09HR8++23WLRoETZs2IDw8HDjJ9u0tDRcvHixwuW4u7vj8ccfNzY7K1euxMSJE6tVy/127tyJ4cOHY8OGDRg+fPgD05OSkpCcnIzRo0c/MG306NE4cOAANBoNxo4diy1btkCj0Rinr127FuPGjQMALF26FF999RW2b9+OrKwsDB8+HEOGDClzaCcmJgarV69Gfn4+WrVq9VD13Ovrr7/G1q1bkZ6ejtTUVHTu3Bn9+/dHdnY2Ro8ejenTpxvnfeutt3Dw4EEcOHAA165dQ2BgIEaNGvXIGYgeBRsBIgA3btwAAHh7e1c4z86dO9GwYUPMnDkTtra26NWrF0aPHo3Vq1cDKN0gvfTSS2jXrh0cHBzwf//3f0hNTcWRI0eqleHFF19E8+bN4eDggDFjxuDo0aMASj9JZ2dn46233oKdnR06deqE5557rsLlHD58GNeuXcO8efPg4OCAkJAQzJgxAzExMdX8bvxl4sSJWLVqFTQaDbZs2YKxY8fWeBmrVq3CCy+8gB9//BG9evUqd56srCwAQNOmTR+Y1rRpU+j1ety8eRMRERFQq9XYtm0bgNI9HCkpKXjmmWcAANHR0fjggw/QsmVLqFQqvPTSS9BoNMa9GQDwwgsvoFWrVlAqlbCzs6txPfd74YUX4OPjAzc3NzzxxBPw9PTE8OHDoVQq8dxzz+HMmTMoLi6GEAJLlizB4sWL4eXlBTs7O8ybNw8HDx4ss/eHqLaxESACjIcCrl69WuE8aWlp8PPzKzPO398faWlp5U63t7dH06ZNjdOr0qRJE+NrZ2dn5OXlAQDS09Ph5eUFlUplnN6sWbNKczZt2rTMRu7enDXRr18/XLt2DR9++CG6dOlSJmN1LVy4EJGRkQgJCalwnrvf//J2k6enp0OpVKJ+/foAgLFjx2LNmjUAgDVr1uDZZ5+Fk5MTgNJDFZGRkXB3dzf+v3XrVpnaK/vePYzGjRsbXzs5OT0wLIRAYWEhsrKyUFBQgJ49exqzNWnSBHZ2dmwESFJsBIgABAYGws/PD5s2bapwHrVajeTk5DLjkpOToVary51eXFyM9PR04/SH1bRpU2RkZJQ5xJCSklJpzvT0dJSUlJSbsyZsbGwwfvx4LFiw4KEPC/zwww/YsGEDPvroowrnCQwMhK+vb5lzLu7auHEjunXrBkdHRwClV1zExcXh6tWr2LhxY5m9FD4+Pvjmm2+Qk5Nj/F9YWFjmkI+NjTR/9jw9PeHk5ITDhw+XyafRaHhJIkmKjQARAIVCgc8++wwLFizAZ599huzsbAClx64nT56MK1eu4IknnsD169exZMkS6HQ67N+/H+vXrzcen46MjMTnn3+OhIQEaLVazJ07F97e3ujYseMjZevcuTPc3d0xf/58lJSU4I8//sDmzZsrnL9jx45o3Lgx3nnnHWi1Wpw5cwafffYZxo8f/1Drf/nll7Fr1y4MGTLkod7fvHlz/PLLL1i6dCnmz59f7jwKhQKffPIJ5s+fj6+++gr5+fnIycnBwoULsWnTJixatMg4r4+PD7p3747JkyfDzs4Offr0MU6LiorCO++8g8TERADA7du3sW3bNuPeFSnZ2Nhg+vTpePXVV417ALKzs/H1119LnIysHRsBojueeuop/PDDD9i5cydatGgBd3d3jBgxAq1bt4aXlxc8PDzwww8/YN26dfD09MTUqVOxdOlSdO/eHQAwbtw4zJw5E0899RSaNGmCkydPYvv27WV26T8MW1tbbNu2DTt27ICHhwdef/11REZGwt7evsL5d+zYgaNHj6JJkyYYOnQoXnnllXJPxKuO+vXro3///rC1tX3oGnx9ffHLL7/gq6++wocffljuPM888wy2bNmCVatWoUmTJmjWrBn27t2Ln3/+ucyVCkDp9/qnn35CZGQkFAqFcfyMGTMwYcIEDB8+HPXq1UObNm2wYcOGh85tavPnz0eXLl3Qt29fuLq6okOHDti1a5fUscjKKYQQQuoQRFQz06ZNg8FgwJdffil1FCKycNwjQGQB9u/fj9TUVBgMBsTFxWH9+vUYOXKk1LGISAYebZ8lEdWKS5cuYdSoUbh16xbUajUWLFiAgQMHSh2LiGSAhwaIiIisGA8NEBERWTE2AkRERFbMYs4RqFev3iPfmOV+hYWFxjuSyZXca2R9lk/uNcq9PkD+NcqhvrS0NNy+fbvcaRbTCKjV6ho/KawqcXFx6Nevn0mXWdfIvUbWZ/nkXqPc6wPkX6Mc6gsKCqpwGg8NEBERWTE2AkRERFaMjQAREZEVYyNARERkxdgIEBERWTGzNgKpqano168f2rRpg+DgYMyZM8c4bfbs2QgICEBgYCC2bNlizhhERERUAbNePqhSqbBw4UJERESguLgY/fr1w7Zt2+Ds7IxDhw4hMTERGRkZ6NKlCwYNGgQXFxdzxiEiIqL7mLUR8PLygpeXFwDAzs4O4eHhSElJQUJCAiZMmAClUglvb29069YNu3btwvDhw80Zh4jIpIQQEAIQAAx3XhvuPL7l7muNTiCvqAQGAUAAAn+9Rwhx52vp+Dv/jO9/YF5Rdtpfr//KY8x2Zxrume+voXuH75t+z9NnBO4ZuG/avcPJtw04nZb7wPzlvueB6eW8p4plPKjqR+ZUtYzKJp+/ZUC95Js1XubDcrC1QYja3TwLL0et3VDo5s2biI2Nxa5dux7Y6Ddr1gypqam1FYWITEBvENDq9NCWGFCsN6BYZ4BWp0exTqBEXzquRHfnq750XMk9r3V6A4r1AueSdTj38wWU6A3QGwRK9AJ6gwE6g4BOL6Az/DWsN9wZ1gvoRemw3iBgEH9NMxhKp939qjcAhrvz3LOx1hvKvjaI0o2SQZS+vnfDfvd1mWHUYEPw8y5z/ijqhsMHpE5gXvG/1dqqmjdwxs+v9a619dVKI1BcXIwRI0Zg1qxZaN26dbkd4P2io6MRHR1tHM7MzERcXJxJcxUUFJh8mXWN3GtkfRXTGwSK9IBGBxTpBYrufNXqAK2+9HWxHtDqBbT60nHFeoFiA1Csv/PfUDpPiaF0Womh9LXOAOhN+WnofGK5o20Upf+Vd/7b3PO19L/intelJz3dfa0oM04BlQKwuzOsUAAK5V/zKe59H/4ap1AoSk+kuvv1zrLuTr+bUQEA97/3zuuSkhLY29oapwMK4zTcM9/d1yhv+N7l3//6vvfc6975cd/85a3v/vHlvUlhrOIvWq0W9vb2leaoTLnTK8tU3WXUcKaKJhcVFcHBwaE6azAJB2Vxrf5dM3sjoNfrMXr0aISFheHVV18FAPj4+JTZA5CSkoKOHTuWeV9UVBSioqKMw0FBQSa/xaMcbhtZFbnXKPf6ftq9B+0e64ZbhSW4VViMnMIS5BQWI0dTgtx7/t/WlOB2kQ55RSXIu/O1qMRQo3XZKhVwsFXC0VYJRzslHO2VcLVVwkFlAwdbJRxsbWCv+uurvcoGdqrS13Z3XtupbGCv/Ou1rdIGtkoF7JQ2sL0zrLJRwE5V+tVWaYPffzuE3j17wFapgPLOOJVN6WtFeVsVCyP331FA/jXKvT6zNwJTp06Fq6srPv74Y+O44cOH48MPP8T48eORkZGBAwcOYPny5eaOQlQnlOgNuJGnxfU8LTJvF+F6nhZZeVpk5d/9X4ybBcXIztfidpEOqOSTga1SATdHW9RztEU9B1t4uzuinoMtXB1UcLFXweXuV3sVnOxVcLFXwtlOBWd7FRztSl872inhZKeErVKaq4mTHBRo6GovybqJyMyNwMGDB7Fy5Uq0bdsW4eHhAIBJkybhpZdewu7duxEYGAgbGxssXrwYrq6u5oxCVCuEEMguKEbaLQ2u3tLgak4h0nOKcC1Xg2u5RbiWW4SsfG25x5ZtlQp4OtujgasdmtV3QpiPO/KzMxDWJgAeTnbwcLKFu5Md3J1s4e5kCzdHWzjaKmXxqZmIpGPWRqBbt24Vng+waNEiLFq0yJyrJzILIQQybhfh0o0CXMoqwJWsAly5WYiU7EKk3CyEpkRfZn6FAmjkag8vN0dE+HqgiZsDGtdzQON69mjk6oBGrvZo6GoPN0fbBzbqcXFx6NcnoDbLIyIrYzGPISaqbQaDQOqtQiRl5iMpMw9JmXm4cD0fl24UPLCxb1zPHr71nfFkiBd8PJyg9nCEt4cj1B6OaFzPQbLd7kREVWEjQARAq9Pj3LU8nE2/jYRruUhIv41zGXkoLP5rg+9ir0JAIxc80c4L/g2d0aKhM5o3cIGvpxMcbJUSpicienhsBMjqCCFwOasAx1JycDI1ByfTcvDntdsouXM9nMpGgYBGLni8bRO0buKKwMal/73cHHg8nohkh40AyZ5Ob8DZ9Ns4fDkb8cm3cPTKLWQXFAMo3ei3auKKER18EKp2Q1tvN7Rs7AJ7FT/hE5F1YCNAsiOEQFJmPg5eyMKhi9k4fDkbeUU6AICHky0i/OojwtcDEX4eCG7qxt36RGTV2AiQLOQVleDghWz8knQd+xJv4FpuEQCgvrMderZsiC4tPNHZ3xMtGjpz9z4R0T3YCJDFupGnxS9pOqxZeQSHLmahRC9gq1Qgwrc+JnT1Q4+WDdG6iStsbLjhJyKqCBsBsihZ+VrsPHUNO06lI/7KLQgB1HfOxfBwNfq2aYRuAQ3gYs9fayKi6uJfTKrzCrQ6/HQ2A7En0nHwQhb0BoGmbg6Y0NUPDbXpmDa8P5T81E9E9FDYCFCdJITAidQcfP1HKrafTEdBsR7uTrYY9ZgPhoV7o0MzD9jYKBAXl8EmgIjoEbARoDolX6vDlqNp2HgkBecy8qC0UaBf60b4W4QPegY2hJ2Kd+gjIjIlNgJUJ6RkF2L1b8nY/Ecq8rQ6+Hk64Y3HW+PZDt5o5Fp7zwEnIrI2bARIUsdTbmHpvovY/WcmhAB6BTbEpO7N0SOgAc/2JyKqBWwESBK/X8rG53sv4MCFLDjY2mBMp2aY0LU5Ahq5SB2NiMiqsBGgWvXbxWx8sjsJR5JvwsVehRd7t8Ck7s3RwMVe6mhERFaJjQDVioT021j44zn8knQD9RxU+Ef/lpjQ1Q/uTnZSRyMismpsBMisUm8WYvHuJMSeuAp7lQ1e6N0C03u1gJujrdTRiIgIbATITIpK9Pjil4tYuu8idAaBUY/5YFa/QDRx4xUARER1CRsBMrm4PzPx/vYEpNwsRPeABnhvaBACGrlKHYuIiMrBRoBM5lquBm/HnsGeP6+jqZsDloxpj8Ftm/Bpf0REdRgbAXpkQgh8E5+GD3ckQFOix4u9W2BG3wA42fHXi4ioruNfanok6TkazNl6Gr8k3UAbr3r498gQBDd1kzoWERFVk1lv3D5r1iyo1WqoVH/1G1qtFuPGjUO7du3Qpk0bzJgxAwaDwZwxyEy2nbiKQZ/8ioMXsvBy/0B8N6MbmwAiIgtj1kZg5MiRiI+PLzNuzZo10Ol0OH36NM6cOYMTJ07ghx9+MGcMMrHCYh1e/99JzNp0At4ejvhuRnfM6t8Stko+EIiIyNKY9dBA9+7dyx1fUFAAnU6H4uJiaLVaNGnSxJwxyIQSM/IQteEYLlzPx5hOzfD2U0FwsFVKHYuIiB5SrX+EGzduHOzs7NC4cWM0btwYAwcORIcOHWo7Bj2Eb+JTMfTzA8jMLUL06Pb41zPt2AQQEVk4hRBCmHslKpUKOp0OALB//34sXrwYmzZtgkajQe/evbFkyRJ07dq1zHuio6MRHR1tHM7MzMTmzZtNmqugoADOzs4mXWZdY4oa9QaBb87rsDtFD796CkwPsUVDx7pxGEDuP0O51wfIv0a51wfIv0Y51Ddz5kwkJCSUP1HUAqVSaXz94osviuXLlxuH//nPf4qPP/64ymW0adPG5Ln27Nlj8mXWNY9aY05BsYhc8bvwfWOHmLXxmNAU60yUzDTk/jOUe31CyL9GudcnhPxrlEN9lW1Da/1jna+vL+Li4mAwGKDVarF//34EBQXVdgyqhgvX8/B09AEcuJCF2YNb45PnwngogIhIZszaCEybNg1qtRp6vR5qtRpRUVGIiooCALRt2xZhYWHo3bs3Hn/8cXPGoIfwR/JNDF9yCFn5xfhqfASm92rBOwQSEcmQWa8aWLZsWbnjN23aZM7V0iPak5CJqA3H4OFkh9WTOqJVEz4ngIhIrnhnQSpjc3wq5mw9DV9PJ6yd3Ane7o5SRyIiIjNiI0BGX/xyEQt+OIdQtRtWTeyI+s52UkciIiIzYyNAAID/xp3H4t1J6NGyAb6I7ABne/5qEBFZA/61J0T/fAGLdyehb+tGWBrZHvYqXhlARGQt6sZdYUgyS/ddxEc/JaJ3q4ZsAoiIrBAbASv25a+XsPDHc+gZ2BBfRHZgE0BEZIXYCFiptb9fwb++/xM9WjbA8rEdeKMgIiIrxUbACv10NgPvbDuDx/w8sHxsBJsAIiIrxkbAyhy9chMvbTyOFg1d8OW4CDjasQkgIrJmbASsyMUb+Zi8Oh5ujraImfgY3J14nwAiImvHywetxPXbRRi/8gh0eoENUzpC7eEkdSQiIqoDuEfAChSV6PH8mnhk3i7CsrEdENS0ntSRiIiojuAeAZkTQuDt2DM4mZaLBcPboVtAA6kjERFRHcJGQOZ+TtPjm3NpGNOpGUZ1bCZ1HCIiqmN4aEDGjly+iU2JOkT4euDdIcFSxyEiojqIjYBMXcvV4MX1R+FqByyJbA87FX/URET0IG4dZKhYZ8D0dcdwW6NDVKgdGrk6SB2JiIjqKDYCMrR4dxJOpubg3aFB8Hfjj5iIiCrGrYTMHLyQhWW/XsTjwU0wmicHEhFRFdgIyMjNgmK8/PUJNHZ1wIJn20GhUEgdiYiI6jhePigTQgi8/r9TuJGvxYYpnXn7YCIiqhbuEZCJdb9fwZ4/MxHVOwBdWnhKHYeIiCyEWRuBWbNmQa1WQ6Uqu+Ph8OHD6NixI4KDgxEcHIz09HRzxpC9C9fzMW/nnwjzcces/i2ljkNERBbErI3AyJEjER8fX2ZcXl4exo0bh9WrV+Ps2bM4dOgQPDw8zBlD1gwGgTe2nIJCAXz6XBhsldzJQ0RE1WfWcwS6d+/+wLgNGzbgiSeeQJs2bQAAbm5u5owge2t/v4KjV25h7pNt4NfAWeo4RERkYWr94+O5c+eg1WrRv39/hIeH491334UQorZjyELarUIs/PEcQn3cMbFbc6njEBGRBVKISrbC+/fvR0xMDPbv349r167B0dERbdu2xZAhQzBx4kS4u7tXayUqlQo6nQ4AMHPmTOzbtw/79++Hg4MDnn76aURGRmLs2LFl3hMdHY3o6GjjcGZmJjZv3vwQJVasoKAAzs6W+SlaCIFPj5fgz5sGvNPZDmqX8ns6S66xOlif5ZN7jXKvD5B/jXKob+bMmUhISCh/oqjA448/LsaNGydiY2PF1atXRUlJicjLyxNHjx4VixYtEl27dhU7duyo6O1lKJVK4+uFCxeKV155xTgcHR0tZs2aVeUy2rRpU6111cSePXtMvszasuVoqvB9Y4dYvCux0vksucbqYH2WT+41yr0+IeRfoxzqq2wbWuE5AsuXL4ePj0+ZcS4uLmjfvj3at2+P11577aHO9h8+fDgiIyOh1Wpha2uLuLg4PP744zVejjW7kafFBzsSENjYBS/2aSF1HCIismAVniPg4OBQ7m6EhIQEZGVlQaFQwNvbu9KFT5s2DWq1Gnq9Hmq1GlFRUQgICMDYsWMRHh6OkJAQNG7cGJMmTXr0SqzI/B/+RK6mBAueDYG9Sil1HCIismAV7hF45ZVXMGnSJAQFBZUZf/36dcyfPx9r166tcuHLli0rd3xUVBSioqJqGJUA4FjKLWw9dhV/79gM7ZvxsksiIno0Fe4RSExMRJ8+fR4Y37t3b5w4ccKcmagCBoPA+9+dhauDCq8NDJQ6DhERyUCFjUBubm6FbyouLjZLGKrc/46l4WRaLl7uHwhPF3up4xARkQxU2Ag0aNAAhw4demD8b7/9hvr165s1FD3odlEJFv14Di0buWBsF1+p4xARkUxUeI7AvHnzMGzYMEyfPh0REREAgPj4eCxbtgwbN26stYBU6rO488jKL8anz4XzNsJERGQyFTYCffr0wa+//ooFCxbg+++/BwC0bdsW+/btM94emGrHxRv5WHUwGYOCG6N7ywZSxyEiIhmp9FkDrVu3RkxMTC1FoYr8384/YWOjwNwng6qemYiIqAYqbAQ0Gg3Wr18PIQQiIyPh6OhYm7nojvjkm4g7dx3TevrDp76T1HGIiEhmKjzYPHbsWDRu3BheXl6IjIyszUx0hxACi35MhKu9CtN78Q6CRERkehU2AhkZGejcuTM6deqEa9eu1WYmuuOXpBs4knwTU3v6w8PZTuo4REQkQxUeGlixYgXmzJkDIQS++OKL2sxEKL150Ec/JcLT2Q4Tu/MRw0REZB4VNgKtW7fGihUrajML3ePHsxk4m34bbz8VBBf7Ss/pJCIiemgVHho4duxYpW8sKipCYmKiyQMRoNMb8O9diWjq5oAxnZpJHYeIiGSswo+a//nPf5Ceno7hw4cjIiICjRo1QlFREc6fP489e/bg0KFD+Pjjj9GqVavazGsVth6/iks3CrDw2XZwsOXTBYmIyHwqbARWr16NCxcuYO3atXjjjTeQnp4OR0dHtG3bFkOGDMFHH30EW1vb2sxqFYp1Bvxnz3k0b+CMZ9urpY5DREQyV+nB54CAALz//vu1lYUAxJ64iqs5Gnz6XBhUvJUwERGZGbc0dYjeIPDFLxfhU98RT4V4SR2HiIisABuBOmTX2QxculGAqT1bcG8AERHVCm5t6gghBJbsu4gGLvYY2YHnBhARUe2oViPw008/4fPPPwcAZGZmIikpyayhrNHBC9k4fTUXk7s355UCRERUa6psBObOnYslS5bgs88+AwAYDAZMmDDB3LmszpJ9F+DqoEJkZ943gIiIak+VjcCOHTsQGxsLJ6fSJ995eXkhPz/f7MGsyYnUHBy6mI1xXXzh6sBLMomIqPZU2Qg4ODhAoVAYh7VaLYQQZg1lbZb8fAH2KhtM7MZnChARUe2qshHo378/3n77bRQWFuK7777DsGHDMGLEiNrIZhUuXM/DroRMjHrMBw1c7KWOQ0REVqbKRuDDDz+Ev78/QkJCsGbNGowYMQLvvvtutRY+a9YsqNVqqFQP3rcoJycH3t7emDJlSs1Ty0jMoWQoFMCUHv5SRyEiIitU6Z0F9Xo9Bg4ciLi4OEycOLHGCx85ciTmzJkDtfrBy+Fef/119O3bt8bLlJPbRSXYeuwq+rVuDJ/6TlLHISIiK1TpHgGlUgmdTgeNRvNQC+/evTuaNGnywPiff/4ZWq0W/fr1e6jlysU38WkoLNZjQlc/qaMQEZGVUogqzvyLjIzEyZMnMXToUDg7OxvHv/nmm9VeiUqlgk6nAwBoNBr06dMH27dvx86dO3HgwAGsWLHigfdER0cjOjraOJyZmYnNmzdXe53VUVBQUKam2mQQAm8dLIbSBviwi12ZEzJNScoaawPrs3xyr1Hu9QHyr1EO9c2cORMJCQnlTqv00AAAtGzZEi1btgQAlJSUPHKY9957D1OnTkXDhg0rnS8qKgpRUVHG4aCgIJPvQYiLi5Nsr8TPiddxXfMHPnw6GP27+JltPVLWWBtYn+WTe41yrw+Qf41yr6/KRqC6JwZW16FDh/D111/jgw8+QH5+PrRaLWxsbLB8+XKTrqeuW30oGa72Kgzno4aJiEhCVTYCH3zwQbnj33nnnYda4f79+42vY2JicODAAatrAi5nFWBf4g1M6OoHZ/sqfwRERERmU+Xlg0ql0vi/pKQE27Ztw+XLl6u18GnTpkGtVkOv10OtVpfZ1W/N1vyWDAAY18VX2iBERGT1qvw4+tZbb5UZnjNnDoYMGVKthS9btqzS6RMmTLC65xYUaHX4X3waegU2hH9DF6njEBGRlavxY4iVSiVSUlLMkcUqbD1+FXlaHS8ZJCKiOqHKPQIDBgwwXtpmMBhw7tw5jB8/3uzB5OrrP1Kg9nBEr8DKr5ogIiKqDVU2AnPnzv1rZpUKfn5+8Pb2NmsouTqbnoszV2/jlQGBsLExz30DiIiIaqLKQwM//fQTevXqhV69eqFbt27w9vau0c2E6C/fxKdBoQCe7cBLBomIqG6oshH48ccfHxj3/fffmyWMnGl1esSeuIruAQ3g7e4odRwiIiIAlRwaWLt2LdasWYOLFy9i4MCBxvG3b99Go0aNaiWcnOxJuI6cwhL8LcJH6ihERERGFTYCPXv2hI+PD1555ZUylxC6uroiNDS0VsLJyeb4VNRzUGFAUGOpoxARERlV2Aj4+vrC19cXx44dq808spSeo8Gv529gbGdfONgqpY5DRERkVOVVA1euXMGcOXNw5swZFBUVGccnJSWZNZicbD2WBiHAwwJERFTnVHmy4MSJEzFy5Ejo9Xps3boVgwYNwujRo2sjmywYDAKb49PQxqsegpvWkzoOERFRGVU2Arm5uXjmmWegVCrRtm1bfPbZZ9i5c2dtZJOFI8k3kXKzEH+LUBtvzERERFRXVHloQKUqnaVBgwY4cuQImjZtiuzsbLMHk4vN8amwU9pgWBhvwkRERHVPlY3ACy+8gFu3buFf//oXRo0ahby8PCxatKg2slk8TbEeP57JQL82jeDhbCd1HCIiogdU2gjo9Xq4uLjAw8MDXbp0waVLl2orlyzsPXcdhcV6PB3WVOooRERE5ar0HAGlUon58+fXVhbZ2X4yHS72KvRuxRswERFR3VTlyYJPPPEEVq9ejeLi4trIIxt5RSXYm3gdA4Ma894BRERUZ1V5jsCCBQug1+sxceJEqFQqCCGgUCjYGFRhd0IminUGDAnlYQEiIqq7qmwESkpKaiOH7Gw/mQ53J1t0C2ggdRQiIqIKVXloACh9FPHnn38OALh+/TrvKliFWwXF2H8+C4PbNoGdqlrfYiIiIklUuZWaO3culixZgs8++wxA6ZUEEyZMMHcui/bDmQzoDAJDQnhYgIiI6rYqG4EdO3YgNjYWTk5OAAAvLy/k5+ebPZgl234yHQ1d7dHJ31PqKERERJWqshFwcHAoc2tcrVYLIYRZQ1my67eL8PvlbDzZzgtKG95SmIiI6rYqG4H+/fvj7bffRmFhIb777jsMGzYMI0aMqNbCZ82aBbVabbxNMQCsX78eoaGhCAkJQUREBPbu3fvw6eugnaevQQjwagEiIrIIVTYCH374Ifz9/RESEoI1a9ZgxIgRePfdd6u18JEjRyI+Pr7MuObNm2Pv3r04deoUVq1ahb///e8wGAwPl74O2n4yHd7ujmjfzF3qKERERFWq8vJBhUKBiRMnYtSoUQAAR0fHai+8e/fuD4zr2rWr8XXbtm2h1WpRUFAAV1fXai+3rkrP0eBYSg6m9fTnkwaJiMgiKEQVB/zPnDmDiRMn4vLlywAAf39/rFy5Em3btq32SlQqFXQ63QPjN27ciCVLlmD//v0PTIuOjkZ0dLRxODMzE5s3b672OqujoKAAzs7OJlteXIoOGxJ1eKujHfzd6sZlg6ausa5hfZZP7jXKvT5A/jXKob6ZM2ciISGh/ImiCh06dBDbtm0zDm/fvl106NChqreVoVQqHxh37Ngx0axZM5GUlFStZbRp06ZG66yOPXv2mHR5Y778XTw2b7fQ6w0mXe6jMHWNdQ3rs3xyr1Hu9Qkh/xrlUF9l29AqP7YWFxdj6NChxuGnnnrqke82mJSUhGeffRabNm1Cy5YtH2lZdUWupgS/X8pGvzaNYcOrBYiIyEJU2QgMGTIEy5Ytg0ajgUajwfLly8s0BjWVlpaGJ598EkuXLkWXLl0eejl1zb7E69AZBAYGN5Y6ChERUbVVebLgokWLoNfr8cILLxjHKZVKLFiwoMqHD02bNg07d+6EXq+HWq3G008/jeLiYly/fh1vvPEG3njjDQBAbGws/Pz8Hr0aCe1KyISznRJdW/AmQkREZDnM+tChZcuWlTv+yy+/fOhl1kVanR77zl1Hr1YNYa/iI4eJiMhyVNkIAKXnCWRkZJQ589/f399soSzNbxezUVCsx8CgJlJHISIiqpEqG4FPP/0UH374IRo2bGi8Nl6hUFR8GYIV2p2QCaWNAn1aNZI6ChERUY1U2Qh8/vnnOH/+POrXr18beSyOwSCwOyETnZrXh5uTrdRxiIiIaqTKqwb8/f3h4uJSG1ks0qmrubiep8XAIF4tQERElqdaVw3069cPPXv2hL29vXH8O++8Y9ZglmJ3QgYAoD8bASIiskBVNgKvvfYaGjVqBAcHB9jY1I3b5tYlu85mIsirHtQeTlJHISIiqrEqG4Fr165hz549tZHF4iRnFeD89XzM6iePuyMSEZH1qfIjfs+ePXH06NHayGJx4s5dBwAM4GEBIiKyUFXuEdi1axe+/PJLeHt7w97eHkIIKBQKJCUl1Ua+Ou2XpBto6GqP4Kb1pI5CRET0UKpsBPbu3VsbOSxOUYkehy9l48kQL+P9FYiIiCxNlYcGfH19odFocOLECfj6+sLT0xOurq61ka1OO3z5JrQ6A3oFNpQ6ChER0UOrshFYsmQJxo4di1dffRUAkJGRgREjRpg9WF33a9INKBRA94AGUkchIiJ6aFU2AsuWLcPBgweNewECAgJw48YNswer635NuoF23m7wdLGvemYiIqI6qsJG4ODBgwAAR0dH2NnZGccbDAbzp6rj0nM0OH89n4cFiIjI4lXYCMycORMAEB4ejhUrVkCn0+HUqVOYMmUK+vTpU2sB66Jfk0r3iPRkI0BERBauwkZACAEA+OSTT5CamgoHBwdMnjwZarUa//73v2stYF30S9INuDqoEO7jLnUUIiKiR1Lh5YOpqamYOnUqgNKmICwsDEDpnQZnzJiB5cuX10rAukanN+DAhSx0a9EAKiVvuUxERJatwkbA2dkZ3bp1q80sFuFEag7yinTo1YqHBYiIyPJV2Ah4enpi/PjxtZnFIvD8ACIikpMqzxGgsn45n4UWDZ3h7e4odRQiIqJHVmEjsH///trMYRFuFhTjVFoOegU2kjoKERGRSVTYCLi4uNRmDotw4EIWhAB6BvJugkREJA9mPe191qxZUKvVUKnKnoowe/ZsBAQEIDAwEFu2bDFnBJP6NekG7FQ26NTcU+ooREREJmHWRmDkyJGIj48vM27Pnj04dOgQEhMT8fPPP+Pll19Gfn6+OWOYzO+XshHh6wFHO6XUUYiIiEzCrI1A9+7d0aRJkzLjtmzZggkTJkCpVMLb2xvdunXDrl27zBnDJFJvFiLtlgad/bk3gIiI5KPCywfNJS0tDcOHDzcON2vWDKmpqQ/MFx0djejoaONwZmYm4uLiTJqloKCg2ss8kK4DANjeSkZcXIpJc5hTTWq0RKzP8sm9RrnXB8i/RrnXV+uNQHUvS4yKikJUVJRxOCgoCP369TNplri4uGov8/vNJ2GvSsekYX1hr7KcQwM1qdESsT7LJ/ca5V4fIP8a5V5frd8j18fHp8wegJSUFKjV6tqOUWO/X8pGhJ+HRTUBREREVan1RmD48OGIiYmBXq/H1atXceDAAQwcOLC2Y9RI6s1CXM3RoDOvFiAiIpkxayMwbdo0qNVq6PV6qNVqREVFYcCAAejcuTMCAwPRu3dvLF68GK6uruaM8ch+u5QNAOjcgo0AERHJi1nPEVi2bFm54xctWoRFixaZc9Um9fulbDjY2iBE7SZ1FCIiIpPic3SrIITA4Us3EeFbn+cHEBGR7LARqELaLU3p+QH+9aWOQkREZHJsBKrw28U75wfwRkJERCRDbASq8PulbDjaKhGidpc6ChERkcmxEaiEEMJ4/wA7Fb9VREQkP9y6VSL1pgbpuUU8LEBERLLFRqASv13KAgCeKEhERLLFRqASv1+6yfMDiIhI1tgIVOLwnfMDbJX8NhERkTxxC1eB9JzS8wMe8+NhASIiki82AhU4euUWAKCDr4fESYiIiMyHjUAFjqXcgo0CCPVxlzoKERGR2bARqMCxK7fQqkk9uNib9blMREREkmIjUA5NsR5n02+jg6+71FGIiIjMio1AOU6l5UBnEDw/gIiIZI+NQDmOptw5UbAZrxggIiJ5YyNQjmNXctDAxR4+9R2ljkJERGRWbATuI4TAsZRbaN/MHQqFQuo4REREZsVG4D7J2YW4WVDM8wOIiMgqsBG4D28kRERE1oSNwH2OXrkFW6UCbb3dpI5CRERkdmwE7nPsyi209XaDg61S6ihERERmJ1kjsG7dOoSEhCAsLAw9evRAYmKiVFGMbheVIOl6Hjo042EBIiKyDpI0AoWFhZg1axb27t2LEydOYMyYMZg7d64UUco4kZIDIYD2PD+AiIishCSNgMFggBAC+fn5AIDc3Fx4eXlJEaUMnihIRETWRiGEEFKseMOGDZg6dSrc3Nzg5uaG3377DW5uf52gFx0djejoaONwZmYmNm/ebNIMBQUFcHZ2Ng5/fLQYmYUGLOrhYNL1SOn+GuWG9Vk+udco9/oA+dcoh/pmzpyJhISE8icKCRQXF4tu3bqJhIQEIYQQb7/9tpg2bVql72nTpo3Jc+zZs8f4Wqc3iOB3fhQzNxwz+XqkdG+NcsT6LJ/ca5R7fULIv0Y51FfZNlSSQwMnTpyAEAJt2rQBAIwaNQqHDh2SIorR+et5yNfq0L6Zu6Q5iIiIapMkjYBarUZiYiKuXr0KANi9ezeCgoKkiGJ0Ki0XABDq4y5pDiIiotqkkmKlXl5eWLBgAQYMGABbW1s0bNgQK1eulCKK0em0XKhsFGjjVU/SHERERLVJkkYAAKZMmYIpU6ZItfoHnLqai8DGrryREBERWRXeWRBAsc6AP6/dRoiatxUmIiLrwkYAQFJmHop1BrRjI0BERFaGjQCA01dLTxQM8XaXNggREVEtYyOA0isG7JQ2CGziInUUIiKiWiXZyYJ1yam0HLT2coW9iicKElHdJKS5CWy11fV8j8pS6lMoFDV+j9U3AkUleiRm5OG5x3ykjkJE9ICSkhKkpqZCq9VKHaVC9evXx7lz56SOYTaWVJ+9vT18fHxga2tb7fdYfSNwLiMPOoNAqNpd6ihERA9ITU2Fq6sr/Pz8HurTXm24ffs26tWT7z1YLKU+IQSys7ORmpoKf3//ar/P6huB02k5AMArBoiozhFCQKvVws/PDzY2dfeULoVCUWebFFOwlPoUCgU8PT2RlZUFIUS1M9fd36xaciotF/YqG7RsxBMFiahusoSNENUND/O7YvWNwOmruQhuWg8qpdV/K4iIypWXl4ewsDCEhYWhWbNm8PT0NA7v379f6njlSk5Oxpo1ayRZt1arxaBBgxAWFobo6Ogy02JiYqBQKPDtt98ax50+fRoKhQLvvfceAOCLL77AihUrAADvvfce5s2b98BrU7LqQwNavUBSZh7GdfGTOgoRUZ3l6uqKEydOACjdkB04cMC4oQJKj6FLQafTQaUqfzN2txEYN26cyZZZXcePH0dBQYHxe3a/0NBQbNiwAc888wwAYP369QgNDTVOnz59+iOtv6as+mNwap6AQQDtvHl+ABFRTWg0GkyfPh0dO3ZEly5d8PnnnwMo3QD7+/tj+vTpCAoKQt++fXH8+HEMGDAA/v7++OSTT4zLUKlUePPNN9GuXTu0b98ep0+fBgAYDAa89dZb6NixI0JCQvDmm2+Wec/cuXMRHh6OrVu3YtWqVejYsSPCw8PRu3dvXL58GQDw8ssv4/DhwwgLC8Ps2bOxb98+9O/f37icefPmGT+Bv/feexgzZgx69uyJAQMGPLD+999/v9zvwdmzZ9GjRw+EhISgZ8+eSExMRGpqKiIjI3Hy5EmEhYXh8OHDD7wvIiIC58+fR15eHoQQ+P777/HUU08Zp1fnk//27dvRsWNHZGRkVDpfdVj1HoHk2wYA4DMGiMgivP6/k0jKzDf5cgMbu2DRiNCqZ7zH/Pnz0b59e3zxxRe4ceMGBg8ejL59+8LJyQnJycmYNGkSvvjiCwwfPhwzZszAnj17UFhYiJYtW+KFF16Ag4MD9Ho91Go1Tp8+jW3btmHixImIj49HTEwMAODIkSMwGAwYNmwYfvjhBwwePBh6vR7+/v44fvw4ACA7OxsTJ04EAGzduhVvvvkmNm7ciE8++QTz5s3Dnj17AAD79u2rtJ4TJ07g8OHDcHFxMT4N9+76n3rqKeP67xUZGYn3338fQ4cOxZYtWzBu3DgcPnwYK1asKLPu8gwbNgxbt25F8+bN0aFDBzg4OECn01Xre79u3TosW7YMP/30Ezw8PKr1nspYfSPgZKeEf0OeKEhEVBPff/89NBoNlixZAoPBgPz8fCQmJiI8PBze3t7o2LEjACAsLAwajQaOjo5wdHREgwYNkJGRAT8/PwDA+PHjAQBPP/00Jk6ciIKCAnz//fc4efIkdu7cCQAoKCjA+fPnjRvi0aNHG3MkJibirbfeQlZWFvR6/UNfXTF06FC4uLgYa7t3/Xl5eWXWDwC5ublIS0vD0KFDAQDPPvssnn/+eRQWFlZrfaNHj8bMmTPh5+eH0aNHl7vnoDzffPMNXF1d8eOPPxrzPiorbwQE2jZ1h9KGZ+QSUd1X00/t5iSEwPr16xEWFlbmOvvk5GTY29sb57OxsXlguKpPvkIIfPTRRxg2bNgD05RKJRwcHIzDY8aMwYYNG9ClSxecPn3aeNz9fiqVCgaDwThcVFRU5lwAZ2fnCtdf3n0EHvVKjsDAQOTm5uKXX37BkiVLqt0IBAUF4ejRozh//jzCw8MfKcNdVnuOQL5Wh4wCwfsHEBE9hMGDB+M///kP9Ho9AOD8+fMPddLg2rVrAQA7duyAv78/nJ2dMXjwYCxduhRFRUUAgPT09AqPhd++fRve3t4AgOXLlxvHu7q6lsnj5+eHc+fOQaPRoLCwEN9//32ltd27/mvXrj2w/nr16sHHxwc7duwAAHz77bcIDAyEk5NTtWtfsGABFi5cCKWy+re3Dw4OxjfffIPnnnsOp06dqvb7KmO1ewTOXs2FAM8PICJ6GHPnzsU///lPhIaGwmAwoFGjRvjf//5Xo2UolUqkpqYiJCQEKpXKeG7A5MmTkZaWhoiICCgUCjg7OyMmJgZNmjR5YBmLFi1Cr169UL9+fTz55JPG8aGhoXBzc0NISAieeOIJLFiwABMmTEC7du3g4+ODdu3aVZjr/vU7ODhg7dq1D6x/7dq1mD59Ot588024u7tj9erVNaq/d+/eNZr/rtDQUGzatAkjR47E1q1bERwc/FDLuUshLORJCkFBQUhISDDZ8r789RL+9f2f2PtqL1mfIxAXF4d+/fpJHcNsWJ/lk3uNj1KfEALnzp1D69at6/RNhR72FrwqlaraJ8hJyVJuMQxU/DtT2TbUavcIDG7XBBlXLsDP07nqmYmIiGTKas8RUHs4oWtTJWx4oiARkSQsYW+ANbDaRoCIiIjYCBAR1XkWcioX1QEP87si2TkCBQUFePHFF/H7779DoVDg5ZdfxrRp06SKQ0RU5ygUCtjb2yM7Oxuenp519oRBIYSsmxVLqU8IgezsbNjb29fod0WyRuDVV19FcHAwVq9eDSEEbty4IVUUIqI6y8fHB6mpqcjKypI6SoWKiorK3ORHbiypPnt7e/j4+NToPZI0Anl5efjuu++QkpICoLTrbdSokRRRiIjqNFtbW/j7+9fpT6R79+5F3759pY5hNpZU38PsNZLkPgInT57EhAkT0KlTJxw5cgS+vr749NNP4evra5wnOjq6zHOcMzMzsXnzZpPmKCgoKHNbSTmSe42sz/LJvUa51wfIv0Y51Ddz5syK78UjJBAfHy8AiN27dwshhPjqq69E3759K31PmzZtTJ5jz549Jl9mXSP3Glmf5ZN7jXKvTwj51yiH+irbhkpy1YBarYanp6fx2dCjRo3C0aNHpYhCRERk1SQ5R6Bx48YIDg7GsWPH0L59e+zevbvKeyWnpaUhKCjIpDlu3bplkmc512Vyr5H1WT651yj3+gD51yiH+tLS0iqcJtmzBhISEjB58mQUFBTA3d0dy5YtQ5s2bWo1g6mfX1AXyb1G1mf55F6j3OsD5F+j3OuT7PLBoKAg/Pbbb1KtnoiIiMA7CxIREVk1q24EoqKipI5gdnKvkfVZPrnXKPf6APnXKPf6JDtHgIiIiKRn1XsEiIiIrB0bASIiIitmtY3Avn37EBwcjICAAEyZMgV6vV7qSCaTmpqKfv36oU2bNggODsacOXOkjmQ2UVFRUKkku/jFrAoKCjB+/Hi0atUKrVu3xrJly6SOZFLr1q1DSEgIwsLC0KNHDyQmJkod6ZHNmjULarX6gd/J2bNnIyAgAIGBgdiyZYtE6R5defWtX78eoaGhCAkJQUREBPbu3SthwkdT0c8PAHJycuDt7Y0pU6ZIkMzMaunuhnWKXq8XLVq0EGfPnhVCCDFy5EgRExMjcSrTSU9PF3/88YcQQgitViu6d+8uYmNjJU5ler/++qsYN26cUCqVUkcxi2nTpomFCxcKIYQwGAwiMzNT4kSmU1BQIOrXry9u3LghhBBi6dKlYsSIERKnenT79+8X165dK/M7uXv3btGjRw+h0+lEWlqa8PHxEXl5eRKmfHjl1Xfw4EGRlZUlhBDi1KlTolGjRkKv10sV8ZGUV99dzz//vIiMjBSTJ0+WIJl5WeUegT/++ANNmzY13qlw8uTJFt2l38/LywsREREAADs7O4SHhxuf9CgXWq0Ws2fPxr///W+po5jF3Sd0vvLKKwDk94ROg8EAIQTy8/MBALm5ufDy8pI41aPr3r07mjRpUmbcli1bMGHCBCiVSnh7e6Nbt27YtWuXRAkfTXn1de3aFZ6engCAtm3bQqvVoqCgQIp4j6y8+gDg559/hlarRb9+/SRIZX7y3KdahbS0tDLPa27WrBlSU1MlTGQ+N2/eRGxsrMX+4anIBx98gMmTJ6Nhw4ZSRzGLS5cuoXHjxpgxY0aFT+i0ZC4uLvj888/Rtm1buLm5wc3NTbY3GEtLS8Pw4cONw3L+e7Np0ya0a9cOrq6uUkcxGY1Ggzlz5mD79u3YuXOn1HHMwir3CAgruWKyuLgYI0aMwKxZs9C6dWup45jMqVOncPjwYUycOFHqKGaj0+lw4sQJjBgxAseOHcOQIUMwadIkqWOZTElJCZYsWYI//vgDV69exYgRI/DGG29IHcssrOXvzfHjxzF79mysXLlS6igm9d5772Hq1Kmy/dABWGkj4OPjU6YjT0lJgVqtljCR6en1eowePRphYWF49dVXpY5jUgcPHkRCQgKaN28OPz8/6PV6+Pn54fbt21JHMxm5P6HzxIkTEEIYny8yatQoHDp0SOJU5mENf2+SkpLw7LPPYtOmTWjZsqXUcUzq0KFD+OCDD+Dn54fXXnsNX3/9NaZOnSp1LNOS9hQFaeh0OtG8efMyJwuuXLlS4lSmNWnSJDFhwgRhMBikjmJ2cj1ZsGfPnuLo0aNCCCFiY2NF165dJU5kOunp6cLT01OkpaUJIYT49NNPxXPPPSdxKtO593dy165dZU4WVKvV4vbt2xKme3T31peamioCAgLEjz/+KGEi06rob8qqVatkebKgVZ4joFQqsWLFCowYMQJarRa9evXC2LFjpY5lMgcPHsTKlSvRtm1bhIeHAwAmTZqEl156SeJkVBNLly4t84TOFStWSB3JZLy8vLBgwQIMGDAAtra2aNiwoSx2KU+bNg07d+6EXq+HWq3G008/jejoaOzevRuBgYGwsbHB4sWLLfYYenn1FRcX4/r163jjjTeMh3diY2Ph5+cnbdiHUNHPT+54i2EiIiIrZpXnCBAREVEpNgJERERWjI0AERGRFWMjQEREZMXYCBAREVkxNgJERERWjI0AkYXIy8tDWFgYwsLC0KxZM3h6eiIsLAyBgYHo06eP2da7d+9ezJw5s0bviY+Px/Tp0x9pvR988MFDvzc/Px/dunWDTqd7pAxE1oD3ESCyQDExMThw4ECt3GSoZ8+eWLNmTa3fIEalUj3Shvztt99GQEAAxo8fb8JURPLDPQJEFi45ORkBAQHG1/7+/pg+fTqCgoLQt29fHD9+HAMGDIC/vz8++eQT4/v27duH7t27o0OHDhg4cGC5T8S7ePEiNBqNsQmIiYnBkCFDMHjwYPj7++OVV17Bd999h65du6JFixb49ddfjcu++5yEffv2oVu3bhg9ejSCgoIwePBgFBUVAQB69+6NAwcOGNenUpXe7PTll1+GXq9HWFiYcTmnTp1C37590aFDB3Tv3h2nT58GAHz77bcICQlBWFgYQkJCcOXKFQDAs88+i9WrV5vs+0wkW9Le4ZiIHsa99zy/fPmyaNGihfG1QqEQhw8fFkII8cwzz4iuXbuKwsJCkZWVJTw8PIRGoxHZ2dmiS5cuIicnRwghxObNm8Xf/va3B9YTExMjpk+fXma9arVa3Lx5U2g0GtG0aVPx8ssvCyGE2Llzp+jdu7cQQoiff/5Z9OvXz/jayclJXLx4UQghxJNPPinWrVsnhBCiV69eYv/+/cbl33uP93tfFxcXi06dOhmfTXDkyBHRsWNHIYQQ7dq1E+np6UIIIQoLC4VGoxFClD5TxM3NTeh0upp+e4msilU+a4BIzry9vdGxY0cAQFhYGDQaDRwdHeHo6IgGDRogIyMDZ86cQWJiInr16gUAMBgM5d7/Pj09HY0aNSozrnfv3vDw8AAAtGrVCoMGDTKu6/Lly+Vmat++Pfz9/QEAjz32WIXzVSQxMRFnz57Fk08+aRx38+ZNY57IyEgMGzYMTz/9NJo1awag9Jkizs7OuHnzpqwfIUv0qNgIEMmMvb298bWNjc0DwzqdDkII9OjRA7GxsZUuy9HRETk5OdVa/t1lV5VJqVQa51OpVDAYDACAkpIS4+v7CSHQokULnDhx4oFp//3vf3H8+HHs3r0bvXr1wrp169CtWzcAQFFRERwdHSutkcja8RwBIivUpUsXHDlyBGfOnAFQuhG+e8z9XsHBwTh//rzZcjRv3hxHjx4FAGzZsgXinnOXnZycUFBQAABo3bo18vLyEBcXB6C0MTh+/DgAICkpCeHh4Xj99dcxYMAAY7OQnZ0NV1dXuLi4mC0/kRywESCyQg0aNMDGjRsxZcoUhIaGIiwszHii37169OiB06dPo6SkxCw5/vnPfyImJgZhYWFITEyEUqk0TpsxYwY6dOiA/v37w9bWFrGxsZg3bx5CQ0MRHByMLVu2AABef/11tG3bFmFhYcjMzERkZCQAIC4uDkOGDDFLbiI54eWDRFSp2bNnIyIiAiNGjJA6So08+eSTWLx4MVq1aiV1FKI6jXsEiKhSc+bMMdseAXPJz8/HqFGj2AQQVQP3CBAREVkx7hEgIiKyYmwEiIiIrBgbASIiIivGRoCIiMiKsREgIiKyYv8PIBQ3O960/oQAAAAASUVORK5CYII=\n"
          },
          "metadata": {}
        }
      ]
    },
    {
      "cell_type": "markdown",
      "id": "logical-acoustic",
      "metadata": {
        "id": "logical-acoustic"
      },
      "source": [
        "### Exercise 2\n",
        "\n",
        "Write an error function that simulates the temperature of the milk and returns the difference between the final temperature and 20 °C.  Use it to estimate the value of `r` for the milk."
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "from scipy.optimize import minimize_scalar"
      ],
      "metadata": {
        "id": "pjqeuQqBRH48"
      },
      "id": "pjqeuQqBRH48",
      "execution_count": 31,
      "outputs": []
    },
    {
      "cell_type": "code",
      "execution_count": 32,
      "id": "agreed-excuse",
      "metadata": {
        "id": "agreed-excuse"
      },
      "outputs": [],
      "source": [
        "# Solution goes here\n",
        "\n",
        "# Error function to be minimized\n",
        "def error_function(r, system):\n",
        "    \"\"\"Estimates the error between simulated final temperature and target temperature.\n",
        "\n",
        "    r: rate of cooling to test\n",
        "    system: System object with initial properties\n",
        "    \"\"\"\n",
        "    # Update the system's r value\n",
        "    system.r = r\n",
        "    # Run simulation with the current value of r\n",
        "    _, temperatures = run_simulation(system, change_func)\n",
        "    # Calculate the error from the target temperature (20°C)\n",
        "    final_temperature = temperatures[-1]\n",
        "    error = final_temperature - 20\n",
        "    return abs(error)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 33,
      "id": "defined-nudist",
      "metadata": {
        "id": "defined-nudist"
      },
      "outputs": [],
      "source": [
        "# Solution goes here\n",
        "\n",
        "# Using minimize_scalar to find the optimal r\n",
        "result = minimize_scalar(error_function, bounds=(0.01, 1), args=(milk,), method='bounded')\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 34,
      "id": "loose-rehabilitation",
      "metadata": {
        "id": "loose-rehabilitation",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 345
        },
        "outputId": "311d47a2-6fdb-48cd-eecf-fcad3874fd08"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Optimal value of r: 0.6219\n",
            "Minimum error achieved: 0.0\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 600x300 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgIAAAElCAYAAACWBIrfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAuJAAALiQE3ycutAAA7jElEQVR4nO3deVxU5f4H8M8wwyagLO4MgrghIA5qmIqC4lriFt5rarmWllup1+VmWdr9aZblTdFcrlKGGWWpSZaKS6hddzLkuqQiIIqCGyAMs5zfH8TkyCrMIud83q8Xr5hzzjzn+514+XznOc95jkwQBAFEREQkSTbWDoCIiIish4UAERGRhLEQICIikjAWAkRERBLGQoCIiEjCWAgQERFJGAsBogr4+Phg+/btAIDY2Fh07drVarGsX78eTZo0gbOzM86cOVPtdtLS0uDs7Iz79+8DAMaOHYs33ngDAJCamgqZTIZ79+6ZIOK/TJ48GXPnzi13/8GDB+Hq6mrSc1bXgAEDsHr16mq99/HP1pRUKhViYmJM3i4RCwGq9Q4fPowBAwbAzc0Nrq6uaN++PZYtW4aioiKTnmfUqFE4evSoSdusKo1Gg+nTpyMuLg55eXkIDg4udUx4eDhkMhn27dtntP3DDz+ETCYzdPbNmjVDXl4e6tWrZ4nQAQCfffYZPvjgAwCmKzZyc3Px5ptvwsvLC46OjmjRogUWLVoErVZb5TYeLYJK7N69G6+//nq1YrLGZ0tUUywEqFbbtWsXBgwYgH79+uHSpUu4d+8evv76a6SkpODGjRvWDs9kbt68icLCQrRr167C49q0aYNNmzYZbdu0aRP8/PzMGZ7FaTQa9OvXD2fOnMHevXuRl5eHuLg4fPvtt3jxxRetHZ7VaTQaa4dAtQgLAaq1BEHA9OnTMXfuXLzxxhuoX78+AMDPzw8xMTHw9vYGAJw8eRLdunWDq6sr/P398dVXXxm1sXz5crRo0QLu7u7o378/rly5Uub5YmJioFKpDK99fHywbNkyPPvss3BxcUFYWBjS09MN+8+dO2fY17NnT8yZMwfh4eHl5vPHH3+gX79+cHd3R4sWLbBixQoAwJkzZwwduVKpRIsWLcptY8SIEdi9e7dhaPrYsWMAgM6dOxuOeZJv5P/973/h6emJ7777rtS+r7/+Gs8++6zh9QsvvIAmTZoYXs+aNQvTpk0DYPzNOyQkxJCLs7MzYmNjDe/ZsGEDvLy84OHhgTlz5pQbV2xsLC5cuIAdO3bAz88PcrkcHTt2xPfff48dO3bg4MGDAIB3330XAwcOxIQJE1C3bl20atUK33//PQDg008/RWxsLFavXg1nZ2cEBAQAKB5ZKfnsSy5ZrF69Gp6ennBzc8OKFStw/vx5dO7cGXXr1sWQIUOQn59f6rO9desWnJ2djX5kMpkhtsuXLyMyMhINGjSAt7c33n//fej1ekOOq1atMnwWb731VoX/n0r+NhcuXIjGjRtjxIgRFR5P9CgWAlRrXbp0CVevXq3wG+C9e/fQv39/jBgxArdv38aaNWvwyiuv4MiRIwCAzZs34+OPP8b27duRmZmJgIAAREZGVnl4+csvv8RXX32F27dvw8nJCW+//TaA4m9kgwYNwoABA5CTk4OlS5di48aN5baj1WoxcOBAtG/fHpmZmfj++++xbNkybNmyBcHBwTh37hwAICMjA5cvXy63HVdXV/Tv399Q7GzcuBHjxo2rUi6Pi4+Px7Bhw7BlyxYMGzas1P7w8HCcOnUKubm5EAQBhw8fhoODA/73v/8BAPbv349evXqVet/x48cNueTl5WHUqFEAiof6U1JScOnSJRw+fBjR0dGGTvNxP//8M55//vlSQ/AtWrRA586dsWfPHsO2n376CSEhIbhz5w4+/vhjvPjii7h8+TKmT5+OUaNG4fXXX0deXp7hM35cbm4uUlNTcfXqVXzzzTeYPXs2Zs+ejW+//Rbp6en4448/sHbt2lLva9iwIfLy8gw/8+bNQ0BAADp06ICHDx8iIiICERERuH79OhITE7F161bDaM7+/fvx1ltvIS4uzjCylZycXGZ8JZKTk6FQKJCWlobNmzdXeCzRo1gIUK11+/ZtAICnp2e5x8THx6NBgwaYNm0abG1tERYWhpEjR+Lzzz8HUFwITJ8+He3atYODgwP+7//+D+np6YbOqjKvv/46mjdvDgcHB4waNQqnTp0CUPxNOicnB2+99Rbs7OzQuXNn/P3vfy+3nWPHjuHGjRt4//334eDggKCgIEydOrVak8PGjRuHTZs2oaCgANu2bcNLL730xG1s2rQJr732Gn766SeEhYWVeUyjRo3QunVrJCYmIikpCd7e3hg4cCAOHDiAO3fuIDk5ucIRkMcJgmDIv23btujatavh83xcdnY2mjZtWua+pk2bGv42AKB169aYNGkSFAoFIiMj0bNnT6NRoap47733YGdnh969e8Pd3R2RkZHw8vJCvXr18Nxzz+H06dMVvj8uLg7R0dHYtWsX6tati/j4eLi5ueGNN96AnZ0dmjVrhhkzZmDLli0Aikc8Ro0ahS5dusDOzg7vvvsunJycKjxHvXr1DH9vderUeaL8SNpYCFCtVXIp4Pr16+Uek5GRAR8fH6Ntvr6+yMjIKHO/vb09mjZtathfmcaNGxt+d3JyQm5uLgAgMzMTTZo0gUKhMOxv1qxZhXE2bdoUdnZ2Zcb5JCIiInDjxg0sXrwYXbp0MYqxqj744AOMHj0aQUFBFR7Xs2dPHDhwAPv370fPnj0RERGBAwcO4MCBAwgKCoKbm1uVz1m3bl2jDuzRz/Nx9evXR2ZmZpn7MjMz0aBBA8PrkktEj76u6G/mcS4uLnB0dDS8rlOnDho1amT0Oi8vr9z3//e//8XkyZPx/fffG/7WUlNTkZycDFdXV8PPrFmzcPPmTUMOj8Zta2trdNmlLJ6enrCx4T/p9OT4V0O1VuvWreHj44OtW7eWe4xSqURqaqrRttTUVCiVyjL3FxUVITMz07C/upo2bYqbN28aXWJIS0urMM7MzEyjSV6PxvkkbGxsMGbMGCxdurTalwV2796NLVu24MMPP6zwuEcLgV69eiE8PByJiYlISEhAz549y42vpvr06YMff/wRDx48MNp+9epVHDt2DH369DFsu3btmtExaWlphlEkc3ecqampGDJkCD777DOj+RReXl7o2LEj7t27Z/h58OCB4fJE06ZNjeLWaDSVTn5lEUDVxb8cqrVkMhlWrlyJpUuXYuXKlcjJyQEAXLx4ERMmTMC1a9fw3HPP4datW1i9ejW0Wi0SExMRGxuLl19+GQAwevRorFq1CikpKVCr1ViwYAE8PT0NE9qq69lnn4WrqyuWLFkCjUaDEydOIC4urtzjQ0JC0KhRI7zzzjtQq9VITk7GypUrMWbMmGqd/80338SePXsQGRlZrfc3b94chw4dwpo1a7BkyZJyjwsLC8Nvv/2GX3/9FaGhoXB1dYVSqURsbGyZ8wMAoEGDBrCxsalwrkNlRo8ejRYtWmDIkCG4cOECdDodTp8+jaFDh2LgwIFGRcjFixexfv16aLVaxMfHY//+/YbLNI0aNcKVK1dgjqexP3jwAAMHDsS0adPwt7/9zWjfwIEDkZWVhdWrV6OwsBA6nQ4XLlwwzIl48cUXERsbi2PHjqGoqAiLFi0yTEgkMjUWAlSrDRw4ELt370Z8fDxatGgBV1dXREVFwc/PD02aNIGbmxt2796NL7/8Eh4eHnj11VexZs0ahIaGAgBefvllTJs2DQMHDkTjxo3x22+/4YcffjAa0q8OW1tb7NixA7t27YKbmxvmzJmD0aNHw97evtzjd+3ahVOnTqFx48YYNGgQZs6ciZEjR1br/O7u7ujduzdsbW2rnYO3tzcOHTqE//znP1i8eHGZx9SvXx/+/v7w9/c3XMOOiIjAw4cP0aNHjzLf4+joiIULF2LAgAFwdXU1XBd/EnZ2dti7dy/atWuHXr16wcnJCVFRURg8eDC+/vpro2P79++P//73v3B3d8eMGTPw5ZdfolWrVgCAiRMn4vr163B3d6/0MsiTOn36NM6dO4clS5YY3TmQmJgIZ2dn7Nu3DwkJCfDx8YGHhwdGjhxpuDTQu3dvLF682HAnhl6vR2BgoEnjIyohE8xRChNRKZMmTYJer8f69eutHYpkvPvuu0hKSjKsDklEpXFEgMhMEhMTkZ6eDr1ej4SEBMTGxmL48OHWDouIyEjNxj+JqFxXrlzBiBEjcPfuXSiVSixduhR9+/a1dlhEREZ4aYCIiEjCeGmAiIhIwlgIEBERSVitmSNQt27dGi/y8riHDx+KfilOsefI/Go/seco9vwA8ecohvwyMjJKLcBVotYUAkqlEikpKSZtMyEhARERESZt82kj9hyZX+0n9hzFnh8g/hzFkJ+/v3+5+3hpgIiISMJYCBAREUkYCwEiIiIJYyFAREQkYSwEiIiIJMyshUB6ejoiIiLQtm1bBAQEYP78+YZ98+bNQ8uWLdG6dWts27bNnGEQERFROcx6+6BCocAHH3yATp06oaioCBEREdixYwecnJxw9OhRXLhwATdv3kSXLl3Qr18/ODs7mzMcIiIieoxZC4EmTZqgSZMmAIqfHx4cHIy0tDSkpKRg7NixkMvl8PT0RLdu3bBnzx4MGzbMnOE8NTQ6PQo1OhRq9FBrdVBr9dDpBaMfvSBAAFD8JAgBJU+EKHkwRFWfEHHxrh4uV++YPomnBPOr/cSeo9jzA8Sfo6Xzc7C1QZDS1WLns9iCQnfu3MH27duxZ8+eUp1+s2bNkJ6ebqlQLOJufhFOXruLk6l3kJqTj+y8ImTnqZGdq0Z+kc6ywZz81bLnszTmV/uJPUex5weIP0cL5te8vhMOzA632PksUggUFRUhKioKM2bMgJ+fH6rywMPo6GhER0cbXmdlZSEhIcGkceXn55u0zQdFAuKvapGSo0dm/l85utoDde1kqGsng6cH4Gwrh50csLWRwVYO2NoAchlgIwNsZDLYyAAZin/w2O+P/KdK1OpC2Ns7mCrFpw7zq/3EnqPY8wPEn6Ol83OQF5m8v6uI2QsBnU6HkSNHQqVSYdasWQAALy8voxGAtLQ0hISEGL1vypQpmDJliuG1v7+/yZd4NNWykXq9gG9OpeP/fjyPB4U6qLxcEdnRHZ183NHJ2w1uTnYmiLZ6xLA0ZkWYX+0n9hzFnh8g/hzFnp/ZC4FXX30VLi4uWL58uWHbsGHDsHjxYowZMwY3b97E4cOHsW7dOnOHYhaXsnLx1vfJOJ56B36NXRAz7hkEN3OzdlhERERVYtZC4MiRI9i4cSMCAwMRHBwMABg/fjymT5+OvXv3onXr1rCxscHHH38MFxcXc4ZiFomXbmN8zAnIbWSYP8AP40Obw1bOpRmIiKj2MGsh0K1bt3LnAyxbtgzLli0z5+nNqlCjw1vfJ8PDyR7fTO4CL/fa/YhKIiKSJn59raY1By8j7c5DvBPpzyKAiIhqLRYC1XAtJx9rDl1G91b1MSCwsbXDISIiqjYWAk9IEAQs3HkOEID3BgVAJnuSm/mIiIieLiwEntCelCwcvHAbr/bwhW8DLolMRES1GwuBJ/CwSItFP6TA09URU3q2tHY4RERENWaxJYbFYEPiVVy/V4D1L3eCo53c2uEQERHVGEcEnkDC+Vto08gFffwbWTsUIiIik2AhUEUFRTqcu34fzzTnqoFERCQeLASq6GzGPWj1Ajp6sxAgIiLxYCFQRSev3QUAdPJ2t3IkREREpsNCoIpOX7uLBi72ULo5WjsUIiIik2EhUAV6vYBTaXfRyduNCwgREZGosBCogivZ+bj3UMP5AUREJDosBKrg1LU7AMBCgIiIRIeFQBWcTL0Le4UNAprWs3YoREREJsVCoApOpd1Fe6Ur7BT8uIiISFzYs1XiTn4RrtzORwdeFiAiIhFiIVCJ04b1A1gIEBGR+LAQqETJQkIcESAiIjFiIVCJ09fuwreBE9yd7KwdChERkcmxEKhAkVaP3zLuoWMzjgYQEZE4sRCowLnM+1Br9ejkw0KAiIjEiYVABU79OT+ACwkREZFYsRCowKlrd+Faxxa+9Z2tHQoREZFZsBCowNmM+1B5ucLGhg8aIiIicWIhUA69XkDWg0J4uvKxw0REJF4sBMrxoFADrV6Ah7O9tUMhIiIyGxYC5cjOKwIANHDm+gFERCReZi0EZsyYAaVSCYVCYdj24MED/P3vf0e7du0QFBSEnTt3mjOEasvJUwMARwSIiEjUzFoIDB8+HCdPnjTa9q9//QtKpRK///47fvnlFyxYsAB5eXnmDKNacvKLRwQ8uKIgERGJmFkLgdDQUDRu3NhoW3JyMvr37w8AcHV1hZ+fH3bv3m3OMKqFIwJERCQFFp8jEBwcjLi4OOj1ely/fh2HDx9Genq6pcOoVMkcgfqcI0BERCImEwRBMPdJFAoFtFotACA/Px8zZ87EsWPH4Onpibp16+LZZ5/FjBkzjN4THR2N6Ohow+usrCzExcWZNK78/Hw4OTmVuW/z/zRIvK7DZxH2sJHV3nUEKspRDJhf7Sf2HMWeHyD+HMWQ37Rp05CSklL2TsEC5HJ5ufsGDRok/PTTT5W20bZtW1OGJAiCIOzbt6/cfZO+OCmE/Guvyc9paRXlKAbMr/YTe45iz08QxJ+jGPKrqA+1+KWB+/fvQ60uvv5+5MgRXL16FX369LF0GJXKyVfDw4nzA4iISNwUlR9SfZMmTUJ8fDx0Oh2USiUGDx6MMWPGYNSoUVAoFGjcuDG++eYb2Ng8fcsZ5OQVwdONqwoSEZG4mbUQWLt2bZnbL126ZM7TmkR2nhrtvVytHQYREZFZPX1fxZ8CRVo9HhRquYYAERGJHguBMtwpWUyIawgQEZHIsRAoQ7ZhMSGOCBARkbixEChDyfLCDTgiQEREIsdCoAzZuRwRICIiaWAhUIacfD5ngIiIpIGFQBly8vjkQSIikgYWAmXIziuCs70CDrZya4dCRERkViwEypCTr+b8ACIikgQWAmXIySviZQEiIpIEFgJlyMlToz4nChIRkQSwEHiMIAjIzi/iHQNERCQJLAQek6vWokirR33OESAiIglgIfAY3jpIRERSwkLgMTl5XEyIiIikg4XAY7JLRgR4aYCIiCSAhcBjSpYX5l0DREQkBSwEHlMyR4CFABERSQELgcfk5KlhIwNcHW2tHQoREZHZsRB4THZ+Edyd7GFjI7N2KERERGZXYSGQmJiICRMmoHXr1nBxcUHDhg3Rq1cvfPLJJ7h3756FQrSs7Fw11xAgIiLJKLcQGDBgADZs2IBBgwbh4MGDuHv3Lq5cuYKPPvoIWq0Wzz//POLj4y0Zq0Xk5BfxjgEiIpIMRXk71q1bBy8vL6Ntzs7O6NChAzp06IDZs2cjMzPT7AFaWk6eGv5N6lo7DCIiIosod0TAwcEBKSkppbanpKQgOzsbMpkMnp6eZg3O0rQ6Pe4+1PCOASIikoxyC4GZM2ciKyur1PZbt27hzTffNGtQ1nLnIRcTIiIiaSm3ELhw4QJ69uxZant4eDiSkpLMGZPV/LWGAAsBIiKShnILgfv375f7pqKiIrMEY21/PXCIlwaIiEgayi0E6tevj6NHj5ba/uuvv8Ld3d2sQVlLyfLCvDRARERSUe5dA++//z6GDBmCyZMno1OnTgCAkydPYu3atfjqq6+q1PiMGTOwbds23Lx5E1qtFgCgVqvxyiuv4MyZM9BqtYiIiMCnn34KGxvrr210O5fPGSAiImkpt/ft2bMnfvnlF6SlpWHRokVYtGgR0tLScPDgQfTq1atKjQ8fPhwnT5402vbFF19Aq9Xi999/R3JyMpKSkrB79+6aZWEiOfmcLEhERNJS7ogAAPj5+SEmJqbajYeGhpa5PT8/H1qtFkVFRVCr1WjcuHG1z2FKOXlq1LGTo45dhR8LERGRaJQ7IlBQUIANGzZg/fr1KCgoMNkJX375ZdjZ2aFRo0Zo1KgR+vbti44dO5qs/ZrIyeOqgkREJC0yQRCEsnZERUVhzJgxkMlk2LRpE7Zt21btkygUCsMcgcTERHz88cfYunUrCgoKEB4ejtWrV6Nr165G74mOjkZ0dLThdVZWFuLi4qodQ1ny8/Ph5ORkeP3+MTVkMuCtEPHMEXg8R7FhfrWf2HMUe36A+HMUQ37Tpk0rc5FAAIBQjm7dugm3bt0Sbt26JXTp0qW8w6pELpcbfn/99deFdevWGV7/4x//EJYvX15pG23btq1RDGXZt2+f0etuSxOECTHHTX4ea3o8R7FhfrWf2HMUe36CIP4cxZBfRX1ouZcGNmzYgPnz52PevHn47LPPTFaVeHt7IyEhAXq9Hmq1GomJifD39zdZ+zWRk1fENQSIiEhSyp0V5+fnhw0bNtSo8UmTJiE+Ph46nQ5KpRKDBw/GsmXLMGHCBAQGBkIQBAwZMgT9+/ev0XlM4WGRFgUaHecIEBGRpJRbCJw+fRodOnQo942FhYW4du0a2rRpU+4xa9euLXP71q1bnyBEyzCsKsg1BIiISELKLQT+/e9/IzMzE8OGDUOnTp3QsGFDFBYW4tKlS9i3bx+OHj2K5cuXV1gI1CYPCjUAgLoOvHWQiIiko9xe7/PPP8cff/yBzZs3Y+7cucjMzISjoyMCAwMRGRmJDz/8ELa2tpaM1azUWj0AwN5WbuVIiIiILKfCr78tW7bEe++9Z6lYrEqt+bMQUFh/qWMiIiJLYa/3p0KtDgDgwBEBIiKSEBYCf+KIABERSRF7vT+p/xwRYCFARERSUqVe7+eff8aqVasAFC/1e/HiRbMGZQ0lIwK8NEBERFJSaSGwYMECrF69GitXrgQA6PV6jB071txxWRxHBIiISIoq7fV27dqF7du3o06dOgCAJk2aIC8vz+yBWVqhhrcPEhGR9FRaCDg4OEAmkxleq9VqCGU/sLBWKxkRcOCIABERSUilvV7v3r3x9ttv4+HDh9i5cyeGDBmCqKgoS8RmUVxQiIiIpKjSQmDx4sXw9fVFUFAQvvjiC0RFRWHhwoWWiM2iCjUcESAiIumpcGVBnU6Hvn37IiEhAePGjbNUTFah1uoht5FBIWchQERE0lFhryeXy6HValFQUGCpeKxGrdHzjgEiIpKcSh+15+XlhZCQEAwaNAhOTk6G7f/85z/NGpilFWp1XEOAiIgkp9JCoFWrVmjVqhUAQKPRmD0ga+GIABERSVGlhYAYJwaWpVCrYyFARESSU2khsGjRojK3v/POOyYPxprUGj0vDRARkeRU+hVYLpcbfjQaDXbs2IGrV69aIjaLUnNEgIiIJKjSEYG33nrL6PX8+fMRGRlptoCspVCjh7NDpR8HERGRqDzxV2C5XI60tDRzxGJVHBEgIiIpqvQrcJ8+fQzPGtDr9Th//jzGjBlj9sAsTa3Vw17BOQJERCQtlRYCCxYs+OtghQI+Pj7w9PQ0a1DWUKjRw8GWIwJERCQtlfZ8P//8M8LCwhAWFoZu3brB09NTdIsJASWXBjgiQERE0lJpIfDTTz+V2vbjjz+aJRhrUmv0sOeIABERSUy5lwY2b96ML774ApcvX0bfvn0N2x88eICGDRtaJDhL0esFFOn0cOCIABERSUy5hUCPHj3g5eWFmTNnGt1C6OLigvbt21skOEsp0ukBgCMCREQkOeUWAt7e3vD29sbp06ctGY9VFGp0AMDbB4mISHIq7fmuXbuGkSNHIigoCK1btzb8VMWMGTOgVCqhUPxVb2zYsAEqlcrwY2dnhx07dlQ/AxNQa4tHBLjEMBERSU2lhcC4ceMwfPhw6HQ6fPfdd+jXrx9GjhxZpcaHDx+OkydPGm2bOHEikpKSkJSUhB9++AFOTk7o169f9aI3EbXmz0sDHBEgIiKJqbTnu3//PoYOHQq5XI7AwECsXLkS8fHxVWo8NDQUjRs3Lnf/li1bMGzYMDg4OFQ9YjMo1BZfGuCIABERSU2lhUDJsH79+vVx/PhxZGRkICcnxyQnj42NxUsvvWSStmqCIwJERCRVMkEQhIoOiImJweDBg3H+/HmMGjUKubm5WLZsGcaNG1flkygUCmi1WqNtv/32GyIjI3Ht2jXDEsaPio6ORnR0tOF1VlYW4uLiqnzOqsjPz4eTkxMu3dVj6ckiTA6yxTONxDUqUJKjWDG/2k/sOYo9P0D8OYohv2nTpiElJaXsnUIFtFqt8M0331R0SJXI5fJS22bPni3Mmzevym20bdu2xnE8bt++fYIgCELixduC99xdwr6UmyY/h7WV5ChWzK/2E3uOYs9PEMSfoxjyq6gPrXAsXC6XY8mSJSavTPR6Pb766qun4rIAULy8MAAuMUxERJJT6UXx5557Dp9//jmKioqeuPFJkyZBqVRCp9NBqVRiypQpAICDBw+iYcOG8Pf3f/KIzaBQwwWFiIhImip9+uDSpUuh0+kwbtw4KBQKCIIAmUxWpcJg7dq1ZW7v1avXU7VQUcmIAJcYJiIiqam0ENBoNJaIw6pKFhTiiAAREUlNlXq+n3/+GatWrQIA3Lp1CxcvXjRrUJbGJYaJiEiqKu35FixYgNWrV2PlypUAAJ1Oh7Fjx5o7LoviEsNERCRVlRYCu3btwvbt21GnTh0AQJMmTZCXl2f2wCyJIwJERCRVlfZ8Dg4ORgv+qNVqCBWvQVTrcESAiIikqtJCoHfv3nj77bfx8OFD7Ny5E0OGDEFUVJQlYrOYkiWG7eQcESAiImmptOdbvHgxfH19ERQUhC+++AJRUVFYuHChJWKzmEKtDnZyG9jYlF7qmIiISMwqvX1QJpNh3LhxGDFiBADA0dHR7EFZmlqj562DREQkSZX2fsnJyXjmmWfg5eUFLy8vhISEIDk52RKxWYxaq+PywkREJEmVFgJjx47F22+/jezsbGRnZ+Odd94R3e2DhRo97xggIiJJqrT3KyoqwqBBgwyvBw4cKLrVBtVaHRx4aYCIiCSo0t4vMjISa9euRUFBAQoKCrBu3TqjwkAM1Bo9Lw0QEZEkVTpZcNmyZdDpdHjttdcM2+RyOZYuXVrlhw897TgiQEREUsWHDqF4QSG3OnbWDoOIiMjiKi0EgOJ5Ajdv3oRWqzVs8/X1NVtQllao0fH2QSIikqRKC4EVK1Zg8eLFaNCggWGpYZlMhpSUFLMHZylqrR4OnCNAREQSVGkhsGrVKly6dAnu7u6WiMcq1FouKERERNJUae/n6+sLZ2dnS8RiNYUaHdcRICIiSarSXQMRERHo0aMH7O3tDdvfeecdswZmSWqtnk8eJCIiSaq0EJg9ezYaNmwIBwcH2NiI71uzVqeHTi9wRICIiCSp0kLgxo0b2LdvnyVisYpCbfEjiLmgEBERSVGlX4N79OiBU6dOWSIWq1BrdADABYWIiEiSKh0R2LNnD9avXw9PT0/Y29tDEATIZDJcvHjREvGZHUcEiIhIyiotBPbv32+JOKyGIwJERCRllfZ+3t7eKCgoQFJSEry9veHh4QEXFxdLxGYRao4IEBGRhFVaCKxevRovvfQSZs2aBQC4efMmoqKizB6YpRT+OSLABYWIiEiKKu391q5diyNHjhhGAVq2bInbt2+bPTBL4YgAERFJWbmFwJEjRwAAjo6OsLP768l8er3e/FFZkKEQ4IgAERFJULm937Rp0wAAwcHB2LBhA7RaLc6ePYuJEyeiZ8+eVWp8xowZUCqVUCiM5yQeO3YMISEhCAgIQEBAADIzM2uQQs0YLg1wQSEiIpKgcns/QRAAAJ988gnS09Ph4OCACRMmQKlU4qOPPqpS48OHD8fJkyeNtuXm5uLll1/G559/jnPnzuHo0aNwc3OrQQo1UzIiwCWGiYhIisq9fTA9PR2vvvoqgOKiQKVSASheaXDq1KlYt25dpY2HhoaW2rZlyxY899xzaNu2LQCgXr161YnbZDgiQEREUlZuIeDk5IRu3bqZ/ITnz5+HWq1G7969kZOTg0GDBuHdd9+FTCYz+bmqgpMFiYhIymRCyTWAx3To0AGnT582yUkUCgW0Wi2A4rkHBw8eRGJiIhwcHDB48GCMHj0aL730ktF7oqOjER0dbXidlZWFuLg4k8RTIj8/H0ey7fH1RS0+CLVHfUfrFCPmlJ+fDycnJ2uHYTbMr/YTe45izw8Qf45iyG/atGlISUkpc1+5IwLl1Ac15uXlhb59+8LV1RUAMHjwYJw6dapUITBlyhRMmTLF8Nrf3x8REREmjSUhIQHNXLyBixfQK6w7GrjYV/6mWiYhIcHkn9vThPnVfmLPUez5AeLPUez5lXthPDEx0SwnHDZsGI4cOQK1Wg29Xo+EhAQEBASY5VxVwSWGiYhIysrt/ZydnWvc+KRJk6BUKqHT6aBUKjFlyhS0bNkSL730EoKDgxEUFIRGjRph/PjxNT5XdXGOABERSVmlDx2qibVr15a5/fFhf2sq1OggkwG2cvHNDyAiIqqM5MfD1Vo9HBRyq921QEREZE2SLwQKNTouL0xERJIl+R5QrdVzMSEiIpIsyfeAaq2eywsTEZFkSb4QKNToOCJARESSJfkekCMCREQkZSwEtBwRICIi6ZJ8D1io0XMxISIikizJFwJqrY7LCxMRkWRJvgfkiAAREUmZ5AsBNe8aICIiCZN8D6jW6mHPuwaIiEiiJF0ICILAlQWJiEjSJN0D/vkEYj5rgIiIJEvSPaDmz0LAgZMFiYhIoiRdCBRxRICIiCRO0j2gRi8A4IgAERFJl7QLAV3xfzkiQEREUiXpHrBkjgAXFCIiIqliIQBwiWEiIpIsSfeARbriOQIcESAiIqmSdCHw16UBSX8MREQkYZLuAf+6NMARASIikiYWAuCIABERSZeke0BNyRwBThYkIiKJknQPyCWGiYhI6lgIgCMCREQkXZLuATkiQEREUmfWQmDGjBlQKpVQKBSGbQcPHoSLiwtUKhVUKhWGDh1qzhAqVMQ5AkREJHGKyg+pvuHDh2P+/PlQKpVG2zt37ox9+/aZ89RVwiWGiYhI6sxaCISGhpqz+RrT6AFbuQxyG5m1QyEiIrIKq4yJnzp1CsHBwejRowd+/vlna4QAoLgQ4GgAERFJmUwQBMHcJ1EoFNBqtQCABw8eAADq1q2Lc+fOoX///vjll1/QvHlzo/dER0cjOjra8DorKwtxcXEmjSv6zENcum+DFeEOJm33aZKfnw8nJydrh2E2zK/2E3uOYs8PEH+OYshv2rRpSElJKXOfxQuBx40YMQLDhw/HCy+8UGEb/v7+5SZRXUOX78YtjT2OzOtl0nafJgkJCYiIiLB2GGbD/Go/seco9vwA8ecohvwq6kMtfmngxo0bKKk9rl+/jl9//RUBAQGWDgNAyaUB3jFARETSZdbJgpMmTUJ8fDx0Oh2USiUGDx6Mtm3bYs2aNbC1tQUAvP/++/Dz8zNnGOXS6AXYO3COABERSZdZC4G1a9eWuX3q1KnmPG2VFemAehwRICIiCZN0L8hLA0REJHWS7gU1esDBlpcGiIhIuiReCAgcESAiIkmTdC+o0QH2HBEgIiIJk3YhoAccOCJAREQSJuleUKPnkweJiEjazHr74NNMq9NDJ/BZA0RUzJyLrFpgAVerE3uOtSU/mezJH6In2UJArS1+BrEDRwSIJE2j0SA9PR1qtdos7bu7u+P8+fNmaftpIfYca1N+9vb28PLyMizaVxWSLwQ4IkAkbenp6XBxcYGPj0+1vk1V5sGDB6hbt67J232aiD3H2pKfIAjIyclBeno6fH19q/w+CRcCOgAcESCSMkEQoFar4ePjAxsb8/xbIJPJzFJgPE3EnmNtyU8mk8HDwwPZ2dkQBKHKMUu2FyzUcESAiIrVhn/kiaqiOn/Lki0ESkYEuKAQEVlbbm4uVCoVVCoVmjVrBg8PD8PrxMREq8V15MgRBAYGIjg4GOnp6Ub7wsPD0aBBA6NHzM+aNQsymQypqakAgOeeew63bt0CAPj4+CAjI6PU7zXRqVMnw++LFi0y2meqc5QlOTkZHTt2RKtWrTBkyBDk5eWVeVx+fj7GjBmDNm3awM/Pz/D8ndjYWLRv3x5BQUHo1KkT9u/fb3hPTEwM/P39YWNjg8OHDxu1t3DhQgQEBCAgIADvvfeeyfKRbC+o1pRMFuSIABFZl4uLC5KSkpCUlIRFixZh6NChhtfdu3e3WlybN2/GG2+8gTNnzsDLy6vUfk9PT/z8888AAL1ejz179hgd9+OPP6Jhw4Zmi+/kyZOG3x8vBKrj0aKmIpMnT8aSJUtw6dIltG7dGsuXLy/zuFmzZiEgIAAXLlzA//73PwwdOhQA0Lx5c+zfvx9nz57Fpk2b8OKLL0KvL+6TOnbsiB07dqBHjx5Gbe3evRt79+7F6dOncfbsWfz222/45ZdfapDtXyRbCBRqOCJARE+vgoICTJ48GSEhIWjXrh1WrVoFAEhNTYWvry8mT54Mf39/9OrVC2fOnEGfPn3g6+uLTz75xNCGQqHAP//5T7Rr1w4dOnTA77//Xua51qxZg8DAQLRr1w7Tpk2DVqtFdHQ04uLisHjxYvTr16/M940cORJbtmwBABw4cABdu3aFnZ2dYX9l38p1Oh3GjBmDOXPmGG3X6/Vo2rQp1Go19Ho93N3dER0dDQB49913DZ+FQlE8ze3NN9+ETqeDSqVC7969De2sW7cOnTp1QuvWrXH06NEyY1AoFFiwYAGCg4Px3XfflRtriaysLKSlpaFv374AgAkTJmDbtm2ljsvNzcXOnTsxc+ZMAMVD9iVFUdeuXeHh4QEACAwMhFqtRn5+PgCgXbt2aNWqVan2kpOT0bNnT9jb20Mul6N37974+uuvK423KiQ8WfDPOQKcLEhEf5rz7W+4mFX2MG916XQ6tG1aD8ui2j/R+5YsWYIOHTrgs88+g1qtRrdu3dCrVy/UqVMHqampGD9+PD777DMMGzYMU6dOxb59+/Dw4UO0atUKr732GhwcHKDT6aBUKvH7779jx44dGDdunNG3aABISkrC8uXLceLECdSrVw/Dhg3D+vXrMWXKFJw4cQK9e/fG6NGjy4yxa9eu2Lx5Mx4+fIjY2FiMGTMGCQkJVcpPrVYjKioKzzzzDP75z38a7bOxscEzzzyDY8eOoV69emjVqhUOHTqEKVOm4ODBg4ZCoMQnn3yClStXIikpyWh7nTp1cPLkSXz77bdYuHAh9u7dWyoOnU4HX19fnDlzBgDw6aefYuPGjUbH6PV6ODo64tixY8jIyDAa9WjWrFmpyyYAcOXKFTRq1AhTp07F8ePH4e3tjRUrVsDb29vouK1bt6Jdu3ZwcXGp8PMKDg7GnDlzMGfOHNjb22PHjh2wt7ev8D1VJdlesGREwIGTBYnoKfTjjz/i3//+N1QqFTp37ow7d+7gwoULAIqH5ENCQgAAKpUKPXr0gKOjIzw8PFC/fn3cvHnT0M6YMWMAAIMHD8aVK1cM3zxLHDp0CEOHDoWbmxtsbGwwduxYHDhwoMpxRkZG4rvvvsPx48dLDWdXJCoqCt27dy9VBJQICwvDoUOHcOjQIbzyyiu4dOkSCgsL8ccffyAgIKBK5xg+fDgA4JlnnsHVq1fLPW7kyJGG36dPn264LFPyc/jwYRw7dgxA1RcW0mq1SEpKQlRUFE6fPo3IyEiMHz/e6JgzZ85g3rx5pQqPsvTu3RujRo1CeHg4+vTpg3bt2kEuN03/xREBjggQ0Z+e9Ft7VVT3HnRBEBAbGwuVSmW0PTU11eiboI2NTanXVb3WDZSeZf6ks85HjhyJHj16YMKECU/03u7du2PPnj2YMmVKmd9sw8PDMWfOHNSrVw8ffPABEhISEBMTg5CQkCqfp6RduVxe7mcil8vh4OBgeF3ZiIBSqTQaAUhLS4NSqSzVrlKphIeHh+FSxYgRIwyXCQDg4sWLeOGFF7B169YyLwWUZdasWZg1axaA4lGQJ1k0qCKS7QW5oBARPc0GDBiAf//739DpikcvL126hAcPHjxxO5s3bwYA7Nq1C76+vnBycjLa36NHD2zfvh337t2DXq/H559/jl69elW5/cDAQMyfPx+vvfbaE8U1Z84chIaG4oUXXkBRUVGp/SqVCikpKbh8+TJatmyJ8PBwLFmyBOHh4WW2V6dOnVKjHdVR2YhA48aN4eXlhT179gAA/vOf/2DYsGGl2mnUqBECAgJw+vRpAMDevXsNIxkZGRl4/vnnsWbNGnTp0qVKcen1ety+fRsAcOPGDaxbtw6TJk2qcb6AhAsBw6UBjggQ0VNowYIFcHZ2Rvv27REYGIhXXnmlzA6zInK5HOnp6QgKCsI777xT5hC0SqXCzJkzERoaiqCgIDRt2hQTJ058ovNMmjTpiVayK7FgwQJ07NgRw4cPh0ajMdpnY2ODDh06IDAwEEDxpYK0tDSEhYWV2dbUqVPRsWNHo8mC5rJmzRrMnTsXrVq1wvnz5zF79mwAxXcxPPfcc0bHTZkyBUFBQVi+fDk2bNgAAHjvvfdw69YtzJ0713CbaMktlzExMVAqlfj1118xbNgwdOjQAUDxXIawsDAEBASgX79++PDDD+Hj42OSfGRCLXmSgr+/P1JSUkzW3n8OX8XiXSlInNMTXu51TNbu0yYhIQERERHWDsNsmF/tZ80cBUHA+fPn4efnZ7ZFhay5PK1CoXiiywTVVVuW4K2u2pRfeX/TFfWhkp0jEOHXENlpl1Df2TSzLomIiGojyY6L+9R3QqdGcjjacY4AEYmTJUYDqPaTbCFARERELASIiKp8bzjR0646f8uSnSNARCSTyWBvb4+cnBx4eHiYZcKgIAiiLzTEnmNtyU8QBOTk5MDe3v6J/pZZCBCRpHl5eSE9PR3Z2dlmab+wsNBowRoxEnuOtSk/e3v7Mh8QVREWAkQkaba2tvD19TXbN779+/c/0QI9tZHYc6xN+VVnVIuFABERqvcP6NPQ9tNC7DmKOT9OFiQiIpIwFgJEREQSVmuWGK5bt26ZT3iqibt378LNzc2kbT5txJ4j86v9xJ6j2PMDxJ+jGPLLyMgo96FVtaYQMAdTP7/gaST2HJlf7Sf2HMWeHyD+HMWeHy8NEBERSRgLASIiIgmTdCEwZcoUa4dgdmLPkfnVfmLPUez5AeLPUez5SXqOABERkdRJekSAiIhI6lgIEBERSZhkC4GDBw8iICAALVu2xMSJE6HT6awdksmkp6cjIiICbdu2RUBAAObPn2/tkMxmypQpUCjEuVJ2fn4+xowZgzZt2sDPzw9r1661dkgm9eWXXyIoKAgqlQrdu3fHhQsXrB1Sjc2YMQNKpbLU3+S8efPQsmVLtG7dGtu2bbNSdDVXVn6xsbFo3749goKC0KlTJ+zfv9+KEdZMef//AODevXvw9PTExIkTrRCZmQkSpNPphBYtWgjnzp0TBEEQhg8fLsTExFg5KtPJzMwUTpw4IQiCIKjVaiE0NFTYvn27laMyvV9++UV4+eWXBblcbu1QzGLSpEnCBx98IAiCIOj1eiErK8vKEZlOfn6+4O7uLty+fVsQBEFYs2aNEBUVZeWoai4xMVG4ceOG0d/k3r17he7duwtarVbIyMgQvLy8hNzcXCtGWX1l5XfkyBEhOztbEARBOHv2rNCwYUNBp9NZK8QaKSu/Eq+88oowevRoYcKECVaIzLwkOSJw4sQJNG3aFP7+/gCACRMm1Ooq/XFNmjRBp06dAAB2dnYIDg5GWlqalaMyLbVajXnz5uGjjz6ydihmkZubi507d2LmzJkAih940rBhQytHZTp6vR6CICAvLw8AcP/+fTRp0sTKUdVcaGgoGjdubLRt27ZtGDt2LORyOTw9PdGtWzfs2bPHShHWTFn5de3aFR4eHgCAwMBAqNVq5OfnWyO8GisrPwA4cOAA1Go1IiIirBCV+YlzTLUSGRkZRs9rbtasGdLT060YkfncuXMH27dvr7X/8JRn0aJFmDBhAho0aGDtUMziypUraNSoEaZOnYrjx4/D29sbK1asgLe3t7VDMwlnZ2esWrUKgYGBqFevHurVq4dff/3V2mGZRUZGBoYNG2Z4LeZ/b7Zu3Yp27drBxcXF2qGYTEFBAebPn48ffvgB8fHx1g7HLCQ5IiBI5I7JoqIiREVFYcaMGfDz87N2OCZz9uxZHDt2DOPGjbN2KGaj1WqRlJSEqKgonD59GpGRkRg/fry1wzIZjUaD1atX48SJE7h+/TqioqIwd+5ca4dlFlL59+bMmTOYN28eNm7caO1QTOrdd9/Fq6++KtovHYBECwEvLy+jijwtLc3kDzSyNp1Oh5EjR0KlUmHWrFnWDsekjhw5gpSUFDRv3hw+Pj7Q6XTw8fEp94EatZFSqYSHhwd69+4NABgxYgROnTpl5ahMJykpCYIgoG3btgCK8zt69KiVozIPKfx7c/HiRbzwwgvYunUrWrVqZe1wTOro0aNYtGgRfHx8MHv2bHz99dd49dVXrR2WaVl3ioJ1aLVaoXnz5kaTBTdu3GjlqExr/PjxwtixYwW9Xm/tUMxOrJMFe/ToIZw6dUoQBEHYvn270LVrVytHZDqZmZmCh4eHkJGRIQiCIKxYsUL4+9//buWoTOfRv8k9e/YYTRZUKpXCgwcPrBhdzT2aX3p6utCyZUvhp59+smJEplXevymbNm0S5WRBSc4RkMvl2LBhA6KioqBWqxEWFoaXXnrJ2mGZzJEjR7Bx40YEBgYiODgYADB+/HhMnz7dypHRk1izZg0mTJiA/Px8uLq6YsOGDdYOyWSaNGmCpUuXok+fPrC1tUWDBg1EMaQ8adIkxMfHQ6fTQalUYvDgwYiOjsbevXvRunVr2NjY4OOPP66119DLyq+oqAi3bt3C3LlzDZd3tm/fDh8fH+sGWw3l/f8TOy4xTEREJGGSnCNARERExVgIEBERSRgLASIiIgljIUBERCRhLASIiIgkjIUAERGRhLEQIKolcnNzoVKpoFKp0KxZM3h4eEClUqF169bo2bOn2c67f/9+TJs27Ynec/LkSUyePLlG5120aFG135uXl4du3bpBq9XWKAYiKeA6AkS1UExMDA4fPmyRRYZ69OiBL774wuILxCgUihp15G+//TZatmyJMWPGmDAqIvHhiABRLZeamoqWLVsafvf19cXkyZPh7++PXr164cyZM+jTpw98fX3xySefGN538OBBhIaGomPHjujbt2+ZT8S7fPkyCgoKDEVATEwMIiMjMWDAAPj6+mLmzJnYuXMnunbtihYtWuCXX34xtF3ynISDBw+iW7duGDlyJPz9/TFgwAAUFhYCAMLDw3H48GHD+RSK4sVO33zzTeh0OqhUKkM7Z8+eRa9evdCxY0eEhobi999/BwB8//33CAoKgkqlQlBQEK5duwYAeOGFF/D555+b7HMmEi3rrnBMRNXx6JrnV69eFVq0aGH4XSaTCceOHRMEQRCGDh0qdO3aVXj48KGQnZ0tuLm5CQUFBUJOTo7QpUsX4d69e4IgCEJcXJzwt7/9rdR5YmJihMmTJxudV6lUCnfu3BEKCgqEpk2bCm+++aYgCIIQHx8vhIeHC4IgCAcOHBAiIiIMv9epU0e4fPmyIAiC8PzzzwtffvmlIAiCEBYWJiQmJhraf3SN90d/LyoqEjp37mx4NsHx48eFkJAQQRAEoV27dkJmZqYgCILw8OFDoaCgQBCE4meK1KtXT9BqtU/68RJJiiSfNUAkZp6enggJCQEAqFQqFBQUwNHREY6Ojqhfvz5u3ryJ5ORkXLhwAWFhYQAAvV5f5vr3mZmZaNiwodG28PBwuLm5AQDatGmDfv36Gc519erVMmPq0KEDfH19AQDPPPNMuceV58KFCzh37hyef/55w7Y7d+4Y4hk9ejSGDBmCwYMHo1mzZgCKnyni5OSEO3fuiPoRskQ1xUKASGTs7e0Nv9vY2JR6rdVqIQgCunfvju3bt1fYlqOjI+7du1el9kvariwmuVxuOE6hUECv1wMANBqN4ffHCYKAFi1aICkpqdS+Tz/9FGfOnMHevXsRFhaGL7/8Et26dQMAFBYWwtHRscIciaSOcwSIJKhLly44fvw4kpOTARR3wiXX3B8VEBCAS5cumS2O5s2b49SpUwCAbdu2QXhk7nKdOnWQn58PAPDz80Nubi4SEhIAFBcGZ86cAQBcvHgRwcHBmDNnDvr06WMoFnJycuDi4gJnZ2ezxU8kBiwEiCSofv36+OqrrzBx4kS0b98eKpXKMNHvUd27d8fvv/8OjUZjljj+8Y9/ICYmBiqVChcuXIBcLjfsmzp1Kjp27IjevXvD1tYW27dvx/vvv4/27dsjICAA27ZtAwDMmTMHgYGBUKlUyMrKwujRowEACQkJiIyMNEvcRGLC2weJqELz5s1Dp06dEBUVZe1Qnsjzzz+Pjz/+GG3atLF2KERPNY4IEFGF5s+fb7YRAXPJy8vDiBEjWAQQVQFHBIiIiCSMIwJEREQSxkKAiIhIwlgIEBERSRgLASIiIgljIUBERCRh/w/PmKGlvAWJ/gAAAABJRU5ErkJggg==\n"
          },
          "metadata": {}
        }
      ],
      "source": [
        "# Solution goes here\n",
        "\n",
        "# Output the result\n",
        "optimal_r = result.x\n",
        "print(f\"Optimal value of r: {optimal_r:.4f}\")\n",
        "print(f\"Minimum error achieved: {result.fun}\")\n",
        "\n",
        "# Optionally, run the simulation with the optimal r and plot again to verify\n",
        "milk.r = optimal_r\n",
        "timesteps, temperatures = run_simulation(milk, change_func)\n",
        "plt.figure(figsize=(8, 4))\n",
        "plt.plot(timesteps, temperatures, label=f'Temp of Milk with r={optimal_r:.4f}')\n",
        "plt.xlabel('Time (minutes)')\n",
        "plt.ylabel('Temperature (°C)')\n",
        "plt.title('Cooling of Milk with Optimized r')\n",
        "plt.legend()\n",
        "plt.grid(True)\n",
        "plt.show()"
      ]
    }
  ],
  "metadata": {
    "celltoolbar": "Tags",
    "kernelspec": {
      "display_name": "Python 3 (ipykernel)",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.10.6"
    },
    "colab": {
      "provenance": []
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}