{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling and Simulation in Python\n",
    "\n",
    "Chapter 7: Thermal systems\n",
    "\n",
    "Copyright 2017 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# If you want the figures to appear in the notebook, \n",
    "# and you want to interact with them, use\n",
    "# %matplotlib notebook\n",
    "\n",
    "# If you want the figures to appear in the notebook, \n",
    "# and you don't want to interact with them, use\n",
    "# %matplotlib inline\n",
    "\n",
    "# If you want the figures to appear in separate windows, use\n",
    "# %matplotlib qt5\n",
    "\n",
    "# tempo switch from one to another, you have to select Kernel->Restart\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The coffee cooling problem.\n",
    "\n",
    "I'll use a `State` object to store the initial temperature.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "init = State(temp=90)\n",
    "init"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And a `System` object to contain the system parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "coffee = System(init=init,\n",
    "                volume=300,\n",
    "                r=0.01,\n",
    "                T_env=22,\n",
    "                t0=0, \n",
    "                t_end=30,\n",
    "                dt=1)\n",
    "coffee"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `update` function implements Newton's law of cooling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def update(state, system):\n",
    "    \"\"\"Update the thermal transfer model.\n",
    "    \n",
    "    state: State (temp)\n",
    "    system: System object\n",
    "    \n",
    "    returns: State (temp)\n",
    "    \"\"\"\n",
    "    unpack(system)\n",
    "    T = state.temp\n",
    "    T += -r * (T - T_env) * dt\n",
    "\n",
    "    return State(temp=T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's how it works."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "update(init, coffee)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can run simulations using the same function from the previous chapter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def run_simulation(system, update_func):\n",
    "    \"\"\"Runs a simulation of the system.\n",
    "    \n",
    "    Add a TimeFrame to the System: results\n",
    "    \n",
    "    system: System object\n",
    "    update_func: function that updates state\n",
    "    \"\"\"\n",
    "    unpack(system)\n",
    "    \n",
    "    frame = TimeFrame(columns=init.index)\n",
    "    frame.loc[t0] = init\n",
    "    ts = linrange(t0, t_end-dt, dt)\n",
    "    \n",
    "    for t in ts:\n",
    "        frame.loc[t+dt] = update_func(frame.loc[t], system)\n",
    "    \n",
    "    system.results = frame"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's how it works."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "run_simulation(coffee, update)\n",
    "coffee.results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what the results look like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(coffee.results.temp, label='coffee')\n",
    "decorate(xlabel='Time (minutes)',\n",
    "         ylabel='Temperature (C)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After running the simulation, we can extract the final temperature from the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def final_temp(system):\n",
    "    \"\"\"Final temperature.\n",
    "    \n",
    "    If system has no results, return initial temp.\n",
    "    \n",
    "    system: System object.\n",
    "    \n",
    "    returns: temperature (degC)\n",
    "    \"\"\"    \n",
    "    if hasattr(system, 'results'):\n",
    "        return system.results.temp[system.t_end]\n",
    "    else:\n",
    "        return system.init.temp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It will be convenient to wrap these steps in a function.  `kwargs` is a collection of whatever keyword arguments are provided; they are passed along as arguments to `System`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def make_system(T_init=90, r=0.01, volume=300, t_end=30):\n",
    "    \"\"\"Runs a simulation with the given parameters.\n",
    "\n",
    "    T_init: initial temperature in degC\n",
    "    r: heat transfer rate, in 1/min\n",
    "    volume: volume of liquid in mL\n",
    "    t_end: end time of simulation\n",
    "    \n",
    "    returns: System object\n",
    "    \"\"\"\n",
    "    init = State(temp=T_init)\n",
    "    \n",
    "    system = System(init=init,\n",
    "                    volume=volume,\n",
    "                    r=r,\n",
    "                    T_env=22, \n",
    "                    t0=0,\n",
    "                    t_end=t_end,\n",
    "                    dt=1)\n",
    "    return system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's how we use it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "coffee = make_system(T_init=90, r=0.01, volume=300, t_end=30)\n",
    "run_simulation(coffee, update)\n",
    "final_temp(coffee)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:**  Simulate the temperature of 50 mL of milk with a starting temperature of 5 degC, in a vessel with the same insulation, for 15 minutes, and plot the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(milk.results.temp, label='milk')\n",
    "decorate(xlabel='Time (minutes)',\n",
    "         ylabel='Temperature (C)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using `fsolve`\n",
    "\n",
    "As a simple example, let's find the roots of this function; that is, the values of `x` that make the result 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def func(x):\n",
    "    return (x-1) * (x-2) * (x-3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`modsim.py` provides `fsolve`, which does some error-checking and then runs `scipy.optimize.fsolve`.  The first argument is the function whose roots we want.  The second argument is an initial guess."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "fsolve(func, x0=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Usually the root we get is the one that's closest to the initial guess."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "fsolve(func, 1.9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "fsolve(func, 2.9)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But not always."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "fsolve(func, 1.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to find the value of `r` that makes the final temperature 70, so we define an \"error function\" that takes `r` as a parameter and returns the difference between the final temperature and the goal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def error_func1(r):\n",
    "    \"\"\"Runs a simulation and returns the `error`.\n",
    "    \n",
    "    r: heat transfer rate, in 1/min\n",
    "    \n",
    "    returns: difference between final temp and 70 C\n",
    "    \"\"\"\n",
    "    system = make_system(T_init=90, r=r, volume=300, t_end=30)\n",
    "    run_simulation(system, update)\n",
    "    return final_temp(system) - 70"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With `r=0.01`, we end up a little too warm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "error_func1(r=0.01)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The return value from `fsolve` is an array with a single element, the estimated value of `r`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "solution = fsolve(error_func1, 0.01)\n",
    "r_coffee = solution[0]\n",
    "r_coffee"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we run the simulation with the estimated value of `r`, the final temperature is 70 C, as expected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "coffee = make_system(T_init=90, r=r_coffee, volume=300, t_end=30)\n",
    "run_simulation(coffee, update)\n",
    "final_temp(coffee)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:**  When you call `fsolve`, it calls `error_func1` several times.  To see how this works, add a print statement to `error_func1` and run `fsolve` again."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Repeat this process to estimate `r_milk`, given that it starts at 5 C and reaches 20 C after 15 minutes.  \n",
    "\n",
    "Before you use `fsolve`, you might want to try a few values for `r_milk` and see how close you can get by trial and error.  Here's an initial guess to get you started:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "r_milk = 0.1\n",
    "milk = make_system(T_init=5, t_end=15, r=r_milk, volume=50)\n",
    "run_simulation(milk, update)\n",
    "final_temp(milk)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mixing liquids"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function takes `System` objects that represent two liquids, computes the temperature of the mixture, and returns a new `System` object that represents the mixture."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def mix(s1, s2):\n",
    "    \"\"\"Simulates the mixture of two liquids.\n",
    "    \n",
    "    s1: System representing coffee\n",
    "    s2: System representing milk\n",
    "    \n",
    "    returns: System representing the mixture\n",
    "    \"\"\"\n",
    "    assert s1.t_end == s2.t_end\n",
    "    \n",
    "    volume = s1.volume + s2.volume\n",
    "    \n",
    "    temp = (s1.volume * final_temp(s1) + \n",
    "            s2.volume * final_temp(s2)) / volume\n",
    "    \n",
    "    mixture = make_system(T_init=temp,\n",
    "                          volume=volume,\n",
    "                          r=s1.r)\n",
    "    \n",
    "    return mixture"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we'll see what happens if we add the milk at the end.  We'll simulate the coffee and the milk separately."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "coffee = make_system(T_init=90, t_end=30, r=r_coffee, volume=300)\n",
    "run_simulation(coffee, update)\n",
    "final_temp(coffee)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "milk = make_system(T_init=5, t_end=30, r=r_milk, volume=50)\n",
    "run_simulation(milk, update)\n",
    "final_temp(milk)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what the results look like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(coffee.results.temp, label='coffee')\n",
    "plot(milk.results.temp, '--', label='milk')\n",
    "decorate(xlabel='Time (minutes)',\n",
    "         ylabel='Temperature (C)',\n",
    "         loc='center left')\n",
    "\n",
    "savefig('chap07-fig01.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what happens when we mix them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "mix_last = mix(coffee, milk)\n",
    "final_temp(mix_last)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's what we get if we add the milk immediately."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "coffee = make_system(T_init=90, r=r_coffee, volume=300)\n",
    "milk = make_system(T_init=5, r=r_milk, volume=50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "mix_first = mix(coffee, milk)\n",
    "mix_first.t_end = 30\n",
    "run_simulation(mix_first, update)\n",
    "final_temp(mix_first)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function takes `t_add`, which is the time when the milk is added, and returns the final temperature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def run_and_mix(t_add, t_total=30):\n",
    "    \"\"\"Simulates two liquids and them mixes them at t_add.\n",
    "    \n",
    "    t_add: time in minutes\n",
    "    t_total: total time to simulate, min\n",
    "    \n",
    "    returns: final temperature\n",
    "    \"\"\"\n",
    "    coffee = make_system(T_init=90, t_end=t_add, \n",
    "                         r=r_coffee, volume=300)\n",
    "    run_simulation(coffee, update)\n",
    "\n",
    "    milk = make_system(T_init=5, t_end=t_add, \n",
    "                       r=r_milk, volume=50)\n",
    "    run_simulation(milk, update)\n",
    "    \n",
    "    mixture = mix(coffee, milk)\n",
    "    mixture.t_end = t_total - t_add\n",
    "    run_simulation(mixture, update)\n",
    "\n",
    "    return final_temp(mixture)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can try it out with a few values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "run_and_mix(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "run_and_mix(15)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "run_and_mix(30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And then sweep a range of values for `t_add`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sweep = SweepSeries()\n",
    "for t_add in linrange(0, 30, 2):\n",
    "    temp = run_and_mix(t_add)\n",
    "    sweep[t_add] = temp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what the result looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(sweep, color='purple')\n",
    "decorate(xlabel='Time added (min)',\n",
    "         ylabel='Final temperature (C)',\n",
    "         legend=False)\n",
    "\n",
    "savefig('chap07-fig02.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:**  Suppose the coffee shop won't let me take milk in a separate container, but I keep a bottle of milk in the refrigerator at my office.  In that case is it better to add the milk at the coffee shop, or wait until I get to the office?\n",
    "\n",
    "Hint: Think about the simplest way to represent the behavior of a refrigerator in this model.  The change you make to test this variation of the problem should be very small!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use the analytic result to compute temperature as a function of time.  The following function is similar to `run_simulation`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def run_analysis(system):\n",
    "    \"\"\"Computes temperature using the analytic solution.\n",
    "    \n",
    "    Adds TimeFrame to `system` as `results`\n",
    "    \n",
    "    system: System object\n",
    "    \"\"\"\n",
    "    unpack(system)\n",
    "    \n",
    "    T_init = init.temp    \n",
    "    ts = linrange(t0, t_end, dt)\n",
    "    \n",
    "    temp_array = T_env + (T_init - T_env) * exp(-r * ts)\n",
    "    temp_series = TimeSeries(temp_array, index=ts)\n",
    "    \n",
    "    system.results = TimeFrame(temp_series, columns=['temp'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's how we run it.  From the analysis, we have the computed value of `r_coffee2`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "r_coffee2 = 0.011610223142273859"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "init = State(temp=90)\n",
    "coffee2 = System(init=init, T_env=22, r=r_coffee2, \n",
    "                 t0=0, t_end=30)\n",
    "run_analysis(coffee2)\n",
    "final_temp(coffee2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can compare to the results from simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "init = State(temp=90)\n",
    "coffee = System(init=init, T_env=22, r=r_coffee, \n",
    "                t0=0, t_end=30, dt=1)\n",
    "run_simulation(coffee, update)\n",
    "final_temp(coffee)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "They are identical except for small roundoff errors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "coffee.results - coffee2.results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}