<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Chapter 15
Modeling and Simulation in Python
Copyright 2021 Allen Downey
License
Step1:  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.
The coffee cooling problem
The coffee cooling problem was discussed by Jearl Walker in 
"The Amateur Scientist", Scientific American, Volume 237, Issue 5, November 1977. Since then it has become a standard example of modeling and simulation.
Here is my version of the problem
Step2:  The values of T_init, volume, and t_end come from the statement of the problem
Step3:  I chose the value of r arbitrarily for now; we will figure out how to estimate it soon.
In addition, make_system sets the temperature of the environment, T_env, and the time step, dt, which we will use use to simulate the cooling process.
Strictly speaking, Newton's law is a differential equation, but over a short period of time we can approximate it with a difference equation
Step4:  We can test it with the initial temperature of the coffee, like this
Step5:  With dt=1 minute, the temperature drops by about 0.7 °C/min, at least for this value of r.
Now here's a version of run_simulation that simulates a series of time steps from t_0 to t_end
Step6:  This function is similar to previous versions of run_simulation.
One difference is that it uses linrange to make an array of values
from t_0 to t_end with time step dt.
We can run it like this
Step7:  The result is a TimeSeries with one row per time step.
Step8:  Here's what the results look like.
Step9:  The temperature after 30 minutes is 72.3 °C, which is a little higher than what's stated in the problem, 70 °C.
Step10:  By trial and error, we could find the value of r where the final temperature is precisely 70 °C.
But it is more efficient to use a root-finding algorithm.
Finding roots
ModSim provides a function called root_scalar that finds the roots of non-linear equations. As a simple example, suppose you want to find the roots of the polynomial 
$$f(x) = (x - 1)(x - 2)(x - 3)$$ 
A root is a value of $x$ that makes $f(x)=0$. Because of the way I wrote the 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.
I'll use this example to demonstrate root_scalar. First, we have to
write a function that evaluates $f$
Step11:  Now we call root_scalar like this
Step12:  The first argument is the function whose roots we want. The second
argument is an interval that contains or "brackets" a root. The result is an object that contains several variables, including root, which stores the root that was found.
Step13:  If we provide a different interval, we find a different root.
Step14:  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".
res = root_scalar(func, bracket=[4, 5])
Now we can use root_scalar to estimate r.
Estimating r
What we want is the value of r that yields a final temperature of
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
Step15:  This is called an "error function" because it returns the
difference between what we got and what we wanted, that is, the error.
With the right value of r, the error is 0.
We can test error_func like this, using the initial guess r=0.01
Step16:  The result is an error of 2.3 °C, because the final temperature with
this value of r is too high.
Step17:  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.
So we can call root_scalar like this
Step18:  The first argument is the error function.
The second argument is the System object, which root_scalar passes as an argument to error_func.
The third argument is an interval that brackets the root.
Here are the results.
Step19:  In this example, r_coffee turns out to be about 0.0115, in units of min$^{-1}$ (inverse minutes).
We can confirm that this value is correct by setting r to the root we found and running the simulation.
Step20:  The final temperature is very close to 70 °C.
Exercises
Exercise
Step21:  Exercise
<ASSISTANT_TASK:>
Python Code:
# install Pint if necessary
try:
    import pint
except ImportError:
    !pip install pint
# download modsim.py if necessary
from os.path import exists
filename = 'modsim.py'
if not exists(filename):
    from urllib.request import urlretrieve
    url = 'https://raw.githubusercontent.com/AllenDowney/ModSim/main/'
    local, _ = urlretrieve(url+filename, filename)
    print('Downloaded ' + local)
# import functions from modsim
from modsim import *
Explanation: Chapter 15
Modeling and Simulation in Python
Copyright 2021 Allen Downey
License: Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International
End of explanation
def make_system(T_init, volume, r, t_end):
    return System(T_init=T_init,
                  T_final=T_init,
                  volume=volume,
                  r=r,
                  t_end=t_end,
                  T_env=22,
                  t_0=0,
                  dt=1)
Explanation: 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.
The coffee cooling problem
The coffee cooling problem was discussed by Jearl Walker in 
"The Amateur Scientist", Scientific American, Volume 237, Issue 5, November 1977. Since then it has become a standard example of modeling and simulation.
Here is my version of the problem:
Suppose I stop on the way to work to pick up a cup of coffee, which I take with 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?
To help answer this question, I made a trial run with the milk and
coffee in separate containers and took some measurements (not really):
When served, the temperature of the coffee is 90 °C. The volume is
    300 mL.
The milk is at an initial temperature of 5 °C, and I take about
    50 mL.
The ambient temperature in my car is 22 °C.
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.
The milk container is not as well insulated. After 15 minutes, it
    warms up to 20 °C, nearly the ambient temperature.
