# N=3 globally coupled phase oscillators, 2 harmonic coupling
# run using xppaut http://www.math.pitt.edu/~bard/xpp/xpp.html
#
# Peter Ashwin, Jan 2018
#
# simulates using phase differences
# phi1=theta1-theta3
# phi2=theta2-theta3
#
# See: https://www.frontiersin.org/articles/10.3389/fams.2016.00007/full
# for more details of this system

# initial conditions
init phi1=2.81 phi2=3.95

# parameters for coupling
par alpha=1.5 r=0.1 beta=0 q=-1
# omega (actually irrelevant for phase diffs)
par om=0.02
# amplitude of iid noise
par eta=0

# coupling function with 2 harmonics
g(xx)=q*sin(xx-alpha)+r*sin(2*xx-beta)

# for noise
wiener w[0..3]

# put into original phases 
theta1=phi1
theta2=phi2
theta3=0

# d/dt for original phases
dtheta1=om+(        g(theta1-theta2)+g(theta1-theta3))/3
dtheta2=om+(g(theta2-theta1)+        g(theta2-theta3))/3
dtheta3=om+(g(theta3-theta1)+g(theta3-theta2)        )/3

# phase diffs
phi1'=dtheta1-dtheta3+eta*(w[0]-w[2])
phi2'=dtheta2-dtheta3+eta*(w[1]-w[2])

# plot defaults
@ dt=0.2, total=2000, xplot=phi1, yplot=phi2
@ xlo=0, ylo=0, xhi=6.28, yhi=6.28, bound=1000, maxstor=100000
@ tor_per=6.28318531, fold=phi1, fold=phi2
done