# 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 numpy as np
#import zeros,float,sum
from math import sqrt
from numba import jit, cuda, njit
# 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 midpoint(F,x,y,a,b,xStop,nSteps):
# Midpoint formulas
	h = (xStop - x)/nSteps
	y0 = y
	y1 = y0 + h*F(x,y0,a,b)
	for i in range(nSteps-1):
		x = x + h
		y2 = y0 + 2.0*h*F(x,y1,a,b)
		y0 = y1
		y1 = y2
	return(0.5*(y1 + y0 + h*F(x,y2,a,b)))
@jit(target_backend='cuda',nopython=True)
def richardson(r,k):
# Richardson’s extrapolation
	for j in range(k-1,0,-1):
		const = (k/(k - 1.0))**(2.0*(k-j))
		r[j] = (const*r[j+1] - r[j])/(const - 1.0)
	return
@jit(target_backend='cuda',nopython=True)
def integrate(F,x,y,a,b,xStop,tol):

	kMax = 51
	n = len(y)
	r = np.zeros((kMax,n),dtype=np.float64)
# Start with two integration steps
	nSteps = 2
	r[1] = midpoint(F,x,y,a,b,xStop,nSteps)
	r_old = r[1].copy()
# Increase the number of integration points by 2
# and refine result by Richardson extrapolation
	for k in range(2,kMax):
		nSteps = 2*k
		r[k] = midpoint(F,x,y,a,b,xStop,nSteps)
		richardson(r,k)
# Compute RMS change in solution
		e = np.sqrt(np.sum((r[1] - r_old)**2)/n)
# Check for convergence
		if e < tol: 
			return r[1]
		r_old = r[1].copy()
	print("Midpoint method did not converge")