HenonHeiles.py

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Jan 30 10:45:00 2019

@author: David Nolte
Introduction to Modern Dynamics, 2nd edition (Oxford University Press, 2019)

Henon-Heiles model of galactic dynamics
"""

import numpy as np
from scipy import integrate
from matplotlib import pyplot as plt

plt.close('all')
print(' ')
print('HenonHeiles.py')

E = 1
epsE = 0.35
pmax = np.sqrt(2*E)

def flow_deriv(x_y_z_w, tspan):
    x, y, z, w = x_y_z_w
    a = z
    b = w
    c = -x - epsE*(2*x*y)
    d = -y - epsE*(x**2 - y**2)
    return [a,b,c,d]
    
px1 = 2*(np.random.random(1)-0.499)*pmax;
py1 = np.sign(np.random.random(1)-0.5)*np.sqrt(2*(E-px1**2/2));
xp1 = 0;
yp1 = 0;

x_y_z_w = [xp1, yp1, px1, py1]   

tspan = np.linspace(1,1000,10000)
x_t = integrate.odeint(flow_deriv, x_y_z_w, tspan)

y1 = x_t[:,0]
y2 = x_t[:,1]
y3 = x_t[:,2]
y4 = x_t[:,3]

plt.figure(1)
lines = plt.plot(y1,y2)
plt.setp(lines,linewidth=0.5)


repnum = 256                # 32
np.random.seed(1)
for reploop in range(2,repnum):
    
    px1 = 2*(np.random.random(1)-0.499)*pmax;
    py1 = np.sign(np.random.random(1)-0.5)*np.sqrt(2*(E-px1**2/2));
    xp1 = 0;
    yp1 = 0;

    x_y_z_w = [xp1, yp1, px1, py1]   

    tspan = np.linspace(1,1000,10000)
    x_t = integrate.odeint(flow_deriv, x_y_z_w, tspan)
    siztmp = np.shape(x_t)
    siz = siztmp[0]
    
    y1 = x_t[:,0]
    y2 = x_t[:,1]
    y3 = x_t[:,2]
    y4 = x_t[:,3]

    py = np.zeros(shape=(2*repnum,))
    yvar = np.zeros(shape=(2*repnum,))
    cnt = -1
    last = y1[1]
    for loop in range(2,siz):
        if (last < 0)and(y1[loop] > 0):
            cnt = cnt+1
            del1 = -y1[loop-1]/(y1[loop] - y1[loop-1])
            py[cnt] = y4[loop-1] + del1*(y4[loop]-y4[loop-1])
            yvar[cnt] = y2[loop-1] + del1*(y2[loop]-y2[loop-1])
            last = y1[loop]
        else:
            last = y1[loop]
 
    
    plt.figure(2)
    lines = plt.plot(yvar,py,'o',ms=1)
    plt.show()