!> \~japanese
!! \brief データセットの補間を行うための手続きを提供するモジュール.
!!
!! \par igmcore_Interpolation モジュールについて
!!
!! \~english
!! \brief This module provides some subroutines for interpolating data.
!!
!! \par About igmcore_Interpolation module
!!
!! \~
!! \since IGMBaseLib 2.0
!!
!! <br><br>
!! \b Copyright (C) GFD Dennou Club, 2011-2012. All rights reserved. <br>
!! \b license ?? <br>
!! \author Yuta Kawai
!!
module igmcore_Interpolation

  !* gtool ***************
  !*

  ! 種類型パラメータ
  ! Kind type parameters
  !
  use dc_types, only: &
    & DP

  !* IGMBaseLib ***********
  !*

  ! 連立一次方程式ソルバ
  ! Simultaneous equation solver
  !
  use igmcore_SimLinearAlgebraEq, only: &
    & LUdecomp_TriDiagMat, sol_SimLinearAlEq_lutri

  ! 宣言文 ; Declaration statements
  !
  implicit none
  private

  type, public :: CubicSpline_SData
    real(DP) :: a, b, c, d, refPos
  end type CubicSpline_SData

  public :: CubicSpline_Interpolation, CubicSpline_eval

contains
function CubicSpline_Interpolation( &
  & grid, data &
  & ) result(S)

  real(DP), intent(in) :: grid(:)
  real(DP), intent(in) :: data(:)
  type(CubicSpline_SData) :: S(size(grid) - 1)

  real(DP) :: h(size(S))
  real(DP) :: u(size(grid))
  real(DP) :: al(size(grid))
  real(DP) :: ad(size(grid))
  real(DP) :: au(size(grid))
  real(DP) :: RHS(size(grid))

  integer :: S_size, grid_size, rhs_size
  integer :: j

  S_size = size(S)
  grid_size = size(grid)
  rhs_size = size(RHS)

  h(1:S_size) = grid(2:S_size+1) - grid(1:S_size)

  do j=2, grid_size-1
    al(j) = h(j-1)
    ad(j) = 2d0 * ( h(j-1) + h(j) )
    au(j) = h(j)

    RHS(j) = 6.0d0 * ( ( data(j+1) - data(j) ) / h(j) - ( data(j) - data(j-1) ) / h(j-1) )
  end do

  u(1) = 0d0
  u(grid_size) = 0d0

  call LUdecomp_TriDiagMat(al(2:grid_size-1), ad(2:grid_size-1), au(2:grid_size-1))
  u(2:grid_size-1) = sol_SimLinearAlEq_lutri(al(2:grid_size-1), ad(2:grid_size-1), au(2:grid_size-1), RHS(2:grid_size-1))

  do j=1, S_size
    S(j)%a = ( u(j+1) - u(j) ) / ( 6.0d0 * h(j) )
    S(j)%b = 0.5d0 * u(j)
    S(j)%c = ( data(j+1) - data(j) ) /  h(j) - h(j) * ( 2d0 * u(j) + u(j+1) ) / 6d0
    S(j)%d = data(j)
    S(j)%refPos = grid(j)
  end do

end function CubicSpline_Interpolation

function CubicSpline_eval(  &
  & csSData, evalPos &
  & ) result(val)

  type(CubicSpline_SData), intent(in) :: csSData
  real(DP), intent(in) :: evalPos
  real(DP) :: val

  real(DP) :: DelX

  DelX = evalPos - csSData%refPos
  val = ( ( csSData%a * DelX + csSData%b ) * DelX + csSData%c ) * DelX + csSData%d
end function CubicSpline_eval

end module igmcore_Interpolation
