#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
%file synchronized_oscillators.m
%Hasty synthetic oscillator model
%from Hasty et al. (2002) Phys. Rev. Lett. 88 148101
%Figure 7.31

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



#assign parameter values

alpha=11
sigma=2
gammax=0.2
gammay=0.012
ay=0.2
D=0.015



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


#assign initial condition.  
#state vector is  [x1 y1 x2 y2] 
Sinit=[0.3963, 2.3346, 0.5578, 1.9317]
 
#declare right-hand-side for original model
def dSdt_original(S,t):
    dS=np.zeros(4) #generate a list to store derivatives
    
    x=S[0]
    y=S[1]
    xx=S[2]
    yy=S[3]
    
    
    dS[0]=(1+x**2+alpha*sigma*x**4)/((1+x**2+sigma*x**4)*(1+y**4)) - gammax*x + D*(xx-x)
    dS[1]=(ay)*((1+x**2+alpha*sigma*x**4)/((1+x**2+sigma*x**4)*(1+y**4))) - gammay*y
    dS[2]=(1+xx**2+alpha*sigma*xx**4)/((1+xx**2+sigma*xx**4)*(1+yy**4)) - gammax*xx + D*(x-xx)
    dS[3]=(ay)*((1+xx**2+alpha*sigma*xx**4)/((1+xx**2+sigma*xx**4)*(1+yy**4))) - gammay*yy

    return dS


S=odeint(dSdt_original, Sinit, times) #run simulation
   
 #produce figure 7.31
 
 
 
plt.figure()
plt.plot(times, S[:,2], label="X_1", linewidth=2)
plt.plot(times, S[:,3], label="Y_1", linewidth=2)
plt.plot(times, S[:,0], label="X_2", linewidth=2)
plt.plot(times, S[:,1], label="Y_2", linewidth=2)
plt.xlabel('Time (arbitrary units)')
plt.ylabel('Concentration (arbitrary units)')
plt.legend()