# Copyright (c) 2025 Jakub Szyndler
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>

import os
os.environ["KMP_DUPLICATE_LIB_OK"]="TRUE"
import numpy as np
import sys
import matplotlib.pyplot as plt
import matplotlib as mpl
from bulStoer import *
import random
import imageio.v2 as imageio
from os import walk
import glob
import subprocess
from mpl_toolkits.mplot3d import Axes3D
from input import *



def find_peryapsis(r_sep):
	i_peryapsis=[]
	for i in range(2,len(r_sep)):
		if r_sep[i-2]>r_sep[i-1] and r_sep[i]>r_sep[i-1]:
			i_peryapsis.append(i-1)
	return(i_peryapsis)
	

def makeGIF(path, name): #this program will be useful for animation
    filenames = next(walk(path), (None, None, []))[2]  # [] if no file
    filenames.sort()
    images = []
    for filename in filenames:
        images.append(imageio.imread(path + '/' + filename))
    imageio.mimsave(name+'.gif', images, fps=30)
@jit(target_backend='cuda',nopython=True)
def make_r(X):
	return(np.sqrt((X[0]-X[3])**2+(X[1]-X[4])**2+(X[2]-X[5])**2))	
@jit(target_backend='cuda',nopython=True)
def make_v(X):
	return(np.sqrt((X[6]-X[9])**2+(X[7]-X[10])**2+(X[8]-X[11])**2))
