Suris's integrable map

\[\begin{split}p_{n+1} & = p_{n} + \frac{\epsilon}{2\pi}\frac{ \sin2\pi q}{1+\epsilon \cos2\pi q} \\ q_{n+1} &= q_{n} + p_{n+1}\end{split}\]

[1] Y. B. Suris, Integrable Mappings of the Standard Type, Funct Anal Its Appl 23, 74 (1989).

[2]J. D. Meiss, Symplectic Maps, Variational Principles, and Transport, Reviews of Modern Physics 64, 795 (1992).

[3] H. E. Lomelı́ and J. D. Meiss, Heteroclinic Orbits and Flux in a Perturbed Integrable Suris Map, Physics Letters A 269, 309 (2000).

import numpy as np
import matplotlib.pyplot as plt

twopi = 2*np.pi
omega=twopi

def Map(q,p,eps):
    pp = p + V(q,eps)
    qq = q + pp
    return qq, pp

def V(q,eps):
    return -np.arctan(
            eps*np.sin(twopi*q) / (1+eps*np.cos(twopi*q) )
        )*2/twopi


def Phi0(x,y):
    return (1-np.cos(omega*(x-y)))#/omega**2
def Phi1(x,y):
    return (np.cos(omega*x) + np.cos(omega*y))#/(2*omega)

def Phi(x,y,eps):
    return np.cos(omega*y) + eps*(np.cos(omega*x) + np.cos(omega*(x-y)))



def GetTraj(inits, k,sample=100,tmax=500):

    res = [[],[]]
    for q,p in zip(inits[0],inits[1]):
        traj = [np.array([])]*2
        for i in range(tmax):
            traj[0]  = np.append(traj[0], q)
            traj[1]  = np.append(traj[1], p)
            q,p = Map(q,p,eps)
            q = q - np.floor(q)
#            p = p - np.floor(p)
        res[0].append(traj[0])
        res[1].append(traj[1])
    return res


import sys
eps = 0.5

#tau = 1
sample=20
tmax = 200
x = np.linspace(0,1,1000,endpoint=False)
y = np.linspace(-1,1,1000)
X,Y = np.meshgrid(x,y)

fig1 = plt.figure(figsize=(6,6))
ax1 = fig1.add_subplot(1,1,1)

init = [np.array([1/2]),np.array([0])]
traj = GetTraj(init, eps,sample=sample,tmax=tmax)
energy = Phi(init[0][0], init[1][0],eps)

ax1.contour(X,Y,Phi(X,Y,eps),[energy],colors=["k"],linewidths=1.5)

q0 = 0.5+np.arccos(eps/2)/twopi
q1 = 1
p0 = np.arccos(eps**2/2 - 1)/twopi
p1 = 0

q = np.linspace(q0, q1, sample)
p = np.linspace(p0, p1, sample)

init = [q,p]
traj = GetTraj(init, eps,sample=sample,tmax=tmax)

energys = []
colors = []
for i in range(sample):
    card = i/sample
    lx, = ax1.plot(traj[0][i],traj[1][i],'o',markersize=4,mec="none", color=plt.cm.Oranges_r(card) )
    energy = Phi(init[0][i], init[1][i],eps)
    energys.append(energy)
    colors.append(lx.get_color())
clbar = ax1.contour(X,Y,Phi(X,Y,eps),energys, colors=colors)

ax1.set_ylim(-0.7,0.7)
ax1.set_xlim(0,1)
ax1.set_title(r"$k=%.1f$" % (eps),fontsize=16,position=(0.18,1))
ax1.set_xlabel(r"$q$",fontsize=25)
ax1.set_ylabel(r"$p$",fontsize=25)

fig1.tight_layout()
#fig1.savefig("SurisMap_k%.3f.png" % (k),transparent=True)
plt.show()