! -*- 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 1 次元周期境界: 実フーリエ変換
!>
!> @details
!> 1 次元周期境界条件の下での流体運動を
!> 実フーリエ変換によるスペクトル法で数値計算するための
!> Fortran90 関数を提供する.
!> 2 次元データの 1 次元に関して同時にスペクトル計算を実行するための
!> 関数も提供しており, 2, 3 次元領域での計算のベースも提供する.
!>
!> このモジュールは内部で ISPACK3/FXPACK のサブルーチンを呼んでいる.
!> スペクトルデータの格納方法については ISPACK3/FXPACK
!> と異なっているので以下のコメントに注意されたい.
!>
!> ## 関数・変数の名前と型について
!>
!> ### 命名法
!>
!> * 関数名の先頭 (`e_`, `g_`, `ae_`, `ag_`) は, 返す値の形を示している.
!>   * `e_`  :: スペクトルデータ,
!>   * `g_`  :: 1 次元格子点データ,
!>   * `ae_` :: 1 次元スペクトルデータが複数並んだ 2 次元データ,
!>   * `ag_` :: 1 次元格子点データが複数並んだ 2 次元データ.
!>
!> * 関数名の間の文字列(`Dx`)は, その関数の作用を表している.
!>
!> * 関数名の最後 (`_e`, `_ae`, `_g`, `_ag`) は,
!>   入力変数の形スペクトルデータおよび格子点データであることを示している.
!>   * `_e`  :: スペクトルデータ
!>   * `_g`  :: 1 次元格子点データ
!>   * `_ae` :: 1 次元スペクトルデータが複数並んだ 2 次元データ
!>   * `_ag` :: 1 次元格子点データが複数並んだ 2 次元データ
!>
!> ### 各データの種類の説明
!>
!> * 以下 `DP  = kind(1.0D)` である
!>
!> * `g` : 1 次元格子点データ.
!>   * 変数の種類と次元は `real(DP), dimension(0:im-1)`.
!>   * `im` は X 座標の格子点数であり, サブルーチン ae_initial() にて
!>     あらかじめ設定しておく.
!>
!> * `e` : スペクトルデータ.
!>   * 変数の種類と次元は `real(DP), dimension(-km:km)`.
!>   * `km` は X 方向の最大波数であり, サブルーチン ae_initial() にて
!>     あらかじめ設定しておく.
!>   * スペクトルデータの格納のされ方は
!>     `0:km` までが \f$\cos(kx)\f$ の係数,
!>     `-km:-1` が  \f$sin(kx)\f$ の係数となっている.
!>
!> * `ag` : 1 次元(X)格子点データの並んだ 2 次元データ.
!>   * 変数の種類と次元は `real(DP), dimension(:,0:im-1)`.
!>     第 2 次元が X 方向を表す.
!>
!> * `ae` : 1 次元スペクトルデータの並んだ 2 次元データ.
!>   * 変数の種類と次元は `real(DP), dimension(:,-km:km)`.
!>     第 2 次元がスペクトルを表す.
!>
!> * `g_` で始まる関数が返す値は 1 次元格子点データに同じ.
!>
!> * `e_` で始まる関数が返す値はスペクトルデータに同じ.
!>
!> * `ag_` で始まる関数が返す値は 1 次元格子点データの並んだ
!>   2 次元データに同じ.
!>
!> * `ae_` で始まる関数が返す値は 1 次元スペクトルデータの並んだ
!>   2 次元データに同じ.
!>
!> * スペクトルデータに対する微分等の作用とは, 対応する格子点データに
!>   微分などを作用させたデータをスペクトル変換したもののことである.
!>
module ae_module
  use dc_message, only: MessageNotify
  use dc_types, only: DP
  use spml_utils, only: pi
  implicit none

  private
  public ae_Initial
  public ag_ae
  public ae_ag
  public g_e
  public e_g
  public ae_Dx_ae
  public e_Dx_e
  public a_Int_ag
  public Int_g
  public a_Avr_ag
  public Avr_g
  public Interpolate_e
  public a_Interpolate_ae
  public g_X
  public g_X_Weight

  !> 格子点座標(X)を格納した1次元配列.
  real(DP), allocatable    :: g_x(:)
  !> 重み座標を格納した 1 次元配列. X方向の格子点の間隔が格納してある.
  real(DP), allocatable    :: g_x_weight(:)

  integer(8)               :: im = 32    ! 格子点の数
  integer(8)               :: km = 10    ! 切断波数
  real(DP)                 :: xl = 1.0D0 ! 領域の大きさ
  integer(8), allocatable  :: it(:)
  real(DP),   allocatable  :: t(:)

  save im, km, it, t, xl, g_X, g_X_Weight

