#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Dec 12 17:37:29 2022
%ile rapid_equilibrium_approximation1.m
%igure 2.13
%imulation and rapid equilibrium approximation

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



#assign parameter values
k0=5
k1=20
km1=12
k2=2
#set time_grid for simulation
t_min=0; t_max=3; dt=0.01
times=np.arange(t_min, t_max+dt, dt) #generate time-grid list

###Original model#####
#set initial conditions for original model: S(1)=a, S(2)=b
S0=[8,4];
#declare right-hand-side for original model
def dSdt_original(S,t):
    dS=[0,0] #generate a list to store derivatives
    dS[0]=k0 - k1*S[0] + km1*S[1]
    dS[1]=k1*S[0] - km1*S[1] - k2*S[1]
    return dS


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




###Approximate model#####
#set initial conditions for approximate model
Sapprox0=[12];



#declare right-hand-side for approximate model
def dSdt_approx(S,t):
    dS=[k0 - (k2*k1)/(km1+k1)*S[0]]
    return dS

S_approx=odeint(dSdt_approx, Sapprox0, times) #run simulation

#plot simulation
plt.figure() #generate figure
plt.plot(times, km1/(k1+km1)*S_approx[:,0], label="a_approx", linewidth=2)
plt.plot(times, k1/(k1+km1)*S_approx[:,0], label="b_approx", linewidth=2)
plt.plot(times, S[:,0], label="a", linewidth=2)
plt.plot(times, S[:,1], label="b", linewidth=2)
plt.xlabel("Time")
plt.ylabel("concentration")
plt.legend()
plt.ylim(1,8)