!===================================================================================================================================
! Copyright (C) 2025 -  Florian Hindenlang <hindenlang@gmail.com>
!
! This file is part of GVEC. GVEC 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.
!
! GVEC 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 v3.0 for more details.
!
! You should have received a copy of the GNU General Public License along with GVEC. If not, see <http://www.gnu.org/licenses/>.
!=================================================================================================================================
#include "defines.h"

!===================================================================================================================================
!>
!!# Module ** cubic spline **
!!
!! Define a cubic spline by points and function values, with boundary conditions and represent it as a B-Spline for evaluation
!!
!===================================================================================================================================
MODULE MODgvec_cubic_spline
  ! MODULES
  USE MODgvec_globals, ONLY: wp,abort
  use sll_m_bsplines, only: sll_c_bsplines
  IMPLICIT NONE

  PUBLIC

  TYPE :: t_cubspl
    !--------------------------------------------------------------------------------------------------------------------------------
    REAL(wp),ALLOCATABLE  :: coefs(:)   !! B-Spline coefficients
    REAL(wp),ALLOCATABLE  :: knots(:)   !! knots (=break points with repeated knots at the end)
    CLASS(sll_c_bsplines),ALLOCATABLE :: bspl !! b-spline class
    CONTAINS

    PROCEDURE :: eval  => cubspl_eval
    PROCEDURE :: eval_at_rho  => cubspl_eval  !! for testing
    FINAL :: cubspl_free

  END TYPE t_cubspl

  INTERFACE t_cubspl
    MODULE PROCEDURE cubspl_new
  END INTERFACE t_cubspl

  INTERFACE interpolate_cubic_spline
    MODULE PROCEDURE interpolate_cubic_spline
  END INTERFACE

  CONTAINS

  !===================================================================================================================================
  !> Interpolation of function values f(x_i)=f_i, i=1,n with a cubic spline, given left and right boundary condition
  !! types of boundary conditions:
  !! 0: not-a-knot
  !! 1: f'(x_boundary)=0
  !! 2: f''(x_boundary)=0
  !!
  !===============================================================================================================================
  FUNCTION cubspl_new(x,f,BC,BC_val) RESULT(sf)
    ! MODULES
    USE sll_m_bsplines, ONLY: sll_s_bsplines_new
    IMPLICIT NONE
    !-----------------------------------------------------------------------------------------------------------------------------------
    ! INPUT VARIABLES
    REAL(wp) , INTENT(IN) :: x(:) !! x positions
    REAL(wp) , INTENT(IN) :: f(:) !! function values at x positions
    INTEGER  , INTENT(IN) :: BC(1:2) !! Boundary condition at x(1)/x(n): =0: not-a-knot, =1: first der. =BC_val(1)/BC_val(2), =2: second der. =BC_val(1)/BC_val(2)
    REAL(wp) , INTENT(IN),OPTIONAL :: BC_val(1:2) !! Boundary value for BC(1:2) >0,
    !-----------------------------------------------------------------------------------------------------------------------------------
    ! OUTPUT VARIABLES
    TYPE(t_cubspl) :: sf !! self
    !-----------------------------------------------------------------------------------------------------------------------------------
    ! LOCAL VARIABLES
    INTEGER, PARAMETER:: deg=3 !! degree of the spline
    INTEGER :: nbreaks
    !===================================================================================================================================

    CALL interpolate_cubic_spline(x,f,sf%coefs,sf%knots,BC,BC_val)
    nbreaks=SIZE(sf%knots)-2*deg
    CALL sll_s_bsplines_new(sf%bspl, deg, .FALSE., sf%knots(1+deg),sf%knots(nbreaks+deg),nbreaks-1,sf%knots(1+deg:nbreaks+deg))
  END FUNCTION cubspl_new


  SUBROUTINE interpolate_cubic_spline(x,f,coefs,knots,BC,BC_val)
    ! MODULES
    USE sll_m_bsplines, ONLY: sll_s_bsplines_new
    USE sll_m_spline_matrix,ONLY: sll_c_spline_matrix,sll_s_spline_matrix_new
    IMPLICIT NONE
    !-----------------------------------------------------------------------------------------------------------------------------------
    ! INPUT VARIABLES

    REAL(wp) , INTENT(IN) :: x(:) !! x positions
    REAL(wp) , INTENT(IN) :: f(:) !! function values at x positions
    INTEGER  , INTENT(IN) :: BC(1:2) !! Boundary condition at x(1)/x(n): =0: not-a-knot, =1: first der. =BC_val(1)/BC_val(2), =2: second der. =BC_val(1)/BC_val(2)
    REAL(wp) , INTENT(IN),OPTIONAL :: BC_val(1:2) !! Boundary value for BC(1:2) >0,
    !-----------------------------------------------------------------------------------------------------------------------------------
    ! OUTPUT VARIABLES
    REAL(wp),ALLOCATABLE,INTENT(INOUT) :: coefs(:)  !! B-Spline coefficients of interpolated cubic spline
    REAL(wp),ALLOCATABLE,INTENT(INOUT) :: knots(:)  !! B-Spline knots of interpolated cubic spline
    !-----------------------------------------------------------------------------------------------------------------------------------
    ! LOCAL VARIABLES
    INTEGER :: innerpts(1:2),nbc_left,nbc_right,i,jmin,s,n,nbreaks,ncoefs
    INTEGER,PARAMETER :: deg=3
    REAL(wp):: base(0:deg),deriv_values(0:deg,0:deg)
    CLASS(sll_c_bsplines),ALLOCATABLE :: bspl
    CLASS(sll_c_spline_matrix),ALLOCATABLE :: solvemat
    !===================================================================================================================================
    n=SIZE(x)
    IF(SIZE(f).NE.n) THEN
      CALL abort(__STAMP__, &
                 "x and f must have the same size for cubic spline interpolation")
    END IF
    innerpts=(/2,n-1/)
    nbreaks=n
    ncoefs=n
    SELECT CASE(BC(1))
    CASE(0) !not-a-knot: leave out second knot
      nbreaks=nbreaks-1
      innerpts(1)=3
      nbc_left=0
    CASE(1,2)
      ncoefs=ncoefs+1
      nbc_left=1
    CASE DEFAULT
      CALL abort(__STAMP__, &
                 "BC_left not implemented")
    END SELECT
    SELECT CASE(BC(2))
    CASE(0) !not-a-knot: leave out second to last knot
      nbreaks=nbreaks-1
      innerpts(2)=n-2
      nbc_right=0
    CASE(1,2)
      ncoefs=ncoefs+1
      nbc_right=1
    CASE DEFAULT
      CALL abort(__STAMP__, &
                 "BC_right not implemented")
    END SELECT
    IF(n+nbc_left+nbc_right.LT.4) THEN
      CALL abort(__STAMP__, &
                 "minimum number of conditions of a cubic spline is 4")
    END IF
    ALLOCATE(coefs(ncoefs),knots(1:nbreaks+2*deg))
    knots = -42e5
    knots(1:1+deg)=x(1) !repreat first knot deg+1 times
    knots(nbreaks+deg:nbreaks+2*deg)=x(n) !repeat last knot deg+1 times
    IF(innerpts(2)-innerpts(1)+1.GT.0)THEN
      knots(2+deg:nbreaks+deg-1)=x(innerpts(1):innerpts(2))
    END IF
    CALL sll_s_bsplines_new(bspl, deg, .FALSE., x(1),x(n),nbreaks-1,knots(1+deg:nbreaks+deg))

    IF(bspl%nBasis.NE.ncoefs) CALL abort(__STAMP__, &
                            'problem with bspl basis in cubic spline')
    CALL sll_s_spline_matrix_new(solvemat , "banded",ncoefs,deg,deg)
    ! Interpolation points
    do i = 1, n
      call bspl % eval_basis( x(i), base, jmin )
      coefs(i+nbc_left)= f(i)  !! used as rhs, then overwritten in solve
      do s = 0, deg
        call solvemat % set_element( i+nbc_left, jmin+s, base(s) )
      end do
    end do

    ! Boundary conditions
    IF(BC(1).GT.0)THEN !deriv=BC(1)
      CALL bspl%eval_basis_and_n_derivs(x(1),BC(1),deriv_values(0:BC(1),0:deg),jmin)
      do s = 0, deg
        call solvemat % set_element( 1, jmin+s, deriv_values(BC(1),s) )
      end do
      IF(PRESENT(BC_val))THEN
        coefs(1)=BC_val(1)! used as rhs, then overwritten in solve
      ELSE
        coefs(1)=0.0_wp ! used as rhs, then overwritten in solve
      END IF
    END IF
    IF(BC(2).GT.0)THEN !deriv=BC(2)
      CALL bspl%eval_basis_and_n_derivs(x(n),BC(2),deriv_values(0:BC(2),0:deg),jmin)
      do s = 0, deg
        call solvemat % set_element( ncoefs, jmin+s, deriv_values(BC(2),s) )
      end do
      IF(PRESENT(BC_val))THEN
        coefs(ncoefs)=BC_val(2)! used as rhs, then overwritten in solve
      ELSE
        coefs(ncoefs)=0.0_wp ! used as rhs, then overwritten in solve
      END IF
    END IF
    CALL solvemat % factorize()
    CALL solvemat%solve_inplace(1,coefs)

  END SUBROUTINE interpolate_cubic_spline