@jit(target_backend='cuda',nopython=True)
def F(t,X,m,eta):
	r = make_r(X) 
	v = make_v(X)
	F = np.zeros(12)
	
	n_12_x = (X[0]-X[3])/r
	n_12_y = (X[1]-X[4])/r
	n_12_z = (X[2]-X[5])/r
	
	n_21_x = -(X[0]-X[3])/r
	n_21_y = -(X[1]-X[4])/r
	n_21_z = -(X[2]-X[5])/r
	
	v_12_x = (X[6]-X[9])
	v_12_y = (X[7]-X[10])
	v_12_z = (X[8]-X[11])
	
	v_21_x = -(X[6]-X[9])
	v_21_y = -(X[7]-X[10])
	v_21_z = -(X[8]-X[11])
	
	v_1 = np.sqrt(X[6]**2 + X[7]**2 + X[8]**2)
	v_2 = np.sqrt(X[9]**2 + X[10]**2 + X[11]**2)
	v_12 = np.sqrt(v_12_x**2 + v_12_y**2 + v_12_z**2)
	v_21 = np.sqrt(v_21_x**2 + v_21_y**2 + v_21_z**2)
	
	n_12_v_1 = n_12_x*X[6] + n_12_y*X[7] + n_12_z*X[8]
	n_12_v_2 = n_12_x*X[9] + n_12_y*X[10] + n_12_z*X[11]
	n_21_v_1 = n_21_x*X[6] + n_21_y*X[7] + n_21_z*X[8]
	n_21_v_2 = n_21_x*X[9] + n_21_y*X[10] + n_21_z*X[11]
	n_12_v_12 = n_12_x*v_12_x + n_12_y*v_12_y + n_12_z*v_12_z
	n_21_v_12 = n_21_x*v_12_x + n_21_y*v_12_y + n_21_z*v_12_z
	n_12_v_21 = n_12_x*v_21_x + n_12_y*v_21_y + n_12_z*v_21_z
	n_21_v_21 = n_21_x*v_21_x + n_21_y*v_21_y + n_21_z*v_21_z
	v_1_v_2 = X[6]*X[9] + X[7]*X[10] + X[8]*X[11]  
	#Now lets define PN approximation based on Pati 2002 et.al and Blanchet 2024
	
	# 0-th order
	
	A_0_1_x = -m_2*n_12_x/r**2 
	A_0_2_x = -m_1*n_21_x/r**2 
	A_0_1_y = -m_2*n_12_y/r**2 
	A_0_2_y = -m_1*n_21_y/r**2
	A_0_1_z = -m_2*n_12_z/r**2 
	A_0_2_z = -m_1*n_21_z/r**2 
	
	# 1-th order

	
	A_1_1_x = ((5*m_1*m_2/r**3 + 4*m_2**2/r**3 +m_2/r**2*(1.5*n_12_v_2**2 - v_1**2 + 4*v_1_v_2 - 2*v_2**2 ) )*n_12_x + \
	(4*n_12_v_1 - 3*n_12_v_2)*m_2*v_12_x/r**2)
	
	A_1_2_x = ((5*m_1*m_2/r**3 + 4*m_1**2/r**3 + m_1/r**2*(1.5*n_21_v_1**2 - v_2**2 + 4*v_1_v_2 - 2*v_1**2 ) )*n_21_x + \
	(4*n_21_v_2 - 3*n_21_v_1)*m_1*v_21_x/r**2)
	
	A_1_1_y = ((5*m_1*m_2/r**3 + 4*m_2**2/r**3 + m_2/r**2*(1.5*n_12_v_2**2 - v_1**2 + 4*v_1_v_2 - 2*v_2**2 ) )*n_12_y + \
	(4*n_12_v_1 - 3*n_12_v_2)*m_2*v_12_y/r**2)
	
	A_1_2_y = ((5*m_1*m_2/r**3 + 4*m_1**2/r**3 + m_1/r**2*(1.5*n_21_v_1**2 - v_2**2 + 4*v_1_v_2 - 2*v_1**2 ) )*n_21_y + \
	(4*n_21_v_2 - 3*n_21_v_1)*m_1*v_21_y/r**2)
	
	A_1_1_z = ((5*m_1*m_2/r**3 + 4*m_2**2/r**3 + m_2/r**2*(1.5*n_12_v_2**2 - v_1**2 + 4*v_1_v_2 - 2*v_2**2 ) )*n_12_z + \
	(4*n_12_v_1 - 3*n_12_v_2)*m_2*v_12_z/r**2)
	
	A_1_2_z = ((5*m_1*m_2/r**3 + 4*m_1**2/r**3 + m_1/r**2*(1.5*n_21_v_1**2 - v_2**2 + 4*v_1_v_2 - 2*v_1**2 ) )*n_21_z + \
	(4*n_21_v_2 - 3*n_21_v_1)*m_1*v_21_z/r**2)
	
	# 2-th order

	
	A_2_1_x = ((-(57/4*m_1**2*m_2 + 69/2*m_1*m_2**2 + 9*m_2**3)/r**4 + \
	m_2*(-15/8*n_12_v_2**4 + 1.5*n_12_v_2**2*v_1**2 - 6*n_12_v_2**2*v_1_v_2 - 9*v_1_v_2**2 + 4.5*(n_12_v_2)*v_2**2 + 4*v_1_v_2*v_2**2 - 2*v_2**4 )/r**2 +\
	(m_1*m_2*(39/2*n_12_v_1**2 - 39*n_12_v_1*n_12_v_2 + 17/2*(n_12_v_2)**2 - 15/4*v_1**2 - 5/2*v_1_v_2 + 5/4*v_2**2 ) + \
	m_2**2*( 2*n_12_v_1**2 - 4*n_12_v_1*n_12_v_2 - 6*n_12_v_2**2 - 8*v_1_v_2 + 4*v_2**2))/r**3)*n_12_x + \
	((-2*m_2**2*(n_12_v_1 + n_12_v_2) + m_1*m_2*(-63*n_12_v_1 + 55*n_12_v_2)/4)/r**3)*v_12_x)
	
	A_2_2_x = ((-(57/4*m_2**2*m_1 + 69/2*m_2*m_1**2 + 9*m_1**3)/r**4 + \
	m_1*(-15/8*n_21_v_1**4 + 1.5*n_21_v_1**2*v_2**2 - 6*n_21_v_1**2*v_1_v_2 - 9*v_1_v_2**2 + 4.5*(n_21_v_1)*v_1**2 + 4*v_1_v_2*v_1**2 - 2*v_1**4 )/r**2 +\
	(m_2*m_1*(39/2*n_21_v_2**2 - 39*n_21_v_2*n_21_v_1 + 17/2*(n_21_v_1)**2 - 15/4*v_2**2 - 5/2*v_1_v_2 + 5/4*v_1**2 ) + \
	m_1**2*( 2*n_21_v_2**2 - 4*n_21_v_2*n_21_v_1 - 6*n_21_v_1**2 - 8*v_1_v_2 + 4*v_1**2))/r**3)*n_21_x + \
	((-2*m_1**2*(n_21_v_2 + n_21_v_1) + m_1*m_2*(-63*n_21_v_2 + 55*n_21_v_1)/4)/r**3)*v_21_x)
	
	A_2_1_y = ((-(57/4*m_1**2*m_2 + 69/2*m_1*m_2**2 + 9*m_2**3)/r**4 + \
	m_2*(-15/8*n_12_v_2**4 + 1.5*n_12_v_2**2*v_1**2 - 6*n_12_v_2**2*v_1_v_2 - 9*v_1_v_2**2 + 4.5*(n_12_v_2)*v_2**2 + 4*v_1_v_2*v_2**2 - 2*v_2**4 )/r**2 +\
	(m_1*m_2*(39/2*n_12_v_1**2 - 39*n_12_v_1*n_12_v_2 + 17/2*(n_12_v_2)**2 - 15/4*v_1**2 - 5/2*v_1_v_2 + 5/4*v_2**2 ) + \
	m_2**2*( 2*n_12_v_1**2 - 4*n_12_v_1*n_12_v_2 - 6*n_12_v_2**2 - 8*v_1_v_2 + 4*v_2**2))/r**3)*n_12_y + \
	((-2*m_2**2*(n_12_v_1 + n_12_v_2) + m_1*m_2*(-63*n_12_v_1 + 55*n_12_v_2)/4)/r**3)*v_12_y)
	
	A_2_2_y = ((-(57/4*m_2**2*m_1 + 69/2*m_2*m_1**2 + 9*m_1**3)/r**4 + \
	m_1*(-15/8*n_21_v_1**4 + 1.5*n_21_v_1**2*v_2**2 - 6*n_21_v_1**2*v_1_v_2 - 9*v_1_v_2**2 + 4.5*(n_21_v_1)*v_1**2 + 4*v_1_v_2*v_1**2 - 2*v_1**4 )/r**2 +\
	(m_2*m_1*(39/2*n_21_v_2**2 - 39*n_21_v_2*n_21_v_1 + 17/2*(n_21_v_1)**2 - 15/4*v_2**2 - 5/2*v_1_v_2 + 5/4*v_1**2 ) + \
	m_1**2*( 2*n_21_v_2**2 - 4*n_21_v_2*n_21_v_1 - 6*n_21_v_1**2 - 8*v_1_v_2 + 4*v_1**2))/r**3)*n_21_y + \
	((-2*m_1**2*(n_21_v_2 + n_21_v_1) + m_1*m_2*(-63*n_21_v_2 + 55*n_21_v_1)/4)/r**3)*v_21_y)
	
	A_2_1_z = ((-(57/4*m_1**2*m_2 + 69/2*m_1*m_2**2 + 9*m_2**3)/r**4 + \
	m_2*(-15/8*n_12_v_2**4 + 1.5*n_12_v_2**2*v_1**2 - 6*n_12_v_2**2*v_1_v_2 - 9*v_1_v_2**2 + 4.5*(n_12_v_2)*v_2**2 + 4*v_1_v_2*v_2**2 - 2*v_2**4 )/r**2 +\
	(m_1*m_2*(39/2*n_12_v_1**2 - 39*n_12_v_1*n_12_v_2 + 17/2*(n_12_v_2)**2 - 15/4*v_1**2 - 5/2*v_1_v_2 + 5/4*v_2**2 ) + \
	m_2**2*( 2*n_12_v_1**2 - 4*n_12_v_1*n_12_v_2 - 6*n_12_v_2**2 - 8*v_1_v_2 + 4*v_2**2))/r**3)*n_12_z + \
	((-2*m_2**2*(n_12_v_1 + n_12_v_2) + m_1*m_2*(-63*n_12_v_1 + 55*n_12_v_2)/4)/r**3)*v_12_z)
	
	A_2_2_z = ((-(57/4*m_2**2*m_1 + 69/2*m_2*m_1**2 + 9*m_1**3)/r**4 + \
	m_1*(-15/8*n_21_v_1**4 + 1.5*n_21_v_1**2*v_2**2 - 6*n_21_v_1**2*v_1_v_2 - 9*v_1_v_2**2 + 4.5*(n_21_v_1)*v_1**2 + 4*v_1_v_2*v_1**2 - 2*v_1**4 )/r**2 +\
	(m_2*m_1*(39/2*n_21_v_2**2 - 39*n_21_v_2*n_21_v_1 + 17/2*(n_21_v_1)**2 - 15/4*v_2**2 - 5/2*v_1_v_2 + 5/4*v_1**2 ) + \
	m_1**2*( 2*n_21_v_2**2 - 4*n_21_v_2*n_21_v_1 - 6*n_21_v_1**2 - 8*v_1_v_2 + 4*v_1**2))/r**3)*n_21_z + \
	((-2*m_1**2*(n_21_v_2 + n_21_v_1) + m_1*m_2*(-63*n_21_v_2 + 55*n_21_v_1)/4)/r**3)*v_21_z)
	
	# 2.5-th order

	
	A_3_1_x = ((m_1*m_2*n_12_v_12*(208*m_2/15 - 24*m_1/5)/r**4 + 12*m_1*m_2*n_12_v_12*v_12**2/5/r**3)*n_12_x + \
	(m_1*m_2*(8*m_1/5 - 32*m_2/5)/r**4-4*G**2*m_1*m_2*v_12**2/5/r**3)*v_12_x)
	
	A_3_2_x = ((m_2*m_1*n_21_v_21*(208*m_1/15 - 24*m_2/5)/r**4 + 12*m_2*m_1*n_21_v_21*v_21**2/5/r**3)*n_21_x + \
	(m_2*m_1*(8*m_2/5 - 32*m_1/5)/r**4-4*G**2*m_2*m_1*v_21**2/5/r**3)*v_21_x)
	
	A_3_1_y = ((m_1*m_2*n_12_v_12*(208*m_2/15 - 24*m_1/5)/r**4 + 12*m_1*m_2*n_12_v_12*v_12**2/5/r**3)*n_12_y + \
	(m_1*m_2*(8*m_1/5 - 32*m_2/5)/r**4-4*m_1*m_2*v_12**2/5/r**3)*v_12_y)
	
	A_3_2_y = ((m_2*m_1*n_21_v_21*(208*m_1/15 - 24*m_2/5)/r**4 + 12*m_2*m_1*n_21_v_21*v_21**2/5/r**3)*n_21_y + \
	(m_2*m_1*(8*m_2/5 - 32*m_1/5)/r**4-4*m_2*m_1*v_21**2/5/r**3)*v_21_y)
	
	A_3_1_z = ((m_1*m_2*n_12_v_12*(208*m_2/15 - 24*m_1/5)/r**4 + 12*m_1*m_2*n_12_v_12*v_12**2/5/r**3)*n_12_z + \
	(m_1*m_2*(8*m_1/5 - 32*m_2/5)/r**4-4*m_1*m_2*v_12**2/5/r**3)*v_12_z)
	
	A_3_2_z = -((m_2*m_1*n_21_v_21*(208*m_1/15 - 24*m_2/5)/r**4 + 12*m_2*m_1*n_21_v_21*v_21**2/5/r**3)*n_21_z + \
	(m_2*m_1*(8*m_2/5 - 32*m_1/5)/r**4-4*m_2*m_1*v_21**2/5/r**3)*v_21_z)
	
	# 3.5-th order

	
	A_4_1_x = ((m_1**3*m_2*(3992/105*n_12_v_1 - 4328/105*n_12_v_2))/r**5 + (m_1*m_2)**2*(-13576/105*n_12_v_1 + 2872/21*n_12_v_2)/r**5 \
	-3172/21*m_1*m_2**3*n_12_v_12/r**6 + m_1**2*m_2*(48*n_12_v_1**3 - 696/5*n_12_v_1**2*n_12_v_2 + 744/5*n_12_v_1*n_12_v_2**2 -288/5*n_12_v_2**3 \
	-4888/105*n_12_v_1*v_1**2 + 5056/105*n_12_v_2*v_1**2 + 2056/21*n_12_v_1*v_1_v_2 - 2224/21*n_12_v_2*v_1_v_2 - 1028/21*n_12_v_1*v_2**2 + 5812/105*n_12_v_2*v_2**2)/r**4 + \
	m_1*m_2**2*(-582/5*n_12_v_1**3 + 1746/5*n_12_v_1**2*n_12_v_2 - 1954/5*n_12_v_1*n_12_v_2**2 + 158*n_12_v_2**3 + 3568/105*n_12_v_12*v_1**2 \
	-2864/35*n_12_v_1*v_1_v_2 + 10048/105*n_12_v_2*v_1_v_2 + 1432/35*n_12_v_1*v_2**2 - 5752/105*n_12_v_2*v_2**2 )/r**4 + \
	m_1*m_2*(-56*(n_12_v_12)**5 + 60*n_12_v_1**3*v_12**2 - 180*n_12_v_1**2*n_12_v_2*v_12**2 + 174*n_12_v_2**2*n_12_v_1*v_12**2-54*n_12_v_2**3*v_12**2 \
	- 246/35*n_12_v_12*v_1**4 + 1068/35*n_12_v_1*v_1**2*v_1_v_2 - 984/35*n_12_v_2*v_1**2*v_1_v_2 - 1068/35*n_12_v_1*v_1_v_2**2 + \
	+ 180/7*n_12_v_2*v_1_v_2**2 - 534/35*n_12_v_1*(v_1*v_2)**2 + 90/7*n_12_v_2*(v_1*v_2)**2 + 984/35*n_12_v_1*v_2**2*v_1_v_2 \
	- 732/35*n_12_v_2*v_2**2*v_1_v_2 - 204/35*n_12_v_1*v_2**4 + 24/7*n_12_v_2*v_2**4)/r**3)*n_12_x +  \
	(-184/21*m_1**3*m_2/r**5 + 6224/105*m_1**2*m_2**2/r**6 + 6338/105*m_1*m_2**3/r**6 + m_1**2*m_2*(52/15*n_12_v_1**2 \
	- 56/15*n_12_v_1*n_12_v_2 - 44/15*n_12_v_2**2 - 132/35*v_1**2 + 152/35*v_1_v_2 - 48/35*v_2**2)/r**4 + \
	m_1*m_2**2*(454/15*n_12_v_1**2 - 372/5*n_12_v_1*n_12_v_2 + 854/15*n_12_v_2**2 - 152/21*v_1**2 + 2864/105*v_1_v_2 - 1768/105*v_2**2)/r**4 \
	+ m_1*m_2*(60*n_12_v_12**4 - 348/5*n_12_v_1**2*v_12**2 + 684/5*n_12_v_1*n_12_v_2*v_12**2 - 66*n_12_v_2**2*v_12**2 \
	+ 344/35*v_1**4 - 1336/35*v_1**2*v_1_v_2 + 1308/35*v_1_v_2**2 + 654/35*v_1**2*v_2**2 - 1252/35*v_2**2*v_1_v_2 + 292/35*v_2**4)/r**3)*v_12_x
	
	A_4_2_x = (((m_2**3*m_1*(3992/105*n_21_v_2 - 4328/105*n_21_v_1))/r**5 + (m_2*m_1)**2*(-13576/105*n_21_v_2 + 2872/21*n_21_v_1)/r**5 \
	-3172/21*m_2*m_1**3*n_21_v_21/r**6 + m_2**2*m_1*(48*n_21_v_2**3 - 696/5*n_21_v_2**2*n_21_v_1 + 744/5*n_21_v_2*n_21_v_1**2 -288/5*n_21_v_1**3 \
	-4888/105*n_21_v_2*v_2**2 + 5056/105*n_21_v_1*v_2**2 + 2056/21*n_21_v_2*v_1_v_2 - 2224/21*n_21_v_1*v_1_v_2 - 1028/21*n_21_v_2*v_1**2 + 5812/105*n_21_v_1*v_1**2)/r**4 + \
	m_2*m_1**2*(-582/5*n_21_v_2**3 + 1746/5*n_21_v_2**2*n_21_v_1 - 1954/5*n_21_v_2*n_21_v_1**2 + 158*n_21_v_1**3 + 3568/105*n_21_v_21*v_2**2 \
	-2864/35*n_21_v_2*v_1_v_2 + 10048/105*n_21_v_1*v_1_v_2 + 1432/35*n_21_v_2*v_1**2 - 5752/105*n_21_v_1*v_1**2 )/r**4 + \
	m_2*m_1*(-56*(n_21_v_21)**5 + 60*n_21_v_2**3*v_21**2 - 180*n_21_v_2**2*n_21_v_1*v_21**2 + 174*n_21_v_1**2*n_21_v_2*v_21**2-54*n_21_v_1**3*v_21**2 \
	- 246/35*n_21_v_21*v_2**4 + 1068/35*n_21_v_2*v_2**2*v_1_v_2 - 984/35*n_21_v_1*v_2**2*v_1_v_2 - 1068/35*n_21_v_2*v_1_v_2**2 + \
	+ 180/7*n_21_v_1*v_1_v_2**2 - 534/35*n_21_v_2*(v_2*v_1)**2 + 90/7*n_21_v_1*(v_2*v_1)**2 + 984/35*n_21_v_2*v_1**2*v_1_v_2 \
	- 732/35*n_21_v_1*v_1**2*v_1_v_2 - 204/35*n_12_v_2*v_1**4 + 24/7*n_21_v_1*v_1**4)/r**3)*n_21_x +  \
	(-184/21*m_2**3*m_1/r**5 + 6224/105*m_2**2*m_1**2/r**6 + 6338/105*m_2*m_1**3/r**6 + m_2**2*m_1*(52/15*n_21_v_2**2 \
	- 56/15*n_21_v_2*n_21_v_1 - 44/15*n_21_v_1**2 - 132/35*v_2**2 + 152/35*v_1_v_2 - 48/35*v_1**2)/r**4 + \
	m_2*m_1**2*(454/15*n_21_v_2**2 - 372/5*n_21_v_2*n_21_v_1 + 854/15*n_21_v_1**2 - 152/21*v_2**2 + 2864/105*v_1_v_2 - 1768/105*v_1**2)/r**4 \
	+ m_2*m_1*(60*n_21_v_21**4 - 348/5*n_21_v_2**2*v_21**2 + 684/5*n_21_v_2*n_21_v_1*v_21**2 - 66*n_21_v_1**2*v_21**2 \
	+ 344/35*v_2**4 - 1336/35*v_2**2*v_1_v_2 + 1308/35*v_1_v_2**2 + 654/35*v_2**2*v_1**2 - 1252/35*v_1**2*v_1_v_2 + 292/35*v_1**4)/r**3)*v_21_x)
	
	A_4_1_y = ((m_1**3*m_2*(3992/105*n_12_v_1 - 4328/105*n_12_v_2))/r**5 + (m_1*m_2)**2*(-13576/105*n_12_v_1 + 2872/21*n_12_v_2)/r**5 \
	-3172/21*m_1*m_2**3*n_12_v_12/r**6 + m_1**2*m_2*(48*n_12_v_1**3 - 696/5*n_12_v_1**2*n_12_v_2 + 744/5*n_12_v_1*n_12_v_2**2 -288/5*n_12_v_2**3 \
	-4888/105*n_12_v_1*v_1**2 + 5056/105*n_12_v_2*v_1**2 + 2056/21*n_12_v_1*v_1_v_2 - 2224/21*n_12_v_2*v_1_v_2 - 1028/21*n_12_v_1*v_2**2 + 5812/105*n_12_v_2*v_2**2)/r**4 + \
	m_1*m_2**2*(-582/5*n_12_v_1**3 + 1746/5*n_12_v_1**2*n_12_v_2 - 1954/5*n_12_v_1*n_12_v_2**2 + 158*n_12_v_2**3 + 3568/105*n_12_v_12*v_1**2 \
	-2864/35*n_12_v_1*v_1_v_2 + 10048/105*n_12_v_2*v_1_v_2 + 1432/35*n_12_v_1*v_2**2 - 5752/105*n_12_v_2*v_2**2 )/r**4 + \
	m_1*m_2*(-56*(n_12_v_12)**5 + 60*n_12_v_1**3*v_12**2 - 180*n_12_v_1**2*n_12_v_2*v_12**2 + 174*n_12_v_2**2*n_12_v_1*v_12**2-54*n_12_v_2**3*v_12**2 \
	- 246/35*n_12_v_12*v_1**4 + 1068/35*n_12_v_1*v_1**2*v_1_v_2 - 984/35*n_12_v_2*v_1**2*v_1_v_2 - 1068/35*n_12_v_1*v_1_v_2**2 + \
	+ 180/7*n_12_v_2*v_1_v_2**2 - 534/35*n_12_v_1*(v_1*v_2)**2 + 90/7*n_12_v_2*(v_1*v_2)**2 + 984/35*n_12_v_1*v_2**2*v_1_v_2 \
	- 732/35*n_12_v_2*v_2**2*v_1_v_2 - 204/35*n_12_v_1*v_2**4 + 24/7*n_12_v_2*v_2**4)/r**3)*n_12_y +  \
	(-184/21*m_1**3*m_2/r**5 + 6224/105*m_1**2*m_2**2/r**6 + 6338/105*m_1*m_2**3/r**6 + m_1**2*m_2*(52/15*n_12_v_1**2 \
	- 56/15*n_12_v_1*n_12_v_2 - 44/15*n_12_v_2**2 - 132/35*v_1**2 + 152/35*v_1_v_2 - 48/35*v_2**2)/r**4 + \
	m_1*m_2**2*(454/15*n_12_v_1**2 - 372/5*n_12_v_1*n_12_v_2 + 854/15*n_12_v_2**2 - 152/21*v_1**2 + 2864/105*v_1_v_2 - 1768/105*v_2**2)/r**4 \
	+ m_1*m_2*(60*n_12_v_12**4 - 348/5*n_12_v_1**2*v_12**2 + 684/5*n_12_v_1*n_12_v_2*v_12**2 - 66*n_12_v_2**2*v_12**2 \
	+ 344/35*v_1**4 - 1336/35*v_1**2*v_1_v_2 + 1308/35*v_1_v_2**2 + 654/35*v_1**2*v_2**2 - 1252/35*v_2**2*v_1_v_2 + 292/35*v_2**4)/r**3)*v_12_y
	
	A_4_2_y = (((m_2**3*m_1*(3992/105*n_21_v_2 - 4328/105*n_21_v_1))/r**5 + (m_2*m_1)**2*(-13576/105*n_21_v_2 + 2872/21*n_21_v_1)/r**5 \
	-3172/21*m_2*m_1**3*n_21_v_21/r**6 + m_2**2*m_1*(48*n_21_v_2**3 - 696/5*n_21_v_2**2*n_21_v_1 + 744/5*n_21_v_2*n_21_v_1**2 -288/5*n_21_v_1**3 \
	-4888/105*n_21_v_2*v_2**2 + 5056/105*n_21_v_1*v_2**2 + 2056/21*n_21_v_2*v_1_v_2 - 2224/21*n_21_v_1*v_1_v_2 - 1028/21*n_21_v_2*v_1**2 + 5812/105*n_21_v_1*v_1**2)/r**4 + \
	m_2*m_1**2*(-582/5*n_21_v_2**3 + 1746/5*n_21_v_2**2*n_21_v_1 - 1954/5*n_21_v_2*n_21_v_1**2 + 158*n_21_v_1**3 + 3568/105*n_21_v_21*v_2**2 \
	-2864/35*n_21_v_2*v_1_v_2 + 10048/105*n_21_v_1*v_1_v_2 + 1432/35*n_21_v_2*v_1**2 - 5752/105*n_21_v_1*v_1**2 )/r**4 + \
	m_2*m_1*(-56*(n_21_v_21)**5 + 60*n_21_v_2**3*v_21**2 - 180*n_21_v_2**2*n_21_v_1*v_21**2 + 174*n_21_v_1**2*n_21_v_2*v_21**2-54*n_21_v_1**3*v_21**2 \
	- 246/35*n_21_v_21*v_2**4 + 1068/35*n_21_v_2*v_2**2*v_1_v_2 - 984/35*n_21_v_1*v_2**2*v_1_v_2 - 1068/35*n_21_v_2*v_1_v_2**2 + \
	+ 180/7*n_21_v_1*v_1_v_2**2 - 534/35*n_21_v_2*(v_2*v_1)**2 + 90/7*n_21_v_1*(v_2*v_1)**2 + 984/35*n_21_v_2*v_1**2*v_1_v_2 \
	- 732/35*n_21_v_1*v_1**2*v_1_v_2 - 204/35*n_12_v_2*v_1**4 + 24/7*n_21_v_1*v_1**4)/r**3)*n_21_y +  \
	(-184/21*m_2**3*m_1/r**5 + 6224/105*m_2**2*m_1**2/r**6 + 6338/105*m_2*m_1**3/r**6 + m_2**2*m_1*(52/15*n_21_v_2**2 \
	- 56/15*n_21_v_2*n_21_v_1 - 44/15*n_21_v_1**2 - 132/35*v_2**2 + 152/35*v_1_v_2 - 48/35*v_1**2)/r**4 + \
	m_2*m_1**2*(454/15*n_21_v_2**2 - 372/5*n_21_v_2*n_21_v_1 + 854/15*n_21_v_1**2 - 152/21*v_2**2 + 2864/105*v_1_v_2 - 1768/105*v_1**2)/r**4 \
	+ m_2*m_1*(60*n_21_v_21**4 - 348/5*n_21_v_2**2*v_21**2 + 684/5*n_21_v_2*n_21_v_1*v_21**2 - 66*n_21_v_1**2*v_21**2 \
	+ 344/35*v_2**4 - 1336/35*v_2**2*v_1_v_2 + 1308/35*v_1_v_2**2 + 654/35*v_2**2*v_1**2 - 1252/35*v_1**2*v_1_v_2 + 292/35*v_1**4)/r**3)*v_21_y)

	A_4_1_z = ((m_1**3*m_2*(3992/105*n_12_v_1 - 4328/105*n_12_v_2))/r**5 + (m_1*m_2)**2*(-13576/105*n_12_v_1 + 2872/21*n_12_v_2)/r**5 \
	-3172/21*m_1*m_2**3*n_12_v_12/r**6 + m_1**2*m_2*(48*n_12_v_1**3 - 696/5*n_12_v_1**2*n_12_v_2 + 744/5*n_12_v_1*n_12_v_2**2 -288/5*n_12_v_2**3 \
	-4888/105*n_12_v_1*v_1**2 + 5056/105*n_12_v_2*v_1**2 + 2056/21*n_12_v_1*v_1_v_2 - 2224/21*n_12_v_2*v_1_v_2 - 1028/21*n_12_v_1*v_2**2 + 5812/105*n_12_v_2*v_2**2)/r**4 + \
	m_1*m_2**2*(-582/5*n_12_v_1**3 + 1746/5*n_12_v_1**2*n_12_v_2 - 1954/5*n_12_v_1*n_12_v_2**2 + 158*n_12_v_2**3 + 3568/105*n_12_v_12*v_1**2 \
	-2864/35*n_12_v_1*v_1_v_2 + 10048/105*n_12_v_2*v_1_v_2 + 1432/35*n_12_v_1*v_2**2 - 5752/105*n_12_v_2*v_2**2 )/r**4 + \
	m_1*m_2*(-56*(n_12_v_12)**5 + 60*n_12_v_1**3*v_12**2 - 180*n_12_v_1**2*n_12_v_2*v_12**2 + 174*n_12_v_2**2*n_12_v_1*v_12**2-54*n_12_v_2**3*v_12**2 \
	- 246/35*n_12_v_12*v_1**4 + 1068/35*n_12_v_1*v_1**2*v_1_v_2 - 984/35*n_12_v_2*v_1**2*v_1_v_2 - 1068/35*n_12_v_1*v_1_v_2**2 + \
	+ 180/7*n_12_v_2*v_1_v_2**2 - 534/35*n_12_v_1*(v_1*v_2)**2 + 90/7*n_12_v_2*(v_1*v_2)**2 + 984/35*n_12_v_1*v_2**2*v_1_v_2 \
	- 732/35*n_12_v_2*v_2**2*v_1_v_2 - 204/35*n_12_v_1*v_2**4 + 24/7*n_12_v_2*v_2**4)/r**3)*n_12_z +  \
	(-184/21*m_1**3*m_2/r**5 + 6224/105*m_1**2*m_2**2/r**6 + 6338/105*m_1*m_2**3/r**6 + m_1**2*m_2*(52/15*n_12_v_1**2 \
	- 56/15*n_12_v_1*n_12_v_2 - 44/15*n_12_v_2**2 - 132/35*v_1**2 + 152/35*v_1_v_2 - 48/35*v_2**2)/r**4 + \
	m_1*m_2**2*(454/15*n_12_v_1**2 - 372/5*n_12_v_1*n_12_v_2 + 854/15*n_12_v_2**2 - 152/21*v_1**2 + 2864/105*v_1_v_2 - 1768/105*v_2**2)/r**4 \
	+ m_1*m_2*(60*n_12_v_12**4 - 348/5*n_12_v_1**2*v_12**2 + 684/5*n_12_v_1*n_12_v_2*v_12**2 - 66*n_12_v_2**2*v_12**2 \
	+ 344/35*v_1**4 - 1336/35*v_1**2*v_1_v_2 + 1308/35*v_1_v_2**2 + 654/35*v_1**2*v_2**2 - 1252/35*v_2**2*v_1_v_2 + 292/35*v_2**4)/r**3)*v_12_z
	
	A_4_2_z = (((m_2**3*m_1*(3992/105*n_21_v_2 - 4328/105*n_21_v_1))/r**5 + (m_2*m_1)**2*(-13576/105*n_21_v_2 + 2872/21*n_21_v_1)/r**5 \
	-3172/21*m_2*m_1**3*n_21_v_21/r**6 + m_2**2*m_1*(48*n_21_v_2**3 - 696/5*n_21_v_2**2*n_21_v_1 + 744/5*n_21_v_2*n_21_v_1**2 -288/5*n_21_v_1**3 \
	-4888/105*n_21_v_2*v_2**2 + 5056/105*n_21_v_1*v_2**2 + 2056/21*n_21_v_2*v_1_v_2 - 2224/21*n_21_v_1*v_1_v_2 - 1028/21*n_21_v_2*v_1**2 + 5812/105*n_21_v_1*v_1**2)/r**4 + \
	m_2*m_1**2*(-582/5*n_21_v_2**3 + 1746/5*n_21_v_2**2*n_21_v_1 - 1954/5*n_21_v_2*n_21_v_1**2 + 158*n_21_v_1**3 + 3568/105*n_21_v_21*v_2**2 \
	-2864/35*n_21_v_2*v_1_v_2 + 10048/105*n_21_v_1*v_1_v_2 + 1432/35*n_21_v_2*v_1**2 - 5752/105*n_21_v_1*v_1**2 )/r**4 + \
	m_2*m_1*(-56*(n_21_v_21)**5 + 60*n_21_v_2**3*v_21**2 - 180*n_21_v_2**2*n_21_v_1*v_21**2 + 174*n_21_v_1**2*n_21_v_2*v_21**2-54*n_21_v_1**3*v_21**2 \
	- 246/35*n_21_v_21*v_2**4 + 1068/35*n_21_v_2*v_2**2*v_1_v_2 - 984/35*n_21_v_1*v_2**2*v_1_v_2 - 1068/35*n_21_v_2*v_1_v_2**2 + \
	+ 180/7*n_21_v_1*v_1_v_2**2 - 534/35*n_21_v_2*(v_2*v_1)**2 + 90/7*n_21_v_1*(v_2*v_1)**2 + 984/35*n_21_v_2*v_1**2*v_1_v_2 \
	- 732/35*n_21_v_1*v_1**2*v_1_v_2 - 204/35*n_12_v_2*v_1**4 + 24/7*n_21_v_1*v_1**4)/r**3)*n_21_z +  \
	(-184/21*m_2**3*m_1/r**5 + 6224/105*m_2**2*m_1**2/r**6 + 6338/105*m_2*m_1**3/r**6 + m_2**2*m_1*(52/15*n_21_v_2**2 \
	- 56/15*n_21_v_2*n_21_v_1 - 44/15*n_21_v_1**2 - 132/35*v_2**2 + 152/35*v_1_v_2 - 48/35*v_1**2)/r**4 + \
	m_2*m_1**2*(454/15*n_21_v_2**2 - 372/5*n_21_v_2*n_21_v_1 + 854/15*n_21_v_1**2 - 152/21*v_2**2 + 2864/105*v_1_v_2 - 1768/105*v_1**2)/r**4 \
	+ m_2*m_1*(60*n_21_v_21**4 - 348/5*n_21_v_2**2*v_21**2 + 684/5*n_21_v_2*n_21_v_1*v_21**2 - 66*n_21_v_1**2*v_21**2 \
	+ 344/35*v_2**4 - 1336/35*v_2**2*v_1_v_2 + 1308/35*v_1_v_2**2 + 654/35*v_2**2*v_1**2 - 1252/35*v_1**2*v_1_v_2 + 292/35*v_1**4)/r**3)*v_21_z)

	#Lets define a differential equations system
	
	# dx_n/dt = v_n
	F[0] = X[6] 
	F[1] = X[7]
	F[2] = X[8]
	F[3] =  X[9]
	F[4] = X[10]
	F[5] = X[11]
	
	#dx^2_n/dt = a
	
	F[6] = A_0_1_x + A_1_1_x + A_2_1_x + A_3_1_x + A_4_1_x
	F[7] = A_0_1_y + A_1_1_y + A_2_1_y + A_3_1_y + A_4_1_y
	F[8] = A_0_1_z + A_1_1_z + A_2_1_z + A_3_1_z + A_4_1_z
	F[9] = A_0_2_x + A_1_2_x + A_2_2_x + A_3_2_x + A_4_2_x
	
	F[10] = A_0_2_y + A_1_2_y + A_2_2_y + A_3_2_y + A_4_2_y
	
	F[11] = A_0_2_z + A_1_2_z + A_2_2_z + A_3_2_z + A_4_2_z
	return(F)