To use this data and answer the question, we have to know something
about temperature and heat, and we have to make some modeling decisions.
Temperature and heat
To understand how coffee cools (and milk warms), we need a model of
temperature and heat. Temperature is a property of an object or a
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 up the object.
When particles in a hot object contact particles in a cold object, the
hot object gets cooler and the cold object gets warmer as energy is
transferred from one to the other. The transferred energy is called
heat; in SI units it is measured in joules (J).
Heat is related to temperature by the following equation (see
http://modsimpy.com/thermass): 
$$Q = C~\Delta T$$ 
where $Q$ is the amount of heat transferred to an object, $\Delta T$ is resulting change in temperature, and $C$ is the thermal mass of the object, which quantifies how much energy it takes to heat or cool it. In SI units, thermal mass is measured in joules per degree Celsius (J/°C).
For objects made primarily from one material, thermal mass can be
computed like this: 
$$C = m c_p$$ 
where $m$ is the mass of the object and $c_p$ is the specific heat capacity of the material (see http://modsimpy.com/specheat).
We can use these equations to estimate the thermal mass of a cup of
coffee. The specific heat capacity of coffee is probably close to that
of water, which is 4.2 J/g/°C. Assuming that the density of coffee is
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.
So when a cup of coffee cools from 90 °C to 70 °C, the change in
temperature, $\Delta T$ is 20 °C, which means that 25 200 J of heat
energy was transferred from the coffee to the surrounding environment
(the cup holder and air in my car).
To give you a sense of how much energy that is, if you were able to
harness all of that heat to do work (which you cannot), you could
use it to lift a cup of coffee from sea level to 8571 m, just shy of the height of Mount Everest, 8848 m.
Assuming that the cup has less mass than the coffee, and is made from a material with lower specific heat, we can ignore the thermal mass of the cup. 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.
Heat transfer
In a situation like the coffee cooling problem, there are three ways
heat transfers from one object to another (see http://modsimpy.com/transfer):
Conduction: When objects at different temperatures come into
    contact, the faster-moving particles of the higher-temperature
    object transfer kinetic energy to the slower-moving particles of the lower-temperature object.
Convection: When particles in a gas or liquid flow from place to
    place, they carry heat energy with them. Fluid flows can be caused
    by external action, like stirring, or by internal differences in
    temperature. For example, you might have heard that hot air rises,
    which is a form of "natural convection\".
Radiation: As the particles in an object move due to thermal energy,
    they emit electromagnetic radiation. The energy carried by this
    radiation depends on the object's temperature and surface properties
    (see http://modsimpy.com/thermrad).
For objects like coffee in a car, the effect of radiation is much
smaller than the effects of conduction and convection, so we will ignore it.
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".
Newton's law of cooling
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:
$$\frac{dT}{dt} = -r (T - T_{env})$$ 
where $T$, the temperature of the object, is a function of time, $t$; $T_{env}$ is the temperature of the environment, and $r$ is a constant that characterizes how quickly heat is transferred between the system and the environment.
Newton's so-called "law " is really a model: it is a good approximation in some conditions and less good in others.
For example, if the primary mechanism of heat transfer is conduction,
Newton's law is "true", which is to say that $r$ is constant over a
wide range of temperatures. And sometimes we can estimate $r$ based on
the material properties and shape of the object.
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.
When radiation makes up a substantial part of heat transfer, Newton's
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
hundred degrees Celsius, say).
However, for a situation like the coffee cooling problem, we expect
Newton's model to be quite good.
Implementation
To get started, let's forget about the milk temporarily and focus on the coffee. 
Here's a function that takes the parameters of the system and makes a System object:
End of explanation
coffee = make_system(T_init=90, volume=300, r=0.01, t_end=30)
Explanation: The values of T_init, volume, and t_end come from the statement of the problem:
End of explanation
def change_func(T, t, system):
    r, T_env, dt = system.r, system.T_env, system.dt    
    return -r * (T - T_env) * dt
Explanation: I chose the value of r arbitrarily for now; we will figure out how to estimate it soon.
In addition, make_system sets the temperature of the environment, T_env, and the time step, dt, which we will use use to simulate the cooling process.
Strictly speaking, Newton's law is a differential equation, but over a short period of time we can approximate it with a difference equation:
$$\Delta T = -r (T - T_{env}) dt$$ 
where $dt$ is the time step and $\Delta T$ is the change in temperature during that time step.
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
environment, $T - T_{env}$. To minimize confusion, I avoid this second
use.
The following function takes the current temperature, T, the current time t, and a System object, and computes the change in temperature during a time step:
End of explanation
change_func(coffee.T_init, 0, coffee)
Explanation: We can test it with the initial temperature of the coffee, like this:
End of explanation
def run_simulation(system, change_func):
    t_array = linrange(system.t_0, system.t_end, system.dt)
    n = len(t_array)
    
    series = TimeSeries(index=t_array)
    series.iloc[0] = system.T_init
    
    for i in range(n-1):
        t = t_array[i]
        T = series.iloc[i]
        series.iloc[i+1] = T + change_func(T, t, system)
    
    system.t_end = t_array[-1]
    system.T_final = series.iloc[-1]
    return series
Explanation: With dt=1 minute, the temperature drops by about 0.7 °C/min, at least for this value of r.
Now here's a version of run_simulation that simulates a series of time steps from t_0 to t_end:
End of explanation
results = run_simulation(coffee, change_func)
Explanation: This function is similar to previous versions of run_simulation.
One difference is that it uses linrange to make an array of values
from t_0 to t_end with time step dt.
We can run it like this:
End of explanation
results.head()
Explanation: The result is a TimeSeries with one row per time step.
End of explanation
results.plot(label='coffee')
decorate(xlabel='Time (minute)',
         ylabel='Temperature (C)',
         title='Coffee Cooling')
Explanation: Here's what the results look like.
End of explanation
coffee.T_final
Explanation: The temperature after 30 minutes is 72.3 °C, which is a little higher than what's stated in the problem, 70 °C.
End of explanation
def func(x):
    return (x-1) * (x-2) * (x-3)
Explanation: By trial and error, we could find the value of r where the final temperature is precisely 70 °C.
But it is more efficient to use a root-finding algorithm.
Finding roots
ModSim provides a function called root_scalar that finds the roots of non-linear equations. As a simple example, suppose you want to find the roots of the polynomial 
$$f(x) = (x - 1)(x - 2)(x - 3)$$ 
A root is a value of $x$ that makes $f(x)=0$. Because of the way I wrote the 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.
I'll use this example to demonstrate root_scalar. First, we have to
write a function that evaluates $f$:
End of explanation
res = root_scalar(func, bracket=[1.5, 2.5])
res
Explanation: Now we call root_scalar like this:
End of explanation
res.root
Explanation: The first argument is the function whose roots we want. The second
argument is an interval that contains or "brackets" a root. The result is an object that contains several variables, including root, which stores the root that was found.
End of explanation
res = root_scalar(func, bracket=[2.5, 3.5])
res.root
Explanation: If we provide a different interval, we find a different root.
End of explanation
def error_func(r, system):
    system.r = r
    results = run_simulation(system, change_func)
    return system.T_final - 70
Explanation: 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".
res = root_scalar(func, bracket=[4, 5])
Now we can use root_scalar to estimate r.
Estimating r
What we want is the value of r that yields a final temperature of
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:
End of explanation
coffee = make_system(T_init=90, volume=300, r=0.01, t_end=30)
error_func(0.01, coffee)
Explanation: This is called an "error function" because it returns the
difference between what we got and what we wanted, that is, the error.
With the right value of r, the error is 0.
We can test error_func like this, using the initial guess r=0.01:
End of explanation
error_func(0.02, coffee)
Explanation: The result is an error of 2.3 °C, because the final temperature with
this value of r is too high.
End of explanation
res = root_scalar(error_func, coffee, bracket=[0.01, 0.02])
res.flag
Explanation: 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.
So we can call root_scalar like this:
End of explanation
r_coffee = res.root
r_coffee
Explanation: The first argument is the error function.
The second argument is the System object, which root_scalar passes as an argument to error_func.
The third argument is an interval that brackets the root.
Here are the results.
End of explanation
coffee.r = res.root
run_simulation(coffee, change_func)
coffee.T_final
Explanation: In this example, r_coffee turns out to be about 0.0115, in units of min$^{-1}$ (inverse minutes).
We can confirm that this value is correct by setting r to the root we found and running the simulation.
End of explanation
# Solution
milk = make_system(T_init=5, t_end=15, r=0.1, volume=50)
results_milk = run_simulation(milk, change_func)
milk.T_final
# Solution
results_milk.plot(color='C1', label='milk')
decorate(xlabel='Time (minutes)',
         ylabel='Temperature (C)')
Explanation: The final temperature is very close to 70 °C.
Exercises
Exercise:  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.
By trial and error, find a value for r that makes the final temperature close to 20 °C.
End of explanation
# Solution
def error_func2(r, system):
    system.r = r
    results = run_simulation(system, change_func)
    return system.T_final - 20
# Solution
root_scalar(error_func2, milk, bracket=[0.1, 0.2])
res
# Solution
run_simulation(milk, change_func)
milk.T_final
Explanation: Exercise: 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.
End of explanation
<END_TASK>