coffee cooling sim.py

#-------------------------------------------------------------------------------
# Purpose: Coffee cooling simulation
# Author:      Pat (Hecuba)
#-------------------------------------------------------------------------------
#!/usr/bin/env python
from pprint import pprint, pformat

## Coffee cooling
# Start at 90 C, target temp is 75 C
# Cream drops temp 5 C immediately
# At what time do we add cream to minimize wait time?
##    Questions
##    How do you know that you have made dt small enough?
##    How confident are you in your value(s?) of c? Why?
##    What are the units of c?
##    How well does the mathematical model (Newton's law of cooling) work in reality?
##    Is it faster to add the cream before, or after?


class obj(object):
    pass

wiki_data = obj()
wiki_data.times = range(0, 47, 2)
wiki_data.black = [82.3, 78.5, 74.3, 70.7, 67.6, 65, 62.5, 60.1, 58.1, 56.1, 54.3, 52.8, 51.2, 49.9, 48.6, 47.2, 46.1, 45, 43.9, 43, 41.9, 41, 40.1, 39.5]
wiki_data.cream = [68.8, 64.8, 62.1, 59.9, 57.7, 55.9, 53.9, 52.3, 50.8, 49.5, 48.1, 46.8, 45.9, 44.8, 43.7, 42.6, 41.7, 40.8, 39.9, 39.3, 38.6, 37.7, 37, 36.4]

optimal_temp = 75
c = 0.03
dt = 1  # min
Tr = 23  # room temp

from pylab import *

# Group Hercules' code for plotting
def plot_coffee(temp_list, time_list, cream_insert_time=None, plot_title=None):
    scatter(time_list, temp_list, label="coffee temp (Celsius)")
    ylabel(r'Coffee temp (Celsius)')
    xlabel(r'Time (minutes)')
    plot(time_list, [optimal_temp] * len(time_list))
    if cream_insert_time != None:
        title('Coffee cooling with cream insertion at time {}, c={}'.format(cream_insert_time, c))
    elif plot_title:
        title(plot_title)
    black_times = array(wiki_data.times) + 4  # offset to match sim
    cream_times = array(wiki_data.times) + 6  # offset to match sim
    scatter(black_times, wiki_data.black, c='k', label='black coffee wiki data')
    scatter(cream_times, wiki_data.cream, c='m', label='cream coffee wiki data')
    legend(loc='upper right')
    show()

time_lists = []
Tcs_per_cream = []  # coffee temp
cream_insertion_times = [0, 5.5, 999] # range(1, 30, 1)  # minutes
for cream_time in cream_insertion_times:
    t_elapsed = 0
    time_list = [0]
    cream_added = False
    Tc = [90]
    # KLUDGE for adding cream at beginning
    if cream_time == 0:
        Tc[0] -= 5
        cream_added = True
    while Tc[-1] > 36:  # last temp in wiki data
        t_elapsed += dt
        time_list.append(t_elapsed)
        Tnext = -c * (Tc[-1] - Tr) * dt + Tc[-1]
        if t_elapsed >= cream_time and not cream_added:
            Tnext -= 5
            cream_added = True
        Tc.append(Tnext)
    time_lists.append(time_list)
    Tcs_per_cream.append(Tc)

for cream,Tcs,times in zip(cream_insertion_times, Tcs_per_cream, time_lists):
##    print "cream added at {} minutes".format(cream)
##    pprint(Tcs)
    plot_coffee(Tcs, times, cream)
##    break