@jit(target_backend='cuda',nopython=True)
def make_delta_T(T_peryapsis):
	delta_T = np.zeros(len(T_peryapsis))
	for i in range(1,len(T_peryapsis)):
		delta_T[i] = T_peryapsis[i] - T_peryapsis[i-1]
	return(delta_T)
if __name__=="__main__":
	# WE operate in center of mass frame
	x_0_1 = -r_0*m_1/M
	x_0_2  = r_0*m_2/M
	y_0_1 = 0
	y_0_2  = 0
	z_0_1 = 0
	z_0_2  = 0
	#r_0 = np.sqrt((x_0_1 - x_0_2)**2 + (y_0_1 - y_0_2)**2 + (z_0_1 - z_0_2)**2)
	# initial velocity, an fraction of c , we assume orbits are circular at the start
	v_0_1_x = 0
	v_0_2_x =  0
	v_0_1_y = m_2*np.sqrt(1/M/r_0)
	v_0_2_y = -m_1*np.sqrt(1/M/r_0)
	v_0_1_z = 0
	v_0_2_z = 0
	v_0_1 = np.sqrt(v_0_1_x**2+v_0_1_y**2+v_0_1_z**2)
	v_0_2 = np.sqrt(v_0_2_x**2+v_0_2_y**2+v_0_2_z**2)
	X = np.array([x_0_1,y_0_1,z_0_1,x_0_2,y_0_2,z_0_2,v_0_1_x,v_0_1_y,v_0_1_z,v_0_2_x,v_0_2_y,v_0_2_z]) #columns 1-3 position of first body (xyz), columns 4-6 position of second body, 6-9 velocities of first body , 10 -12 velocities of second body
	t_start = 0.0
	t_stop = t0 #min([t0,5/256*r_0**4/(M**2*eta)*c/365.25]) # We set the maximum time as a time of coalscence in years 
	print(f"The system contains two bodies of masses m_1={m_1*c**2/G},m_2={m_2*c**2/G},initial speeds |v_1|/c={v_0_1}, |v_2|/c={v_0_2}, with separation r_0 = {r_0} AU, and simulation time: {t_stop/c/Year} years")
	#H = 1.0e-9
	T,Y= bulStoer(F,t_start,X,m_1,m_2,t_stop,H)
	
	x_1 = Y[:,0]
	y_1 = Y[:,1]
	z_1 = Y[:,2]
	
	x_2 = Y[:,3]
	y_2 = Y[:,4]
	z_2 = Y[:,5]
	
	v_x_1 = Y[:,6]
	v_y_1 = Y[:,7]
	v_z_1 = Y[:,8]
	
	v_x_2 = Y[:,9]
	v_y_2 = Y[:,10]
	v_z_2 = Y[:,11]
	
	#Y = None
	
	#CENTER OF MASS
	R_x = (m_1*x_1 + m_2*x_2)/(m_1+m_2)
	R_y = (m_1*y_1 + m_2*y_2)/(m_1+m_2)
	R_z = (m_1*z_1 + m_2*z_2)/(m_1+m_2)
	#separation
	r_sep = np.sqrt((x_1 - x_2)**2 + (y_1 - y_2)**2 + (z_1 - z_2)**2)
	
	#relative velocity
	v_rel_x = v_x_1 - v_x_2
	v_rel_y = v_y_1 - v_y_2
	v_rel_z = v_z_1 - v_z_2
	v_rel = np.sqrt(v_rel_x**2 + v_rel_y**2 + v_rel_z**2)
	i_peryapsis=find_peryapsis(r_sep)
	p_x = x_2-x_1
	p_y = y_2-y_1
	p_z = z_2-z_1
	p_x = p_x[i_peryapsis]
	#print(min(p_x))
	#print(max(p_x))

	p_y = p_y[i_peryapsis]
	p_z = p_z[i_peryapsis]
	T_peryapsis = T[i_peryapsis]
	t_pery = make_delta_T(T_peryapsis)/c/Year*Day
	#Save results
	Q2=np.column_stack((T,x_1,y_1,z_1,v_x_1,v_y_1,v_z_1,x_2,y_2,z_2,v_x_2,v_y_2,v_z_2,R_x,R_y,R_z))
	np.savetxt("plots_and_data/data.txt",Q2,fmt='%.18f',delimiter=',',header='T,x_1,y_1,z_1,v_x_1,v_y_1,v_z_1,x_2,y_2,z_2,v_x_2,v_y_2,v_z_2,R_x,R_y,R_z')
	Q2 = None
	# Now make plots
	fig = plt.figure(figsize=(40, 40))
	fig.patch.set_facecolor('white')
	mpl.rcParams['agg.path.chunksize'] = 10000
	mpl.rcParams.update({'font.size': 40})
	plt.rcParams['text.usetex'] = True
	ax = fig.add_subplot(projection='3d')
	plt.title(r" Center of Mass trajectory ")
	ax.plot(R_x, R_y, R_z, color = "blue", label = r"Center of Mass trajectory")
	ax.scatter(x_1[0], y_1[0], z_1[0],color = "red", marker = 'o', label = r"Star 1 initial postion")
	ax.scatter(x_2[0], y_2[0], z_2[0], color = "green", marker = 'o', label = r"Star 2 initial postion")
	ax.text(x_1[0], y_1[0], z_1[0], r'$m_{1}$ [$M_\odot$] = '+str(m_1*c**2/G), color='red')
	ax.text(x_2[0], y_2[0], z_2[0], r'$m_{2}$ [$M_\odot$] = '+str(m_2*c**2/G), color='green')
	#ax.annotate(,xyz=())
	#ax.annotate(r'$m_{2} [$M_\odot$] = $'+str(m_2),xyz=(x_2[0], y_2[0], z_2[0]))
	ax.set_xlabel(r'$x$ [AU]')
	ax.set_ylabel(r'$y$ [AU]')	
	ax.set_zlabel(r'$z$ [AU]')
	legend= ax.legend(loc='upper right', shadow=True,  markerscale=2, scatterpoints=10,  fontsize=40)
	plt.savefig("plots_and_data/center_of_mass_trajectory.png")
	plt.clf()
	fig = plt.figure(figsize=(40, 40))
	fig.patch.set_facecolor('white')
	mpl.rcParams['agg.path.chunksize'] = 10000
	mpl.rcParams.update({'font.size': 40})
	plt.rcParams['text.usetex'] = True
	ax = fig.add_subplot(projection='3d')
	plt.title(r" Trajectory ")
	ax.plot(x_2, y_2, z_2, color = "green", label = r"First body trajectory")
	ax.plot(x_1, y_1, z_1, color = "red", label = r"Second body trajectory")
	ax.scatter(x_1[0], y_1[0], z_1[0],color = "red", marker = 'o', label = r"Star 1 initial postion")
	ax.scatter(x_2[0], y_2[0], z_2[0], color = "green", marker = 'o', label = r"Star 2 initial postion")
	ax.text(x_1[0], y_1[0], z_1[0], r'$m_{1}$ [$M_\odot$] = '+str(m_1*c**2/G), color='red')
	ax.text(x_2[0], y_2[0], z_2[0], r'$m_{2}$ [$M_\odot$] = '+str(m_2*c**2/G), color='green')
	#ax.annotate(,xyz=())
	#ax.annotate(r'$m_{2} [$M_\odot$] = $'+str(m_2),xyz=(x_2[0], y_2[0], z_2[0]))
	ax.set_xlabel(r'$x$ [AU]')
	ax.set_ylabel(r'$y$ [AU]')	
	ax.set_zlabel(r'$z$ [AU]')
	legend= ax.legend(loc='upper right', shadow=True,  markerscale=2, scatterpoints=10,  fontsize=40)
	plt.savefig("plots_and_data/trajectory.png")
	plt.clf()
	fig = plt.figure(figsize=(40, 40))
	fig.patch.set_facecolor('white')
	mpl.rcParams['agg.path.chunksize'] = 10000
	mpl.rcParams.update({'font.size': 40})
	plt.rcParams['text.usetex'] = True
	ax = fig.add_subplot(projection='3d')
	plt.title(r" Movement of peryapsis during "+str(t_stop/c/Year)+r" years")
	ax.scatter(p_x, p_y, p_z, color = "green", label = r"Position of peryapsis")
	#ax.annotate(,xyz=())
	#ax.annotate(r'$m_{2} [$M_\odot$] = $'+str(m_2),xyz=(x_2[0], y_2[0], z_2[0]))
	ax.set_xlabel(r'$x$ [AU]')
	ax.set_ylabel(r'$y$ [AU]')	
	ax.set_zlabel(r'$z$ [AU]')
	ax.set_xlim(min(p_x), max(p_x))
	legend= ax.legend(loc='upper right', shadow=True,  markerscale=2, scatterpoints=10,  fontsize=40)
	plt.savefig("plots_and_data/peryapsis.png")
	plt.clf()
	
	# fig = plt.figure(figsize=(100,100))
	# fig.patch.set_facecolor('white')
	# mpl.rcParams['agg.path.chunksize'] = 10000
	# mpl.rcParams.update({'font.size': 40})
	# plt.rcParams['text.usetex'] = True
	# plt.title(r"Cumulative shift in periastron shift (s)")
	# #plt.errorbar(L1,rho_g,u_rho_g, linestyle="None")
	# plt.scatter(T_peryapsis, t_pery, lw=3.5, color = "blue", label = r"model")
	# #plt.scatter(phi, rho_g_phi, lw=3.5, color = "red", label = r"$\rho_g$")
	# #plt.plot(logm1, p_m, lw=3.5, ls = ':', color = "cyan", label = r"Model ")
	# #plt.xlim(min(z_cos), max(z_cos))
	# #plt.xscale("log")
	# #plt.yscale("log")
	# ax=plt.gca()
	# ax.tick_params(axis='both', which='major', labelsize=20)
	# ax.tick_params(axis='both', which='minor', labelsize=18)
	# ax.set_xlabel(r'${\Delta}t$ [s]')
	# ax.set_ylabel(r'$T$ [Years]')
	# ax.tick_params(direction='out', length=10, width=2)
	# plt.yticks(fontsize=40)
	# plt.xticks(fontsize=40)	
	# legend= ax.legend(loc='upper right', shadow=True,  markerscale=2, scatterpoints=10,  fontsize=40)
	# plt.savefig("plots_and_data/peri_shift.png")
	# plt.clf()

	
	# fig = plt.figure(figsize=(40,40))
	# fig.patch.set_facecolor('white')
	# mpl.rcParams['agg.path.chunksize'] = 10000
	# mpl.rcParams.update({'font.size': 40})
	# plt.rcParams['text.usetex'] = True
	# plt.title(r"Relative velocities")
	# ax1 = plt.subplot(411)
	# plt.plot(T, v_rel_x,color = "blue", label = r"$\frac{v^x_{rel}}{c}$")
	# plt.tick_params('x', labelsize=6)

	# ax2 = plt.subplot(412, sharex=ax1)
	# plt.plot(T, v_rel_y,color = "red", label = r"$\frac{v^y_{rel}}{c}$")
# # make these tick labels invisible
	# plt.tick_params('x', labelbottom=False)

# # share x and y
	# ax3 = plt.subplot(413, sharex=ax1)
	# plt.plot(T, v_rel_z,color = 'green', label = r"$\frac{v^z_{rel}}{c}$")
	# plt.tick_params('x', labelbottom=False)
	# ax4 = plt.subplot(414, sharex=ax1)
	# plt.plot(T, v_rel_z,color = "purple", label = r"$|\frac{v_{rel}}{c}|$")
	# ax=plt.gca()
	# ax.tick_params(axis='both', which='major', labelsize=20)
	# ax.tick_params(axis='both', which='minor', labelsize=18)
	# ax.set_xlabel(r'$t$ [years]')
	# ax.set_ylabel(r'$r$ [AU]')
	# ax.tick_params(direction='out', length=10, width=2)
	# legend= ax.legend(loc='lower left', shadow=True,  markerscale=2, scatterpoints=10,  fontsize=40)
	# plt.savefig("plots_and_data/relative_velocity.png")
	# plt.clf()