#!/usr/bin/python
# -*- coding: utf-8 -*-

from scipy import *
from scipy.integrate import odeint

#Set the constants. 
D0 = 1.07173*10**12
S0= 5.35866*10**10
kb = 8.6173324*10**(-5)
K = 1./21
tlo = 20
thi = 30
#Formating string for printing the numbers
formatStr = "%8.3f"

#Death rate as a function of temperature
def D(T):
    return D0 * exp(-0.66/(kb * (T + 273.15)))
#Sexy time as a function of temperature
def S(T):
    return S0 * exp(-0.66/(kb * (T + 273.15)))

#Calculate the derivative
def derivative(species, time, temp):
    crow_der = S(temp) * species[0] - K * species[0] * species[1]
    cat_der = K * species[0] * species[1] - D(temp) * species[1]
    return (crow_der, cat_der)

#Set the time intervals and initial values
t = arange(0,16,1*10**(-2))
species_i = array([110,7])

#Go through the temperatures
temps = range(tlo, thi + 1)

print temps
result = empty( (len(t),2,len(temps)) )

#For each temperature, solve the ODE
for i,temp in zip(range(len(temps)), temps):
    result[:,:,i] = odeint(derivative, species_i, t, args=(temp,), printmessg=True)

#Print them out
lines = []
for i in range(len(t)):
    lines.append(formatStr % t[i])
    for j in range(len(temps)):
        lines.append(" ")
        lines.append(formatStr % result[i,0,j])
    for j in range(len(temps)):
        lines.append(" ")
        lines.append(formatStr % result[i,1,j])
    lines.append("\n")

with open("out.txt", "w") as f:
    f.write("".join(lines))