#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Dec 12 17:37:29 2022
%problem5.6.2.py
%Problem 5.6.2


@author: bingalls
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint



#assign parameter values


v=10
k=1
k1=1
k2=1
k3=1
q=1



#set time_grid for simulation
times=np.linspace(0, 40, 1000) #generate time-grid list


#set initial conditions State vector is [s1, s2, s3]
Sinit=[0,0, 0];
#declare right-hand-side for original model
def dSdt_original(S,t):
    dS=[0,0,0] #generate a list to store derivatives
    s1=S[0]
    s2=S[1]
    s3=S[2]
    dS[0]=v/(1+np.power(s3/k,q) )- k1*s1
    dS[1]=k1*s1 -	k2*s2
    dS[2]=k2*s2 -	k3*s3
    return dS



S=odeint(dSdt_original, Sinit, times) #run simulation

#plot simulation
plt.figure() #generate figure
plt.plot(times, S[:,0], label="s1", linewidth=2)
plt.plot(times, S[:,1], label="s2", linewidth=2)
plt.plot(times, S[:,2], label="s3", linewidth=2)
plt.xlabel("Time")
plt.ylabel("concentration")
plt.legend()