contains
  !-----------------------------------------------------------------------
  !> スペクトル変換の格子点数, 波数, 領域の大きさを設定する.
  !>
  !> 他の関数や変数を呼ぶ前に, 最初にこのサブルーチンを呼んで
  !> 初期設定をしなければならない.
  !>
  subroutine ae_Initial(i,k,xmin,xmax)
    !> 格子点の数
    integer, intent(in) :: i
    !> 切断波数
    integer, intent(in) :: k
    !> X座標下端
    real(DP),intent(in) :: xmin
    !> X座標上端
    real(DP),intent(in) :: xmax
    integer(8)          :: ii

    im = i
    km = k
    xl = xmax - xmin

    if ( im <= 0 .or. km <= 0 ) then
       call MessageNotify('E','ae_Initial', &
            'Number of grid points and waves should be positive')
    elseif ( mod(im,2_8) /= 0 ) then
       call MessageNotify('E','ae_Initial', &
            'Number of grid points should be even')
    elseif ( km > im/2 ) then
       call MessageNotify('E','ae_Initial', &
            'KM shoud be less equal to IM/2')
    endif

    allocate(it(im/2),t(im*3/2))

    call fxrini(im,it,t)

    allocate(g_x(0:im-1))
    do ii=0,im-1
       g_X(ii) = xmin + (xl/im)*ii
    enddo

    allocate(g_x_weight(0:im-1))
    g_X_Weight = xl/im

    call MessageNotify(&
         'M','ae_initial','ae_module (2019/09/17) is initialized')

  end subroutine ae_Initial

  !-----------------------------------------------------------------------
  !> スペクトルデータから格子点データへ逆変換する(2 次元データ用)
  !>
  !> スペクトル逆変換の定義は以下のとおり:
  !> * スペクトルデータ: \f$e_k \, (k=-km,...,km)\f$
  !> * 格子点データ:     \f$g_j \, (j=-0,...,im-1)\f$
  !> \f{align*}{
  !>  g_j = e_0
  !>      + 2\sum_{k=1}^{km}e_k\cos(2\pi jk/im) - e_{-k}\sin(2\pi jk/im))
  !>                                               \quad (j=0,1,...,im-1).
  !> \f}
  function ag_ae(ae)
    !> スペクトルデータ
    real(DP), dimension(:,-km:), intent(in)  :: ae
    !> 格子点データ
    real(DP), dimension(size(ae,1),0:im-1)   :: ag_ae

    integer(8) :: nm, k
    nm=size(ae,1)

    if ( size(ae,2) < 2*km+1 ) then
       call MessageNotify('E','ag_ae', &
            'The Fourier dimension of input data too small.')
    elseif ( size(ae,2) > 2*km+1 ) then
       call MessageNotify('W','ag_ae', &
            'The Fourier dimension of input data too large.')
    endif

    ag_ae       = 0.0D0
    ag_ae(:,0)  = ae(:,0)
    do k=1,min(km,im/2-1)
       ag_ae(:,2*k)   = ae(:,k)
       ag_ae(:,2*k+1) = ae(:,-k)
    enddo
    if ( km == im/2 ) then
      ag_ae(:,1) = ae(:,km)
    end if

    call fxrtba(nm,im,ag_ae,it,t)

  end function ag_ae

  !-----------------------------------------------------------------------
  !> スペクトルデータから格子点データへ逆変換する(1 次元データ用)
  !>
  !>  * スペクトルデータ: \f$e_k \, (k=-km,...,km)\f$
  !>  * 格子点データ:  \f$g_j \, (j=-0,...,im-1)\f$
  !> \f{align*}{
  !>  g_j = e_0
  !>      + 2\sum_{k=1}^{km}(e_k\cos(2\pi jk/im) - e_{-k}\sin(2\pi jk/im))
  !>                                                  \quad (j=0,1,...,im-1).
  !> \f}
  function g_e(e)
    !> スペクトルデータ
    real(DP), dimension(-km:km), intent(in) :: e
    !> 格子点データ
    real(DP), dimension(0:im-1)             :: g_e

    real(DP), dimension(1,size(e))  :: ae_work
    real(DP), dimension(1,0:im-1)   :: ag_work

    ae_work(1,:) = e
    ag_work = ag_ae(ae_work)
    g_e = ag_work(1,:)

  end function g_e

  !-----------------------------------------------------------------------
  !> 格子点データからスペクトルデータへ正変換する(2 次元データ用)
  !>
  !> スペクトル正変換の定義は以下のとおり:
  !> * 格子点データ     \f$g_j \, (j=-0,...,im-1)\f$
  !> * スペクトルデータ \f$e_k \, (k=-km,...,km)\f$
  !> \f{align*}{
  !>   e_0 & = (1/im)\sum_{j=0}^{im-1} g_j \\
  !>   e_k & = (1/im)\sum_{j=0}^{im-1} g_j \cos(2\pi jk/im)
  !>          \quad (k=1,2,...,km) \\
  !>   e_{-k} & = - (1/im)\sum_{j=0}^{im-1} g_j \sin(2\pi jk/im)
  !>          \quad (k=1,2,...,km)
  !> \f}
  function ae_ag(ag)
    !> 格子点データ
    real(DP), dimension(:,:), intent(in)   :: ag
    !> スペクトルデータ
    real(DP), dimension(size(ag,1),-km:km) :: ae_ag

    real(DP), dimension(size(ag,1),0:im-1) :: ag_work
    integer(8) :: nm, k

    nm = size(ag,1)
    if ( size(ag,2) < im ) then
       call MessageNotify('E','ae_ag', &
            'The Grid points of input data too small.')
    elseif ( size(ag,2) > im ) then
       call MessageNotify('W','ae_ag', &
            'The Grid points of input data too large.')
    endif
    ag_work = ag

    call fxrtfa(nm,im,ag_work,it,t)

    do k=1,min(km,im/2-1)
       ae_ag(:, k) = ag_work(:,2*k)
       ae_ag(:,-k) = ag_work(:,2*k+1)
    enddo
    ae_ag(:,0) = ag_work(:,0)
    if ( km == im/2 ) then
      ae_ag(:,km) = ag_work(:,1)
    end if

  end function ae_ag

  !-----------------------------------------------------------------------
  !> 格子点データからスペクトルデータへ正変換する(1 次元データ用)
  !>
  !> スペクトル正変換の定義は以下のとおり:
  !> * 格子点データ:     \f$g_j \, (j=-0,...,im-1)\f$
  !> * スペクトルデータ: \f$e_k \, (k=-km,...,km)\f$
  !> \f{align*}{
  !>  e_0    & = (1/im)\sum_{j=0}^{im-1} g_j \\
  !>  e_k    & = (1/im)\sum_{j=0}^{im-1} g_j \cos(2\pi jk/im)
  !>             \quad (k=1,2,...,km) \\
  !>  e_{-k} & = - (1/im)\sum_{j=0}^{im-1} g_j \sin(2\pi jk/im)
  !>              \quad (k=1,2,...,km)
  !> \f}
  function e_g(g)
    !> 格子点データ
    real(DP), dimension(0:im-1), intent(in)  :: g
    !> スペクトルデータ
    real(DP), dimension(-km:km)              :: e_g

    real(DP), dimension(1,size(g))        :: ag_work
    real(DP), dimension(1,-km:km)         :: ae_work

    ag_work(1,:) = g
    ae_work = ae_ag(ag_work)
    e_g = ae_work(1,:)

  end function e_g

  !-----------------------------------------------------------------------
  !> 入力スペクトルデータに X 微分を作用する(2 次元データ).
  !>
  !> スペクトルデータの X 微分とは, 対応する格子点データに X 微分を
  !> 作用させたデータのスペクトル変換のことである.
  function ae_Dx_ae(ae)
    !> 入力スペクトルデータ
    real(DP), dimension(:,-km:), intent(in)  :: ae
    !> 入力スペクトルデータの X 微分
    real(DP), dimension(size(ae,1),-km:km)   :: ae_dx_ae

    integer(8) ::  k

    if ( size(ae,2) < 2*km+1 ) then
       call MessageNotify('W','ae_Dx_ae', &
            'The Fourier dimension of input data too small.')
    elseif ( size(ae,2) > 2*km+1 ) then
       call MessageNotify('W','ae_Dx_ae', &
            'The Fourier dimension of input data too large.')
    endif

    do k=-km,km
       ae_Dx_ae(:,k) = -(2*pi*k/xl)*ae(:,-k)
    enddo
  end function ae_dx_ae

  !-----------------------------------------------------------------------
  !> 入力スペクトルデータに X 微分を作用する(1 次元データ).
  !>
  !> スペクトルデータの X 微分とは, 対応する格子点データに X 微分を
  !> 作用させたデータのスペクトル変換のことである.
  function e_Dx_e(e)
    !> 入力スペクトルデータ
    real(DP), dimension(-km:km), intent(in)     :: e
    !> 入力スペクトルデータの X 微分
    real(DP), dimension(-km:km)                 :: e_Dx_e

    real(DP), dimension(1,-km:km)               :: ae_work

    ae_work(1,:) = e
    ae_work = ae_Dx_ae(ae_work)
    e_Dx_e = ae_work(1,:)

  end function e_Dx_e

  !-----------------------------------------------------------------------
  !> 1 次元格子点データが並んだ 2 次元配列の積分
  function a_Int_ag(ag)
    !> 格子点データ
    real(DP), dimension(:,0:), intent(in)     :: ag
    !> 積分結果
    real(DP), dimension(size(ag,1))           :: a_Int_ag

    integer(8) :: i

    if ( size(ag,2) < im ) then
       call MessageNotify('E','ae_Int_ag', &
            'The Grid points of input data too small.')
    elseif ( size(ag,2) > im ) then
       call MessageNotify('W','ae_Int_ag', &
            'The Grid points of input data too large.')
    endif

    a_Int_ag = 0.0D0
    do i=0,im-1
       a_Int_ag(:) = a_Int_ag(:) + ag(:,i)*g_X_Weight(i)
    enddo
  end function a_Int_ag

  !-----------------------------------------------------------------------
  !> 1 次元格子点データの積分
  function Int_g(g)
    !> 格子点データ
    real(DP), dimension(0:im-1), intent(in)   :: g
    !> 積分結果
    real(DP)                                  :: Int_g

    Int_g = sum(g*g_X_Weight)

  end function Int_g

  !-----------------------------------------------------------------------
  !> 1 次元格子点データが並んだ 2 次元配列の平均
  function a_Avr_ag(ag)
    !> 格子点データ
    real(DP), dimension(:,0:), intent(in)     :: ag
    !> 平均結果
    real(DP), dimension(size(ag,1))           :: a_Avr_ag

    a_Avr_ag = a_Int_ag(ag)/sum(g_X_Weight)

  end function a_Avr_ag

  !-----------------------------------------------------------------------
  !> 1 次元格子点データの平均
  function Avr_g(g)
    !> 格子点データ
    real(DP), dimension(0:im-1), intent(in)   :: g
    !> 平均値
    real(DP)                                  :: Avr_g

    Avr_g = Int_g(g)/sum(g_X_Weight)

  end function Avr_g

  !-----------------------------------------------------------------------
  !> 補間計算(1次元)
  function Interpolate_e(e_data,xval)
    !> 入力スペクトルデータ
    real(DP), dimension(-km:km), intent(IN) :: e_data
    !> 補間する点の座標
    real(DP), intent(in)                    :: xval
    !> 補間した結果の値
    real(DP)                                :: Interpolate_e

    integer(8) :: k

    Interpolate_e = e_data(0)
    do k=1,km
       Interpolate_e = Interpolate_e                             &
            &  + 2.0D0 * (   e_data(k)  * cos(2*PI*k/xl*xval)    &
            &              - e_data(-k) * sin(2*PI*k/xl*xval) )
    end do

  end function Interpolate_e

  !-----------------------------------------------------------------------
  !> 補間計算(2次元)
  function a_Interpolate_ae(ae_data,xval)
    !> 入力スペクトルデータ
    real(DP), dimension(:,-km:), intent(IN)   :: ae_data
    !> 補間する点の座標
    real(DP), intent(in)                      :: xval
    !> 補間した結果の値
    real(DP), dimension(size(ae_data,1))      :: a_Interpolate_ae

    integer(8) :: k

    a_Interpolate_ae = ae_data(:,0)

    do k=1,km
       a_Interpolate_ae = a_Interpolate_ae                          &
            &  + 2.0D0 * (   ae_data(:,k)  * cos(2*PI*k/xl*xval)    &
            &              - ae_data(:,-k) * sin(2*PI*k/xl*xval) )
    end do

  end function a_Interpolate_ae

end module ae_module
