! -*- mode: f90; coding: utf-8 -*-
!----------------------------------------------------------------------
! Copyright (c) 2019 SPMODEL Development Group. All rights reserved.
!----------------------------------------------------------------------
!> @author Shin-ichi Takehiro, Youhei SASAKI
!> @copyright Copyright (C) SPMODEL Development Group, 2019.
!>            License is MIT/X11. see [COPYRIGHT](@ref COPYRIGHT) in detail
!> @brief
!> @ref w_module の下部モジュール: 補間計算
!> @note
!> 補間計算の方法については doc/w_module.tex を参照のこと.
!> @warning
!> * 本モジュールを直接呼ぶ事は想定されていないことに注意されたい．
!>   ユーザは @ref w_module を呼ぶこと．
!> * 変換する格子点データ, スペクトルデータの配列の大きさは決めうちである
!----------------------------------------------------------------------
!履歴  2019/03/16  竹広真一
!      2019/10/04  佐々木洋平
!----------------------------------------------------------------------
module w_module_interpolate_mint
  use dc_types, only: DP
  use dc_message, only : MessageNotify
  use spml_utils, only: gammaln
  use w_module_base_mint, only : mm,nn, mi, l_nm
  implicit none
  private
  public Interpolate_w

  interface Interpolate_w
    module procedure Interpolate_array00_w
  end interface Interpolate_w

  interface alpha
    module procedure alpha_array0
  end interface alpha

  interface Pmm
    module procedure Pmm_array0
  end interface Pmm

contains

  !----------------------------------------------------------
  !> 緯度 alat, 経度 alon における関数値を
  !> その球面調和変換係数 w_data から補間計算する
  !----------------------------------------------------------
  function Interpolate_array00_w(w_data,alon,alat)
    !> スペクトルデータ
    real(DP), intent(IN) :: w_data((mm/mi+1)*(2*(nn+1)-mm/mi*mi))
    !> 補間する位置(経度)
    real(DP), intent(IN) :: alon
    !> 補間する位置(緯度)
    real(DP), intent(IN) :: alat
    !> 補間した値
    real(DP)             :: Interpolate_array00_w

    real(DP)   :: mu
    real(DP)   :: y0
    real(DP)   :: y1
    real(DP)   :: y2
    real(DP)   :: AnmPnm
    integer(8) :: k
    integer(8) :: m

    mu = sin(alat)
    Interpolate_array00_w = 0.0D0

    !---- Σa_n^0 P_n^0 の計算
    y2 = 0.0D0
    y1 = 0.0D0
    do k=nn,1,-1
      y0 =   alpha(k,0_8,mu) * y1  &
        &  + beta(k+1_8,0_8) * y2  &
        &  + w_data(l_nm(int(k),0))
      y2 = y1
      y1 = y0
    enddo
    Interpolate_array00_w =  &
      & ( beta(1_8,0_8) * y2 + mu*sqrt(3.0D0) * y1 &
      &   + w_data(l_nm(0,0))  ) * Pmm(0_8,mu)

    !----  Σ Re[s_n^m] P_n^m exp(imλ)の計算
    do m=mi,mm,mi
      y2 = 0.0D0
      y1 = 0.0D0
      do k=nn,m+1,-1
        y0 =   alpha(k,m,mu) * y1  &
          &  + beta(k+1_8,m) * y2  &
          &  + w_data(l_nm(int(k),int(m)))
        y2 = y1
        y1 = y0
      enddo

      AnmPnm = ( w_data(l_nm(int(m),int(m))) &
        &        + beta(m+1_8,m)*y2 &
        &        + mu*sqrt(2.0D0*m+3) * y1 ) * Pmm(m,mu)

      Interpolate_array00_w = Interpolate_array00_w &
        &                     + AnmPnm * 2 * cos(m*alon)
      end do

      !----  Σ Im[s_n^m] P_n^m exp(imλ)の計算
      do m=mi,mm,mi
        y2 = 0.0D0
        y1 = 0.0D0
        do k=nn,m+1,-1
          y0 =   alpha(k,m,mu) * y1 &
            &  + beta(k+1_8,m) * y2 &
            &  + w_data(l_nm(int(k),-int(m)))
          y2 = y1
          y1 = y0
        enddo

        AnmPnm =( w_data(l_nm(int(m),-int(m)))  &
          &       + beta(m+1_8,m) * y2          &
          &       + mu*dsqrt(2.0D0*m+3) * y1 ) * Pmm(m,mu)


        Interpolate_array00_w = Interpolate_array00_w &
          &                     - AnmPnm * 2 * sin(m*alon)
      end do

    end function Interpolate_array00_w

    !----------------------------------------------------------
    !> 漸化式の \f$P_{n}^{m}\f$ の係数
    !----------------------------------------------------------
    function alpha_array0(n,m,x)
      integer(8), intent(IN)  :: n
      integer(8), intent(IN)  :: m
      real(DP), intent(IN) :: x
      real(DP)             :: alpha_array0

      alpha_array0 = dsqrt( (2.0D0*n+3)*(2.0D0*n+1)/((n-m+1)*(n+m+1)) ) * x
    end function alpha_array0

    !----------------------------------------------------------
    !> 漸化式の \f$P_{n-1}^m\f$ の係数
    !----------------------------------------------------------
    function beta(n,m)
      integer(8), intent(IN) :: n
      integer(8), intent(IN) :: m
      real(DP)               :: beta

      beta = - dsqrt( (2.0D0*n+3)*(n+m)*(n-m)/((2*n-1)*(n+m+1)*(n-m+1)) )
    end function beta

    !----------------------------------------------------------
    !> ルジャンドル陪函数 P_m^m(x) の計算
    !>
    !>         sqrt( factrl(2*m+1) )/(2.0D0**m * factrl(m)) &
    !>          * (1-x**2)**(m/2.0D0)
    !>
    !> 係数が発散しないように対数で計算する
    !----------------------------------------------------------
    function Pmm_array0(m,x)
      !> ルジャンドル陪函数の次数
      integer(8),  intent(IN) :: m
      !> 引き数
      real(DP), intent(IN) :: x
      !> ルジャンドル陪函数の値
      real(DP)             :: Pmm_array0

      if ( m < 0 ) call MessageNotify('E', &
        & 'Pmm in w_Intepolate_module',    &
        & 'order m should be larger equal to zero')

      Pmm_array0 = 1.0D0
      if ( m > 0 )then
        Pmm_array0 = exp(gammaln(2.0D0*m+2.0D0)/2.0D0        &
          &              - m*dlog(2.0D0) - gammaln(m+1.0D0)) &
          &          *(1-x**2)**(m/2.0D0)
      endif
    end function Pmm_array0

end module w_module_interpolate_mint