Skip to content

Instantly share code, notes, and snippets.

@pckujawa
Created February 1, 2013 16:20
Show Gist options
  • Save pckujawa/4692310 to your computer and use it in GitHub Desktop.
Save pckujawa/4692310 to your computer and use it in GitHub Desktop.
#-------------------------------------------------------------------------------
# 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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment