# This is a bifurcation diagram program for the replicator equations

from pylab import *
from MADSEngine import *
from math import sqrt

# Define a rudimentary Poincare secion
# This is the cut about the center of the simplex
# by a horizontal line
Mid = CoorTrans(array([1./3.,1./3.,1./3.]))[1]

# Define parameter range
parameter_low = .0
parameter_high = .05

# Define parameters and set initial conditions
dt = .01
nTransients = 90000 
nIterates = 100000

ex = 0.5
ey = -.5
alphaX = 0.02 
alphaY = 0.0

betaX = .5
betaY = .5

u = -log(array([0.2,0.6,0.2]))
v = -log(array([0.1,0.5,0.4]))

# This sets how dense the bifurcation diagram will be
nSteps = input("How dense should the diagram be?\nWarning: this program is very slow. Unless you have a super computer, pick 2 or 3:  ")

# Sweep the control parameter over the desired range
parameterInc = (parameter_high-parameter_low)/float(nSteps)

for r in arange(parameter_low,parameter_high,parameterInc):
	  XAgent, YAgent = MADS_Simulation(ex,ey,r,alphaY,betaX,betaY,u,v,nIterates,dt)
	  
	  parameter_sweep = []
	  X = []

	  for i in range(nTransients, nIterates-1):

	     if (CoorTrans(XAgent[i])[1] >= Mid) and (CoorTrans(XAgent[i+1])[1] <= Mid):
		  X.append(CoorTrans(XAgent[i])[0])
		  parameter_sweep.append(r)
	     elif (CoorTrans(XAgent[i])[1] <= Mid) and (CoorTrans(XAgent[i+1])[1]  >= Mid):
		  X.append(CoorTrans(XAgent[i])[0])
		  parameter_sweep.append(r)
	  	  
          plot(parameter_sweep,X,'ro')
	  
	  
show()