# 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/>

#from numpy import array,sqrt
from midpoint import *
from numba import jit, cuda, njit
import numpy as np
# Loosely based on Jaan Kiusalaas. 2010. Numerical Methods in Engineering with Python (2nd. ed.). Cambridge University Press, USA.
#@jit(target_backend='cuda',nopython=True)
def bulStoer(F,x,y,a,b,xStop,H,tol=1.0e-11):
	X = []
	Y = []
	X.append(x)
	Y.append(y)
	R_x = (a*y[0] + b*y[3])/(a+b)
	R_y = (a*y[1] + b*y[4])/(a+b)
	R_z = (a*y[2] + b*y[5])/(a+b)
	R_1 = np.sqrt((R_x-y[0])**2 + (R_y-y[1])**2 + (R_z-y[2])**2)
	R_2 = np.sqrt((R_x-y[3])**2 + (R_y-y[4])**2 + (R_z-y[5])**2)

	z_1 = ((a+b)*np.sqrt(y[6]**2+y[7]**2+y[8]**2)/R_1)**(2/3)
	z_2 = ((a+b)*np.sqrt(y[10]**2+y[11]**2+y[9]**2)/R_2)**(2/3)
	nu = a*b/(a+b)**2
	z_isco = 3/14/nu*(1-np.sqrt(1-14*nu/9))#3/14/nu*(1-np.sqrt(1-14*nu/9))
	r=np.sqrt((y[0]-y[3])**2+(y[1]-y[4])**2+(y[2]-y[5])**2)
	while x < xStop:
		if z_1+z_2 >=z_isco and r<=12*(a+b):
			break
		#H = min(H,xStop - x)
		H = min(H,xStop - x)
		y = integrate(F,x,y,a,b,x + H,tol) # Midpoint method
		x = x + H
		X.append(x)
		Y.append(y)# Midpoint method
		#y1 = make_list_to_array(y)
		#Y = np.column_stack((Y, y1))
		# Center of mass frame
		R_x = (a*y[0] + b*y[3])/(a+b)
		R_y = (a*y[1] + b*y[4])/(a+b)
		R_z = (a*y[2] + b*y[5])/(a+b)
		R_1 = np.sqrt((R_x-y[0])**2 + (R_y-y[1])**2 + (R_z-y[2])**2)
		R_2 = np.sqrt((R_x-y[3])**2 + (R_y-y[4])**2 + (R_z-y[5])**2)

		z_1 = ((a+b)*np.sqrt(y[6]**2+y[7]**2+y[8]**2)/R_1)**(2/3)
		z_2 = ((a+b)*np.sqrt(y[10]**2+y[11]**2+y[9]**2)/R_2)**(2/3)
		#nu = a*b/(a+b)**2
		r = np.sqrt((y[0]-y[3])**2+(y[1]-y[4])**2+(y[2]-y[5])**2)
		#print(x,xStop)
	#print(Y)
	return np.array(X),np.array(Y)