Laboratory 10: Cup of tea

Prerequisits: curve fitting, data analysis, inverting function (root find)

scipy.optimize.curve_fit

Background

Enjoying a good cup of tea or coffee is an important part of our academic activities. Here, we study how a hot cup of tea or coffee cools down once the liquid has been poured into the cup.

A dedicated research unit in a little known university in central Java, Indonesia, has conducted temperature measurements for a cup of tea and made the data available:

The specialised equipment takes one temperature reading (in degree Celsius) every 10 seconds, although the measurements are somewhat noisy. Unfortunately, something went wrong in the beginning of the measurement, so the data for the first minute are missing.

Assessed part

Create a file lab10.py and include a copy of the model() function definition from the training part (Do not use from training10 import model -- although this is an elegant solution in general, it will not work when we test your submission automatically).

Then add the following functions to lab10.py:

1. read2coldata(filename) which opens a text file with two columns of data. The columns have to be separated by white space. The function should return a tuple of two numpy-arrays where the first contains all the data from the first column in the file, and the second all the data in the second column.

   Example: for a data file testdata.txt which contains

   1.5  4
   8    5
   16   6
   17   6.2
   
   

   we expect this behaviour:

   In [ ]: read2coldata('testdata.txt')
   Out[ ]: (array([  1.5,   8. ,  16. ,  17. ]), array([ 4. ,  5. ,  6. ,  6.2]))
   
   In [ ]: a, b = read2coldata('testdata.txt')
   
   In [ ]: a
   Out[ ]: array([  1.5,   8. ,  16. ,  17. ])
   
   In [ ]: b
   Out[ ]: array([ 4. ,  5. ,  6. ,  6.2])
   

2. A function extract_parameters(ts, Ts) which expects a numpy array ts with time values and a numpy array Ts of the same length as ts with corresponding temperature values. The function should estimate and return a tuple of the three model parameters Ti, Ta and c (in this order) by fitting the model function as in equation (1) to the data ts and Ts.

   Hints:

   • The curve_fit function may need some initial guesses (through the optional function parameter p0) for the model parameters to be able to find a good fit.

   • You are strongly encouraged to plot your fitted curve together with the raw data to check that the fit is reasonable.

     You can use a function like this:

     def plot(ts, Ts, Ti, Ta, c):
         """Input Parameters:
     
           ts : Data for time (ts)
                     (numpy array)
           Ts : data for temperature (Ts)
                     (numpy arrays)
           Ti : model parameter Ti for Initial Temperature
                     (number)
           Ta : model parameter Ta for Ambient Temperature
                     (number)
           c  : model parameter c for the time constant
                     (number)
     
         This function will create plot that shows the model fit together
         with the data.
     
         Function returns None.
         """
     
         pylab.plot(ts, Ts, 'o', label='data')
         fTs = model(ts, Ti, Ta, c)
         pylab.plot(ts, fTs, label='fitted')
         pylab.legend()
         pylab.savefig('testcompare.pdf')  # or pylab.show()
     

3. A function sixty_degree_time(Ti, Ta, c) which expects the model paramaters Ti (initial temperature), Ta (ambient temperature) and c the cooling rate time constant. The function should return an estimate of the number of seconds after which the temperature of the drink has cooled down to 60 degree Celsius (60 degree is a temperature that is generally considered low enough not to damage tissue).

   You have at least two different possible ways of obtaining this number of seconds for a given set of model parameters (Ti, Ta, c). One involves a root finding algorithm. You can assume that Ti > 60 degree Celsius, and that Ti > Ta.

To test your functions, you can put them together like this in the ipython or in a separate file:

from lab10 import read2coldata, extract_parameters, sixty_degree_time

ts, Ts = read2coldata('time_temp.txt')
Ti, Ta, c = extract_parameters(ts, Ts)
print("Model parameters are Ti={} C, Tf={}C,".format(Ti, Ta))
print("                     time constant={}s".format(c))
waittime = sixty_degree_time(Ti, Ta, c)
print("The drink reaches 60 degrees after {:.2f} seconds = {:.2f} minutes"
      .format(waittime, waittime / 60.))

If you want to include the testing in your lab10.py file, you need to put it inside the if __name__ == "___main__" statement, such as

if __name__ == "__main__":
    ts, Ts = read2coldata('time_temp.txt')
    Ti, Ta, c = extract_parameters(ts, Ts)
    print("Model parameters are Ti={} C, Tf={}C,".format(Ti, Ta))
    print("                     time constant={}s".format(c))
    waittime = sixty_degree_time(Ti, Ta, c)
    print("The drink reaches 60 degrees after {:.2f} seconds = {:.2f} minutes"
      .format(waittime, waittime / 60.))

Written like this, the test will only be executed if you execute lab10.py as the main file (for example by pressing F5 when editing lab10.py in Spyder), but not if the file is imported from elsewhere (which is what the testing system does).

Can you now answer the three research questions, i.e. (No need to submit these answers.)

• what was the initial temperature of the tea in the cup at time t=0s?
• how quickly does the tea cool down? In particular: after what time is it safe to drink it (we assume that 60C are a safe temperature).
• what will be the final temperature of the tea if we wait infinitely long (presumably this will be the room temperature in this particular lab in Java).

When you have completed your work, check and test the code carefully (including documentation strings). Submit your file by email with subject line lab 10.