!===================================================================================================================================
! Copyright (c) 2025 GVEC Contributors, Max Planck Institute for Plasma Physics
! License: MIT
!===================================================================================================================================
#include "defines.h"

!===================================================================================================================================
!>
!!# Module ** rProfile **
!!
!! Defines a 1-D profile in rho^2 via B-Spline knots and coefficients
!===================================================================================================================================
MODULE MODgvec_rProfile_bspl
  ! MODULES
  USE MODgvec_Globals ,ONLY: wp, abort
  USE MODgvec_rProfile_base, ONLY: c_rProfile
  USE sll_m_bsplines  ,ONLY: sll_s_bsplines_new, sll_c_bsplines
  IMPLICIT NONE

  PUBLIC

  TYPE, EXTENDS(c_rProfile) :: t_rProfile_bspl

    INTEGER               :: n_knots !! number of knots, including repeated edge knots
    !INTEGER               :: n_coefs !! number of B-Spline coefficients, part of abstract type
    INTEGER               :: deg = 0
    REAL(wp), ALLOCATABLE :: knots(:)   !! knot values, includinng edge knots
    !REAL(wp), ALLOCATABLE :: coefs(:)   !! B-Spline coefficients, part of abstract type
    CLASS(sll_c_bsplines),ALLOCATABLE :: bspl !! b-spline class

    CONTAINS
    PROCEDURE :: eval_at_rho2        => bsplProfile_eval_at_rho2
    PROCEDURE :: antiderivative      => bsplProfile_antiderivative
    FINAL :: bsplProfile_free

  END TYPE t_rProfile_bspl

  INTERFACE t_rProfile_bspl
      MODULE PROCEDURE bsplProfile_new
  END INTERFACE t_rProfile_bspl


  CONTAINS

  !===================================================================================================================================
  !> initialize the rProfile of type bspline
  !!
  !===================================================================================================================================
  FUNCTION bsplProfile_new(knots,  coefs) RESULT(sf)
  ! MODULES
  !-----------------------------------------------------------------------------------------------------------------------------------
  ! INPUT VARIABLES
    REAL(wp),    INTENT(IN) :: knots(:)  !! knots of the B-Spline with repeated start and end points
    REAL(wp),    INTENT(IN) :: coefs(:)  !! B-Spline coefficients
  !-----------------------------------------------------------------------------------------------------------------------------------
  ! OUTPUT VARIABLES
    TYPE(t_rProfile_bspl) :: sf !! self
  !-----------------------------------------------------------------------------------------------------------------------------------
  ! LOCAL VARIABLES
    INTEGER :: n_knots, n_coefs
  !===================================================================================================================================
    n_knots=SIZE(knots)
    n_coefs=SIZE(coefs)
    sf%deg   = COUNT((ABS(knots-knots(1)).LE.1E-12))-1 ! multiplicity of the first knot determines the degree
    IF(COUNT((ABS(knots-knots(n_knots)).LE.1E-12)).NE.sf%deg+1) THEN
      CALL abort(__STAMP__, &
                 "The Bspline knot sequence need the same multiplicity at the beginning and end (=degree+1).")
    END IF
    IF (n_coefs .NE. n_knots-sf%deg-1) THEN
      CALL abort(__STAMP__, &
                 "Number of Bspline coeffcients must be number of knots - (degree+1)!")
    END IF
    sf%n_knots = n_knots
    sf%n_coefs = n_coefs
    ALLOCATE(sf%knots(1:n_knots), sf%coefs(1:n_coefs))
    sf%knots = knots
    sf%coefs = coefs
    IF (sf%deg>0) THEN
      CALL sll_s_bsplines_new(sf%bspl, sf%deg, .FALSE., &
                              sf%knots(1),sf%knots(n_knots),&
                              size(sf%knots(sf%deg+1:n_knots-sf%deg))-1 , & ! number of knots handed to the library
                              sf%knots(sf%deg+1:n_knots-sf%deg)) ! remove repeated edge knots
    END IF
  END FUNCTION bsplProfile_new

  !===================================================================================================================================
  !> evaluate the n-th derivative of the bsplProfile at position s
  !!
  !===================================================================================================================================
  FUNCTION bsplProfile_eval_at_rho2( sf, rho2, deriv ) RESULT(profile_prime_value)
  ! MODULES
  !-----------------------------------------------------------------------------------------------------------------------------------
  ! INPUT VARIABLES
    CLASS(t_rProfile_bspl), INTENT(IN)  :: sf !! self
    REAL(wp)              , INTENT(IN)  :: rho2 !! evaluation point in the toroidal flux coordinate (rho2=phi/phi_edge= rhopos^2)
    INTEGER , OPTIONAL    , INTENT(IN)  :: deriv !! derivative of bspline(rho^2) in rho^2
  !-----------------------------------------------------------------------------------------------------------------------------------
  ! OUTPUT VARIABLES
    REAL(wp)                         :: profile_prime_value
  !-----------------------------------------------------------------------------------------------------------------------------------
  ! LOCAL VARIABLES
    REAL(wp)           :: base_eval(0:sf%deg,0:sf%deg) !! value and derivatives of the (deg+1) B-splines that contribute at s_pos
    INTEGER            :: first_non_zero_bspl !! index offset for the coefficients
    INTEGER             :: deriv_case
  !===================================================================================================================================
    IF (PRESENT(deriv)) THEN
      deriv_case = deriv
    ELSE
      deriv_case = 0
    END IF
    IF(deriv_case.LE.sf%deg)THEN
      CALL sf%bspl%eval_basis_and_n_derivs(rho2,deriv_case,base_eval(0:deriv_case,:),first_non_zero_bspl)
      profile_prime_value =  SUM(sf%coefs(first_non_zero_bspl:first_non_zero_bspl+sf%deg)*base_eval(deriv_case,:))
    ELSE
      profile_prime_value = 0.0_wp
    END IF
  END FUNCTION bsplProfile_eval_at_rho2

  !===================================================================================================================================
  !> get the exact spline antiderivative, with respect to rho2
  !! the knotspan is increased by an extra multiplicity on both ends, and the new coefficients are computed as
  !! beta(i) = beta(i-1) + alpha(i)*(t(i+degree+1)-t(i))/(degree+1)
  !! From deBoor, "A practical guide to Splines", p.128
  !!
  !===================================================================================================================================
  FUNCTION bsplProfile_antiderivative(sf) RESULT(antideriv)
  ! MODULES
  !-----------------------------------------------------------------------------------------------------------------------------------
  ! INPUT VARIABLES
    CLASS(t_rProfile_bspl), INTENT(IN)  :: sf !! self
  !-----------------------------------------------------------------------------------------------------------------------------------
  ! OUTPUT VARIABLES
    CLASS(c_rProfile),ALLOCATABLE :: antideriv
  !-----------------------------------------------------------------------------------------------------------------------------------
  ! LOCAL VARIABLES
    REAL(wp) :: coefs(sf%n_coefs+1), knots(sf%n_knots+2), intermid
    INTEGER :: i, n_coefs, n_knots, deg
  !===================================================================================================================================
    coefs = 0.0_wp
    knots = -42.0_wp

    n_coefs = sf%n_coefs+1
    n_knots = sf%n_knots+2
    deg = sf%deg+1
    ! increase multiplicity at the edges
    knots(2:n_knots-1) = sf%knots
    knots(1) = sf%knots(1)
    knots(n_knots) = sf%knots(sf%n_knots)

    DO i=1,sf%n_coefs
      intermid = sf%coefs(i)*(sf%knots(i+deg)-sf%knots(i))/deg
      coefs(i+1) = coefs(i) + intermid
    END DO
    antideriv = t_rProfile_bspl(knots,  coefs)
  END FUNCTION bsplProfile_antiderivative

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

END MODULE MODgvec_rProfile_bspl