!===================================================================================================================================
!> evaluate the n-th derivative of the bsplProfile at position s
!!
!===================================================================================================================================
  FUNCTION cubspl_eval( sf, xpos, deriv ) RESULT(y)
    ! MODULES
    !-----------------------------------------------------------------------------------------------------------------------------------
    ! INPUT VARIABLES
      CLASS(t_cubspl), INTENT(IN)  :: sf !! self
      REAL(wp)       , INTENT(IN)  :: xpos(:) !! position
      INTEGER        , INTENT(IN)  :: deriv !! derivative (=0: no derivative)

    !-----------------------------------------------------------------------------------------------------------------------------------
    ! OUTPUT VARIABLES
      REAL(wp)                     :: y(size(xpos))
    !-----------------------------------------------------------------------------------------------------------------------------------
    ! LOCAL VARIABLES
      REAL(wp) :: deriv_values(0:deriv,0:sf%bspl%degree) !! values of the (deg+1) B-splines that contribute at s_pos
      INTEGER  :: i,first_non_zero_bspl !! index offset for the coefficients
    !===================================================================================================================================
      IF(deriv.EQ.0)THEN
        DO i=1,size(xpos)
          CALL sf%bspl%eval_basis(xpos(i),deriv_values(0,0:sf%bspl%degree),first_non_zero_bspl)
          y(i) =  SUM(sf%coefs(first_non_zero_bspl:first_non_zero_bspl+sf%bspl%degree)*deriv_values(0,0:sf%bspl%degree))
        END DO
      ELSEIF(deriv.LE.sf%bspl%degree) THEN
        DO i=1,size(xpos)
          CALL sf%bspl%eval_basis_and_n_derivs(xpos(i),deriv,deriv_values,first_non_zero_bspl)
          y(i) =  SUM(sf%coefs(first_non_zero_bspl:first_non_zero_bspl+sf%bspl%degree)*deriv_values(deriv,0:sf%bspl%degree))
        END DO
      ELSE
        y = 0.0_wp
      END IF
    END FUNCTION cubspl_eval

    !===================================================================================================================================
    !> finalize the type rProfile
    !!
    !===================================================================================================================================
    SUBROUTINE cubspl_free(sf)
    ! MODULES
    !-----------------------------------------------------------------------------------------------------------------------------------
    ! INPUT VARIABLES
    !-----------------------------------------------------------------------------------------------------------------------------------
    ! OUTPUT VARIABLES
      TYPE(t_cubspl), INTENT(INOUT) :: sf !! self
    !-----------------------------------------------------------------------------------------------------------------------------------
    !===================================================================================================================================
      IF (ALLOCATED(sf%bspl)) CALL sf%bspl%free()
    END SUBROUTINE cubspl_free


END MODULE MODgvec_cubic_spline