Three body Problem

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Mar 12 13:56:29 2020

@author: gabriel
"""

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
 
#Inital Conditions
y0 = [-0.5000001,#x1
      0,#y1
      0.5 ,#x2
      0,#y2
      0,#x3
      0,#y3
     
      0,#vx1
      1,#vy1
      0,#vx2
      -1,#vy2
      0,#vx3
      0] #vy3

#Definition of the Function 
def Planets(t,y):
    f = np.zeros(12)
    f[0] = y[6]
    f[1] = y[7]
    f[2] = y[8]
    f[3] = y[9]
    f[4] = y[10]
    f[5] = y[11]
    
    f[6] = -(y[0]-y[2])/(((y[0]-y[2])**2+(y[1]-y[3])**2)**(3/2)) \
           -(y[0]-y[4])/(((y[0]-y[4])**2+(y[1]-y[5])**2)**(3/2))
           
    f[7] = -(y[1]-y[3])/(((y[0]-y[2])**2+(y[1]-y[3])**2)**(3/2)) \
           -(y[1]-y[5])/(((y[0]-y[4])**2+(y[1]-y[5])**2)**(3/2))
                     
    f[8] = -(y[2]-y[0])/(((y[2]-y[0])**2+(y[3]-y[1])**2)**(3/2)) \
           -(y[2]-y[4])/(((y[2]-y[4])**2+(y[3]-y[5])**2)**(3/2))
           
    f[9] = -(y[3]-y[1])/(((y[2]-y[0])**2+(y[3]-y[1])**2)**(3/2)) \
           -(y[3]-y[5])/(((y[2]-y[4])**2+(y[3]-y[5])**2)**(3/2))
                     
    f[10]= -(y[4]-y[0])/(((y[4]-y[0])**2+(y[5]-y[1])**2)**(3/2)) \
           -(y[4]-y[2])/(((y[4]-y[2])**2+(y[5]-y[3])**2)**(3/2))
          
    f[11]= -(y[5]-y[1])/(((y[4]-y[0])**2+(y[5]-y[1])**2)**(3/2)) \
           -(y[5]-y[3])/(((y[4]-y[2])**2+(y[5]-y[3])**2)**(3/2))    
           
    return(f)
    
N = 10000
T = 0.001 #N*T=10

t = np.linspace(0,N*T,N)
solution = solve_ivp(Planets,[0,800],y0,t_eval=t,rtol=1e-15)

#Evolution in Position with respect to Time
plt.plot(solution.y[0],solution.y[1],'-g') #Positions planet 1
plt.plot(solution.y[2],solution.y[3],'-r') #Positions planet 2
plt.plot(solution.y[4],solution.y[5],'-b') #Positions planet 3
plt.ylabel("Position(y)")
plt.xlabel("Position(x)")
plt.show()

plt.plot(t,solution.y[6])
plt.ylabel("Position (x)")
plt.xlabel("Time")
plt.show()

#Zeropadding
zeroN = 200000
totalN = zeroN + N
zeropadded_y = np.zeros(totalN)
zeropadded_y[:N] = (solution.y[6]) #Zero-pad vector position in x for planet 1

#Fourier Transform
yf = np.fft.rfft((zeropadded_y))
tf = np.linspace(0.0, 1.0//(2*T), totalN/2+1)

#Plot of the FFT
plt.plot(tf,abs(yf)**2) #Plotting the absolute value of the transform
plt.ylabel("|FFT|^2")
plt.xlabel("frequency")
#plt.ylim(0,1000000)
plt.xlim(0,10)
plt.show()




#THE KEY FOR MEANINFUL DATA IS TO MODIFY X AND Y LIMITS!!!
#plt.ylim(0,20000)



#Theoretical Book
#w0 = 6.44 
#T0 = (2*np.pi)/w0
#n = 2**13
#Tmax = 8*T0
#dt = Tmax/n
#tdomain = np.linspace(0,dt,Tmax-dt)