# We start with a cup of very hot black coffee.
# If we want to drink it as soon as possible,
# when should we add the cold milk - sooner or later?
#
# this program lets you experiment, using Newton's cooling law
# 2/2004, Wolfgang.Urban@schule.at

# good old Runge Kutta 4(4)
def rk4(t,x,f,h):
    k1 = h*f(t,x)
    k2 = h*f(t+h/2.0,x+k1/2.0)
    k3 = h*f(t+h/2.0,x+k2/2.0)
    k4 = h*f(t+h,x+k3)
    tnew = t+h
    xnew = x+(k1+2.0*k2+2.0*k3+k4)/6.0
    return tnew,xnew


# we model our cup of coffee
class CupOfCoffee:
    def __init__(self,amount=100,temp=90,roomtemp=20):
        self.c_amount = float(amount)
        self.c_temp = float(temp)
        self.roomtemp = float(roomtemp)
        self.cooling_const = 0.5            # found by experiment
        self.time = 0

    # mixing temperature and new amount
    def add_milk(self,amount=10,temp=5):
        self.c_temp = float(self.c_amount*self.c_temp+amount*temp)/   \
                      float(self.c_amount+amount)
        self.c_amount += float(amount)


    # how long does it take to reach this temperature
    def wait_for_temp(self,temp):
        if temp>self.c_temp:
            raise "coffee is already cooler!"
        if temp<self.roomtemp:
            raise "room is warmer than coffee!"

        h = 0.1
        while 1:
            if self.c_temp<=temp: return self.time
            self.time,self.c_temp = rk4(self.time,self.c_temp,self.newton_cool,h)
            

    # return coffee's temperature
    def gettemp(self): return self.c_temp

    # let time pass for some seconds
    def wait(self,seconds):
        h = 0.1
        for tt in range(int(seconds*10)):
            self.time,self.c_temp = rk4(self.time,self.c_temp,self.newton_cool,h)

    # Newton's cooling Law T' = c(T-T0)
    def newton_cool(self,t,x):
        return -self.cooling_const/self.c_amount*(x-self.roomtemp)


# demo function

def coffee_with_milk(add_milk_after=60):
    
    # 100 ml coffee of 90 degrees in a room with 20 degrees
    c = CupOfCoffee(100,90,20)

    # wait for 60 seconds
    c.wait(add_milk_after)

    ## print c.time,c.gettemp(), "degrees before milk"

    # 10 ml milk of 5 degrees added
    c.add_milk(10,5)

    ## print c.time,c.gettemp(), "degrees after milk"

    # when will this mixture have 30 degrees?
    print c.wait_for_temp(30),"seconds for",c.gettemp(),"degrees"