1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:
70:
71:
72:
73:
74:
75:
76:
77:
78:
79:
80:
81:
82:
83:
84:
85:
86:
87:
88:
89:
90:

import matplotlib.pyplot as plt
import sympy as sy
import numpy as np
from scipy.integrate import odeint


"Poincare Surface of Section for z=0"


def f(Vector_z,t, epsilon):                                                 
    mu, nu, impuls_mu, impuls_nu = Vector_z
    return [impuls_mu, impuls_nu, 2*epsilon*mu - 3*mu*(nu**4) - 6*(mu**3)*nu**2, 
            2*epsilon*nu - 3*nu*(mu**4) - 6*(nu**3)*(mu**2)]

def get_trajectory(startVector, times, epsilon):
    return odeint(f, startVector, times, args=(epsilon,))

def zweiterImpuls(mu, nu, impuls):
    return (-np.sqrt(4 - impuls**2 + 2*epsilon*(mu**2 +nu**2) 
            - 3*(mu**2)*(nu**2)*(mu**2 +nu**2)))


if __name__ == '__main__':


    epsilon = -5.44667/2                  
    nu_PunktePoincare = []
    p_nu_PunktePoincare = []
    times = np.linspace(0,5,2000)
    
    #Anfangswerte festlegen
    mu = 0.43903
    nu = 0.43903
    impuls_mu = 0.9635 
    impuls_nu = zweiterImpuls(mu, nu, impuls_mu)

    startVector = [mu, nu, impuls_mu, impuls_nu]
    
    AllVectors =    [startVector, 
                    [0, 2/np.sqrt(-2*epsilon), 0, zweiterImpuls(0, 2/np.sqrt(-2*epsilon), 0)], 
                    [0, 1.7/np.sqrt(-2*epsilon), 0.4, zweiterImpuls(0, 1.7/np.sqrt(-2*epsilon), 0.4)],
                    [0, 1.7/np.sqrt(-2*epsilon), 0.8, zweiterImpuls(0, 1.7/np.sqrt(-2*epsilon), 0.8)],
                    [0, 1.4/np.sqrt(-2*epsilon), 0, zweiterImpuls(0, 1.4/np.sqrt(-2*epsilon), 0)],
                    [0, 1.4/np.sqrt(-2*epsilon), 0.4, zweiterImpuls(0, 1.4/np.sqrt(-2*epsilon), 0.4)],
                    [0, 1.4/np.sqrt(-2*epsilon), 0.8, zweiterImpuls(0, 1.4/np.sqrt(-2*epsilon), 0.8)], 
                    [0, 1.4/np.sqrt(-2*epsilon), 1.2, zweiterImpuls(0, 1.4/np.sqrt(-2*epsilon), 1.2)],
                    [0, 1.1/np.sqrt(-2*epsilon), 0, zweiterImpuls(0, 1.1/np.sqrt(-2*epsilon), 0)],
                    [0, 1.1/np.sqrt(-2*epsilon), 0.4, zweiterImpuls(0, 1.1/np.sqrt(-2*epsilon), 0.4)],
                    [0, 1.1/np.sqrt(-2*epsilon), 0.8, zweiterImpuls(0, 1.1/np.sqrt(-2*epsilon), 0.8)],  
                    [0, 1.1/np.sqrt(-2*epsilon), 1.2, zweiterImpuls(0, 1.1/np.sqrt(-2*epsilon), 1.2)], 
                    [0, 1.1/np.sqrt(-2*epsilon), 1.6, zweiterImpuls(0, 1.1/np.sqrt(-2*epsilon), 1.6)],
                    [0, 0.7/np.sqrt(-2*epsilon), 0, zweiterImpuls(0, 0.7/np.sqrt(-2*epsilon), 0)],
                    [0, 0.7/np.sqrt(-2*epsilon), 0.4, zweiterImpuls(0, 0.7/np.sqrt(-2*epsilon), 0.4)],
                    [0, 0.7/np.sqrt(-2*epsilon), 0.8, zweiterImpuls(0, 0.7/np.sqrt(-2*epsilon), 0.8)],
                    [0, 0.7/np.sqrt(-2*epsilon), 1.2, zweiterImpuls(0, 0.7/np.sqrt(-2*epsilon), 1.2)], 
                    [0, 0.7/np.sqrt(-2*epsilon), 1.8, zweiterImpuls(0, 0.7/np.sqrt(-2*epsilon), 1.8)],
                    [0, 0.4/np.sqrt(-2*epsilon), 0, zweiterImpuls(0, 0.4/np.sqrt(-2*epsilon), 0)],
                    [0, 0.4/np.sqrt(-2*epsilon), 0.4, zweiterImpuls(0, 0.4/np.sqrt(-2*epsilon), 0.4)],
                    [0, 0.4/np.sqrt(-2*epsilon), 0.8, zweiterImpuls(0, 0.4/np.sqrt(-2*epsilon), 0.8)],
                    [0, 0.4/np.sqrt(-2*epsilon), 1.2, zweiterImpuls(0, 0.4/np.sqrt(-2*epsilon), 1.2)], 
                    [0, 0.4/np.sqrt(-2*epsilon), 1.8, zweiterImpuls(0, 0.4/np.sqrt(-2*epsilon), 1.8)],
                    [0, 0.1/np.sqrt(-2*epsilon), 0, zweiterImpuls(0, 0.1/np.sqrt(-2*epsilon), 0)],
                    [0, 0.1/np.sqrt(-2*epsilon), 0.4, zweiterImpuls(0, 0.1/np.sqrt(-2*epsilon), 0.4)],
                    [0, 0.1/np.sqrt(-2*epsilon), 0.8, zweiterImpuls(0, 0.1/np.sqrt(-2*epsilon), 0.8)],
                    [0, 0.1/np.sqrt(-2*epsilon), 1.2, zweiterImpuls(0, 0.1/np.sqrt(-2*epsilon), 1.2)], 
                    [0, 0.1/np.sqrt(-2*epsilon), 1.8, zweiterImpuls(0, 0.1/np.sqrt(-2*epsilon), 1.8)]]
    
    
    for element in AllVectors:
        liste = get_trajectory(element, times, epsilon)
    

    # mu wird Null gesetzt, also kleiner als 0.003
    for element in liste:
        if(abs(element[0]) < 0.003):
            nu_PunktePoincare.append(element[1])
            p_nu_PunktePoincare.append(element[3])
            
    #nu noch skalieren, fürs spätere Plotten
    scaledPoincare_nuPoints = []
    for element in nu_PunktePoincare:
        temporary = element*np.sqrt(-2*epsilon)
        scaledPoincare_nuPoints.append(temporary)

    plt.plot(scaledPoincare_nuPoints, p_nu_PunktePoincare, 'b.')
    plt.xlabel(r'$\sqrt{-2 \epsilon} \nu$', fontsize=20)
    plt.ylabel(r'$p_{\nu}$', fontsize=20)
    plt.show()