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

# In[1]:


import numpy as np
import matplotlib.pyplot as plt
get_ipython().run_line_magic('matplotlib', 'inline')
from scipy.integrate import odeint

k = 0.001 # try this or try k = 1 or try k = 2

from IPython.html.widgets import interact
from IPython.display import clear_output, display, HTML
def firstorder(k, C_A0):
    def f(C_list, time_array):
        C_A = C_list[0]
        C_B = C_list[1]
        dcadt = -k*C_A
        dcbdt = k*C_A
        return [dcadt, dcbdt]

#     k = 0.001
#     C_A0 = 10 
    C_B0 = 0 
    time_start = 0
    time_finish = 1000
    N_points = 1000
    time_array = np.linspace(time_start, time_finish, N_points)

    C_num_list = odeint(f, [C_A0, C_B0], time_array)

    C_A_num = C_num_list[:,0]
    C_B_num = C_num_list[:,1]

    # print(C_num_list)
    plt.plot(time_array, C_A_num, 'r--',label='A')
    plt.plot(time_array, C_B_num, 'b--',label='B')
    plt.xlabel('t, s')
    plt.ylabel('C, mol/L')
    plt.legend()
    plt.show()

interact(firstorder, k = (0.0, 0.1, 0.01), C_A0 = (1., 10., 0.5))


# In[19]:


import numpy as np
import matplotlib.pyplot as plt
get_ipython().run_line_magic('matplotlib', 'inline')
from scipy.integrate import odeint
plt.style.use('ggplot')

# from IPython.html.widgets import interact
# from IPython.display import clear_output, display, HTML
def systemodes(k1,k2):
    def f(C_list, time_array):
        C_A = C_list[0]
        C_B = C_list[1]
        C_C = C_list[2]
        C_D = C_list[3]
        r1 = k1*C_A*C_B
        r2 = k2*C_B*C_C
        dCadt = -r1
        dCbdt = -r1 -r2
        dCcdt = r1 - r2
        dCddt = r2
        return [dCadt, dCbdt, dCcdt, dCddt]

    time_start = 0
    time_finish = 3
    N_points = 100
    time_array = np.linspace(time_start, time_finish, N_points)
  
    C_0_list = [1, 1, 0, 0] # [C_A0, C_B0, C_C0, C_D0]
    C_num_list = odeint(f, C_0_list, time_array)
    
    
    C_A_num = C_num_list[:,0]
    C_B_num = C_num_list[:,1]
    C_C_num = C_num_list[:,2]
    C_D_num = C_num_list[:,3]

    # print(C_num_list)
    plt.plot(time_array, C_A_num, 'r--',label='A')
    plt.plot(time_array, C_B_num, 'b--',label='B')
    plt.plot(time_array, C_C_num, 'g--',label='C')
    plt.plot(time_array, C_D_num, 'k--',label='D')
    plt.xlabel('t, s')
    plt.ylabel('C, mol/L')
    plt.legend()
    plt.show()

    
interact(systemodes, k1=(0., 1, 0.1),k2=(0., 2, 0.1))
#tips and tricks of an old python enthusiast

# systemod es(k1=1., k2=2.1)
# cmd - ]
# cntrl - ]

# shift-enter -run cell
# autocomplete - tab
# doc shift-tab 
# comment stuff : cmd - /
# comment stuff control - /

# k1 = 1.
# k2 = 1.5


    


# In[2]:


import numpy as np
import matplotlib.pyplot as plt
get_ipython().run_line_magic('matplotlib', 'inline')
from scipy.integrate import odeint
plt.style.use('ggplot')


P = 0 # 10
# d = 0.0001
B = 0.0095
G = 0.001
A = 0.001

def russians_are_not_scared_of_zombies(d):

    def f(C_list, time_array):
        Si = C_list[0]
        Zi = C_list[1]
        Ri = C_list[2]
        dsdt = P - B*Si*Zi - d*Si
        dzdt = B*Si*Zi + G*Ri - A*Si*Zi
        drdt = d*Si + A*Si*Zi - G*Ri
        return [dsdt, dzdt, drdt]



    S0 = 500
    Z0 = 0
    R0 = 0
    C_list_0 = [S0, Z0, R0]

    time_start = 0
    time_finish = 10.
    N_points = 100
    time_array = np.linspace(time_start, time_finish, N_points)

    C_all_num = odeint(f, C_list_0, time_array)


    C_S_num = C_all_num[:,0]
    C_Z_num = C_all_num[:,1]
    C_R_num = C_all_num[:,2]

    plt.plot(time_array, C_S_num, 'r--',label='living')
    plt.plot(time_array, C_Z_num, 'b--',label='zombie')
    plt.plot(time_array, C_R_num, 'g--',label='dead')
    plt.xlabel('days after the outbreak')
    plt.ylabel('population')
    plt.legend()
    plt.show()

mymax = 0.1
mymin = 0.00000
interact(russians_are_not_scared_of_zombies, d=(mymin, mymax, (mymax-mymin)/10))


# # Problem 2:
# 
# Calculate the functions $T(t), C_{A}(t)$ in a mixer:
# 
# Stirred Tank reactor:
# 
# $T_{in}, C_{A0} \to T_f, C_A$
# mixer
# $q = \dot{V} = 100 m^3/hr$
# $V = 100 m^3$
# 
# 
# energy balance:
# 
# $V \dfrac{dT}{dt} = q (T_{f} - T(t))$
# 
# mass balance:
# 
# $V \dfrac{dC_A}{dt} = q (C_{f} - C_{A}(t))$
# 
# You need to give the program parameters 
# $C_{A0}=0, C_{Af} = 1$
# $T_0 = 350$, $T_f = 300$
# 
# you can call the python function:
# 
# ```python
# Tf, Cf = ...
# def mixer(X_list, t):
#     ...
#     return [dTdt, dCa/dt]
# ```

# In[ ]: