! -*- mode: f90; coding: utf-8 -*-
!--
!------------------------------------------------------------------------
! Copyright (c) 2017 SPMODEL Development Group. All rights reserved.
!------------------------------------------------------------------------
!++
module asc_module
  !
  != asc_module
  !
  ! Copyright :: 2017 Shin-ichi Takehiro, Youhei SASAKI,
  !              SPMODEL Development Group
  ! License :: MIT/X11. Please see COPYRIGHT[link:../COPYRIGHT]
  !
  !== 概要
  !
  ! spml/asc_module モジュールは 1 次元周期境界条件の下での流体運動を実
  ! フーリエ変換によるスペクトル法で数値計算するための Fortran90 関数
  ! を提供する.
  !
  ! 2 次元データの 1 次元に関して同時にスペクトル計算を実行するための
  ! 関数も提供しており, 2, 3 次元領域での計算のベースも提供する.
  !
  ! このモジュールは内部で ISPACK/FTPACK の Fortran77 サブルーチンを呼
  ! んでいる. スペクトルデータの格納方法については ISPACK/FTPACK と異なっ
  ! ているので以下のコメントに注意されたい.
  !
  !== 関数・変数の名前と型について
  !
  !=== 命名法
  !
  ! * 関数名の先頭 (s_, c_, g_, as_, ac_, ag_) は, 返す値の形を示している.
  !   s_  :: スペクトルデータ(sin 変換)
  !   c_  :: スペクトルデータ(cos 変換)
  !   g_  :: 1 次元格子点データ
  !   as_ :: 1 次元スペクトルデータ(sin 変換)が複数並んだ 2 次元データ
  !   ac_ :: 1 次元スペクトルデータ(sin 変換)が複数並んだ 2 次元データ
  !   ag_ :: 1 次元格子点データが複数並んだ 2 次元データ.
  !
  ! * 関数名の間の文字列(Dx)は, その関数の作用を表している.
  !
  ! * 関数名の最後 (_s, _c, _as, _ac, _g, _ag) は, 入力変数の形が
  !   スペクトルデータおよび格子点データであることを示している.
  !   _s  :: スペクトルデータ(sin 変換)
  !   _c  :: スペクトルデータ(cos 変換)
  !   _g  :: 1 次元格子点データ
  !   _as :: 1 次元スペクトルデータ(sin 変換)が複数並んだ 2 次元データ
  !   _ac :: 1 次元スペクトルデータ(cos 変換)が複数並んだ 2 次元データ
  !   _ag :: 1 次元格子点データが複数並んだ 2 次元データ
  !
  !=== 各データの種類の説明
  !
  ! * 以下, DP  = kind(1.0D) である
  ! * g : 1 次元格子点データ.
  !   * 変数の種類と次元は real(DP), dimension(0:im).
  !   * im は X 座標の格子点数であり, サブルーチン ae_Initial にて
  !     あらかじめ設定しておく.
  !
  ! * s : スペクトルデータ(sin 変換).
  !   * 変数の種類と次元は real(8), dimension(km).
  !   * km は X 方向の最大波数であり, サブルーチン asc_Initial にて
  !     あらかじめ設定しておく.
  !
  ! * c : スペクトルデータ(cos 変換).
  !   * 変数の種類と次元は real(8), dimension(0:km).
  !   * km は X 方向の最大波数であり, サブルーチン asc_Initial にて
  !     あらかじめ設定しておく.
  !
  ! * ag : 1 次元(X)格子点データの並んだ 2 次元データ.
  !   * 変数の種類と次元は real(8), dimension(:,0:im).
  !     第 2 次元が X 方向を表す.
  !
  ! * as : 1 次元スペクトルデータ(sin 変換)の並んだ 2 次元データ.
  !   * 変数の種類と次元は real(8), dimension(:,km).
  !     第 2 次元がスペクトルを表す.
  ! * ac : 1 次元スペクトルデータ(cos 変換)の並んだ 2 次元データ.
  !   * 変数の種類と次元は real(8), dimension(:,0:km).
  !     第 2 次元がスペクトルを表す.
  !
  ! * g_ で始まる関数が返す値は 1 次元格子点データに同じ.
  !
  ! * s_ で始まる関数が返す値はスペクトルデータ(sin 変換)に同じ.
  !
  ! * c_ で始まる関数が返す値はスペクトルデータ(cos 変換)に同じ.
  !
  ! * ag_ で始まる関数が返す値は 1 次元格子点データの並んだ
  !   2 次元データに同じ.
  !
  ! * as_ で始まる関数が返す値は 1 次元スペクトルデータ(sin 変換)の並んだ
  !   2 次元データに同じ.
  !
  ! * ac_ で始まる関数が返す値は 1 次元スペクトルデータ(cos 変換)の並んだ
  !   2 次元データに同じ.
  !
  ! * スペクトルデータに対する微分等の作用とは, 対応する格子点データに
  !   微分などを作用させたデータをスペクトル変換したもののことである.
  !
  !== 変数・手続き群の要約
  !
  !==== 初期化
  !
  ! asc_Initial       :: スペクトル変換の格子点数, 波数, 領域の大きさの設定
  !
  !==== 座標変数
  !
  ! g_X              :: 格子点座標(X)を格納した 1 次元配列
  ! g_X_Weight       :: 重み座標を格納した 1 次元配列
  !
  !==== 基本変換
  !
  ! g_s, ag_as       :: スペクトルデータ(sin 変換)から格子データへの変換
  ! g_c, ag_ac       :: スペクトルデータ(cos 変換)から格子データへの変換
  ! s_g, as_ag       :: 格子データからスペクトルデータ(sin 変換)への変換
  ! c_g, ac_ag       :: 格子データからスペクトルデータ(sin 変換)への変換
  !
  !==== 微分
  !
  ! c_Dx_s, ac_Dx_as :: スペクトルデータ(sin 変換)に X 微分を作用させる
  ! s_Dx_c, as_Dx_ac :: スペクトルデータ(cos 変換)に X 微分を作用させる
  !
  !==== 積分・平均
  !
  ! a_Int_ag, a_Avr_ag :: 1 次元格子点データの並んだ 2 次元配列の積分および平均
  ! Int_g, Avr_g       :: 1 次元格子点データの積分および平均
  !
  !==== 補間
  !
  ! Interpolate_s      :: スペクトルデータからの特定の座標の値の補間(1 次元)
  ! Interpolate_c      :: スペクトルデータからの特定の座標の値の補間(1 次元)
  ! a_Interpolate_as   :: スペクトルデータからの特定の座標の値の補間(2 次元)
  ! a_Interpolate_ac   :: スペクトルデータからの特定の座標の値の補間(2 次元)
  !
  !
  use dc_message
  use dc_types, only: DP
  implicit none
  private
  public asc_Initial                      ! 初期化ルーチン
  public ag_as, as_ag, g_s, s_g           ! 基本変換
  public ag_ac, ac_ag, g_c, c_g           ! 基本変換
  public ac_Dx_as, c_Dx_s                 ! 微分
  public as_Dx_ac, s_Dx_c                 ! 微分
  public a_Int_ag, Int_g, a_Avr_ag, Avr_g ! 積分・平均
  public Interpolate_s, a_Interpolate_as  ! 補間
  public Interpolate_c, a_Interpolate_ac  ! 補間
  public g_X, g_X_Weight                  ! 座標変数
  integer            :: im=32             ! 格子点の数
  integer            :: km=10             ! 切断波数
  real(DP)           :: xl=1.0            ! 領域の大きさ
  integer, dimension(5)               :: its, itc
  real(DP), dimension(:),allocatable  :: ts, tc
  real(DP), parameter                 :: pi=3.1415926535897932385D0
  real(DP), allocatable :: g_x(:)         ! 格子点座標(X)を格納した 1 次元配列.
  real(DP), allocatable :: g_x_weight(:)  ! 重み座標を格納した 1 次元配列.
                                          ! X 方向の格子点の間隔が格納してある.
  save im, km, its, ts, itc, tc, xl, g_X, g_X_Weight
  contains
    !--------------- 初期化 -----------------
    subroutine asc_Initial(i,k,xmin,xmax)
      !
      ! スペクトル変換の格子点数, 波数, 領域の大きさを設定する.
      !
      ! 他の関数や変数を呼ぶ前に, 最初にこのサブルーチンを呼んで
      ! 初期設定をしなければならない.
      !
      integer,intent(in)  :: i              ! 格子点の数
      integer,intent(in)  :: k              ! 切断波数
      real(DP),intent(in) :: xmin, xmax     ! X 座標範囲
      integer :: ii
      im = i
      km = k
      xl = xmax-xmin
      if ( im <= 0 .or. km <= 0 ) then
         call MessageNotify('E','asc_Initial', &
              'Number of grid points and waves should be positive')
      elseif ( mod(im,2) /= 0 ) then
         call MessageNotify('E','asc_Initial', &
              'Number of grid points should be even')
      endif
      allocate(ts(5*im/2),tc(im*3))
      call fttcti(im,itc,tc)
      call fttsti(im,its,ts)
      allocate(g_x(0:im))
      do ii=0,im
         g_X(ii) = xmin + xl/im*ii
      enddo
      allocate(g_x_weight(0:im))
      g_X_Weight = xl/im
      g_X_Weight(0)  = xl/(2*im)
      g_X_Weight(im) = xl/(2*im)
      call MessageNotify('M','asc_initial','asc_module for OpenMP (2017/05/18) is initialized')
    end subroutine asc_Initial
    !--------------- 基本変換(逆変換) -----------------
    function ag_as(as)
      !
      ! スペクトルデータ(sin 変換)から格子点データへ逆変換する(2 次元データ用)
      !
      ! sin スペクトル逆変換の定義は以下のとおり.
      !
      !   スペクトルデータ s_k (k=1,...,km)
      !   格子点データ     g_j (j=0,...,im)
      !
      !  g_j = \sum_{k=1}^{km}(s_k\sin(\pi jk/im)   (j=0,1,...,im).
      !
      real(DP), dimension(:,:), intent(in)   :: as    !(in)  スペクトルデータ
      real(DP), dimension(size(as,1),0:im)   :: ag_as !(out) 格子点データ
      real(DP), dimension(size(as,1),im)   :: y
      integer :: n, nm, k
      nm=size(as,1)
      ag_as = 0.0D0
      do k=1,km
         ag_as(:,k)   = as(:,k)
      enddo
!$omp parallel do default(shared)
      do n=1,nm
        call fttstb(1,im,ag_as(n,1:im),y(n,:),its,ts)
      enddo
!$omp end parallel do
    end function ag_as
    function g_s(s)
      !
      ! スペクトルデータ(sin 変換)から格子点データへ逆変換する(1 次元データ用)
      !
      ! sin スペクトル逆変換の定義は以下のとおり.
      !
      !   スペクトルデータ s_k (k=1,...,km)
      !   格子点データ     g_j (j=0,...,im)
      !
      !  g_j = \sum_{k=1}^{km}(s_k\sin(\pi jk/im)   (j=0,1,...,im).
      !
      real(DP), dimension(0:im)             :: g_s  !(out) 格子点データ
      real(DP), dimension(km), intent(in)   :: s    !(in)  スペクトルデータ
      real(DP), dimension(1,size(s))  :: as_work
      real(DP), dimension(1,0:im)     :: ag_work
      as_work(1,:) = s
      ag_work = ag_as(as_work)
      g_s = ag_work(1,:)
    end function g_s
    function ag_ac(ac)
      !
      ! スペクトルデータ(cos 変換)から格子点データへ逆変換する(2 次元データ用)
      !
      ! cos スペクトル逆変換の定義は以下のとおり.
      !
      !   スペクトルデータ c_k (k=0,...,km)
      !   格子点データ     g_j (j=0,...,im)
      !
      !  g_j = c_0/2 + c_im(-1)**j
      !      + \sum_{k=1}^{im-1}(c_k\cos(\pi jk/im)   (j=0,1,...,im).
      !
      real(DP), dimension(:,0:), intent(in)   :: ac    !(in)  スペクトルデータ
      real(DP), dimension(size(ac,1),0:im)    :: ag_ac !(out) 格子点データ
      real(DP), dimension(size(ac,1),im)      :: y
      integer :: n, nm, k
      nm=size(ac,1)
      ag_ac = 0.0D0
      do k=0,km
         ag_ac(:,k)   = ac(:,k)
      enddo
!$omp parallel do default(shared)
      do n=1,nm
        call fttctb(1,im,ag_ac(n,:),y(n,:),itc,tc)
      enddo
!$omp end parallel do
    end function ag_ac
    function g_c(c)
      !
      ! スペクトルデータ(cos 変換)から格子点データへ逆変換する(1 次元データ用)
      !
      ! cos スペクトル逆変換の定義は以下のとおり.
      !
      !   スペクトルデータ c_k (k=0,...,km)
      !   格子点データ     g_j (j=0,...,im)
      !
      !  g_j = c_0/2 + c_im(-1)**j
      !      + \sum_{k=1}^{im-1}(c_k\cos(\pi jk/im)   (j=0,1,...,im).
      !
      real(DP), dimension(0:im)             :: g_c  !(out) 格子点データ
      real(DP), dimension(0:km), intent(in) :: c    !(in)  スペクトルデータ
      real(DP), dimension(1,size(c))  :: ac_work
      real(DP), dimension(1,0:im)     :: ag_work
      ac_work(1,:) = c
      ag_work = ag_ac(ac_work)
      g_c = ag_work(1,:)
    end function g_c
    !--------------- 基本変換(正変換) -----------------
    function as_ag(ag)
      !
      ! 格子点データからスペクトルデータ(sin 変換)へ正変換する(2 次元データ用)
      !
      ! スペクトル正変換の定義は以下のとおり.
      !
      !   格子点データ     g_j (j=0,...,im)
      !   スペクトルデータ s_k (k=1,...,km)
      !
      !  s_k = (2/im)\sum_{j=1}^{im-1} g_j \sin(\pi jk/im) ]    (k=1,...,km)
      !
      real(DP), dimension(:,:), intent(in)   :: ag     !(in)  格子点データ
      real(DP), dimension(size(ag,1),km)     :: as_ag  !(out) スペクトルデータ
      real(DP), dimension(size(ag,1),im)     :: y
      real(DP), dimension(size(ag,1),0:im)   :: ag_work
      integer :: n, nm, k
      nm = size(ag,1)
      ag_work = ag
!$omp parallel do default(shared)
      do n=1,nm
        call fttstf(1,im,ag_work(n,1:im),y(n,:),its,ts)
      enddo
!$omp end parallel do
      do k=1,km
         as_ag(:,k) = ag_work(:,k)
      enddo
    end function as_ag
    function s_g(g)
      !
      ! 格子点データからスペクトルデータ(sin 変換)へ正変換する(2 次元データ用)
      !
      ! スペクトル正変換の定義は以下のとおり.
      !
      !   格子点データ     g_j (j=0,...,im)
      !   スペクトルデータ s_k (k=1,...,km)
      !
      !  s_k = (2/im)\sum_{j=1}^{im-1} g_j \sin(\pi jk/im) ]    (k=1,...,km)
      !
      real(DP), dimension(km)                  :: s_g   !(out) スペクトルデータ
      real(DP), dimension(0:im), intent(in)    :: g     !(in)  格子点データ
      real(DP), dimension(1,size(g))           :: ag_work
      real(DP), dimension(1,km)                :: as_work
      ag_work(1,:) = g
      as_work = as_ag(ag_work)
      s_g = as_work(1,:)
    end function s_g
    function ac_ag(ag)
      !
      ! 格子点データからスペクトルデータ(cos 変換)へ正変換する(2 次元データ用)
      !
      ! スペクトル正変換の定義は以下のとおり.
      !
      !   格子点データ     g_j (j=0,...,im)
      !   スペクトルデータ c_k (k=0,...,km)
      !
      !  cs_k = (2/im)[ (1/2)g_0 + (1/2)g_im*(-1)^k
      !       + \sum_{j=1}^{im-1} g_j \cos(\pi jk/im) ]    (k=0,1,...,km)
      !
      real(DP), dimension(:,:), intent(in)   :: ag     !(in)  格子点データ
      real(DP), dimension(size(ag,1),0:km)   :: ac_ag  !(out) スペクトルデータ
      real(DP), dimension(size(ag,1),im)     :: y
      real(DP), dimension(size(ag,1),0:im)   :: ag_work
      integer :: n, nm, k
      nm = size(ag,1)
      ag_work = ag
!$omp parallel do default(shared)
      do n=1,nm
        call fttctf(1,im,ag_work(n,:),y(n,:),itc,tc)
      enddo
!$omp end parallel do
      do k=0,km
         ac_ag(:,k) = ag_work(:,k)
      enddo
    end function ac_ag
    function c_g(g)
      !
      ! 格子点データからスペクトルデータ(cos 変換)へ正変換する(1 次元データ用)
      !
      ! スペクトル正変換の定義は以下のとおり.
      !
      !   格子点データ     g_j (j=0,...,im)
      !   スペクトルデータ c_k (k=0,...,km)
      !
      !  cs_k = (2/im)[ (1/2)g_0 + (1/2)g_im*(-1)^k
      !       + \sum_{j=1}^{im-1} g_j \cos(\pi jk/im) ]    (k=0,1,...,km)
      !
      real(DP), dimension(0:km)                :: c_g   !(out) スペクトルデータ
      real(DP), dimension(0:im), intent(in)    :: g     !(in)  格子点データ
      real(DP), dimension(1,size(g))           :: ag_work
      real(DP), dimension(1,0:km)              :: ac_work
      ag_work(1,:) = g
      ac_work = ac_ag(ag_work)
      c_g = ac_work(1,:)
    end function c_g
  !--------------- 微分計算 -----------------
    function ac_Dx_as(as)
      !
      ! 入力スペクトルデータに X 微分を作用する(2 次元データ).
      !
      ! スペクトルデータの X 微分とは, 対応する格子点データに X 微分を
      ! 作用させたデータのスペクトル変換のことである.
      !
      !
      real(DP), dimension(:,:), intent(in)  :: as
      !(in)  入力スペクトルデータ
      real(DP), dimension(size(as,1),0:km)  :: ac_dx_as
      !(out) 入力スペクトルデータの X 微分
      integer ::  k
      ac_Dx_as(:,0) = 0.0D0
      do k=1,km
         ac_Dx_as(:,k) = (pi*k/xl)*as(:,k)
      enddo
    end function ac_Dx_as
    function c_Dx_s(s)
      !
      ! 入力スペクトルデータに X 微分を作用する(1 次元データ).
      !
      ! スペクトルデータの X 微分とは, 対応する格子点データに X 微分を
      ! 作用させたデータのスペクトル変換のことである.
      !
      !
      real(DP), dimension(km), intent(in)     :: s
      !(in)  入力スペクトルデータ
      real(DP), dimension(0:km)               :: c_Dx_s
      !(out) 入力スペクトルデータの X 微分
      real(DP), dimension(1,km)               :: as_work
      real(DP), dimension(1,0:km)             :: ac_work
      as_work(1,:) = s
      ac_work = ac_Dx_as(as_work)
      c_Dx_s= ac_work(1,:)
    end function c_Dx_s
    function as_Dx_ac(ac)
      !
      ! 入力スペクトルデータに X 微分を作用する(2 次元データ).
      !
      ! スペクトルデータの X 微分とは, 対応する格子点データに X 微分を
      ! 作用させたデータのスペクトル変換のことである.
      !
      !
      real(DP), dimension(:,0:), intent(in) :: ac
      !(in)  入力スペクトルデータ
      real(DP), dimension(size(ac,1),km)    :: as_dx_ac
      !(out) 入力スペクトルデータの X 微分
      integer ::  k
      do k=1,km
         as_Dx_ac(:,k) = -(pi*k/xl)*ac(:,k)
      enddo
    end function as_Dx_ac
    function s_Dx_c(c)
      !
      ! 入力スペクトルデータに X 微分を作用する(1 次元データ).
      !
      ! スペクトルデータの X 微分とは, 対応する格子点データに X 微分を
      ! 作用させたデータのスペクトル変換のことである.
      !
      !
      real(DP), dimension(0:km), intent(in)   :: c
      !(in)  入力スペクトルデータ
      real(DP), dimension(km)                 :: s_Dx_c
      !(out) 入力スペクトルデータの X 微分
      real(DP), dimension(1,km)               :: as_work
      real(DP), dimension(1,0:km)             :: ac_work
      ac_work(1,:) = c
      as_work = as_Dx_ac(ac_work)
      s_Dx_c = as_work(1,:)
    end function s_Dx_c
  !--------------- 積分計算 -----------------
    function a_Int_ag(ag)
      !
      ! 1 次元格子点データが並んだ 2 次元配列の積分
      !
      real(DP), dimension(:,0:), intent(in)     :: ag        !(in) 格子点データ
      real(DP), dimension(size(ag,1))           :: a_Int_ag  !(out) 積分結果
      integer :: i
      a_Int_ag = 0.0D0
      do i=0,im
         a_Int_ag(:) = a_Int_ag(:) + ag(:,i)*g_X_Weight(i)
      enddo
    end function a_Int_ag
    function Int_g(g)
      !
      ! 1 次元格子点データの積分
      !
      real(DP), dimension(0:im), intent(in)   :: g      !(in) 格子点データ
      real(DP)                                :: Int_g  !(out) 積分結果
      Int_g = sum(g*g_X_Weight)
    end function Int_g
    function a_Avr_ag(ag)
      !
      ! 1 次元格子点データが並んだ 2 次元配列の平均
      !
      real(DP), dimension(:,0:), intent(in)     :: ag        !(in) 格子点データ
      real(DP), dimension(size(ag,1))           :: a_Avr_ag  !(out) 平均結果
      a_Avr_ag = a_Int_ag(ag)/sum(g_X_Weight)
    end function a_Avr_ag
    function Avr_g(g)
      !
      ! 1 次元格子点データの平均
      !
      real(DP), dimension(0:im), intent(in)   :: g
      real(DP)                                :: Avr_g
      Avr_g = Int_g(g)/sum(g_X_Weight)
    end function Avr_g
  !--------------- 補間計算 -----------------
    function Interpolate_s(s_data,xval)
      real(DP), dimension(km), intent(IN) :: s_data
      !(in) 入力フーリエデータ
      real(DP), intent(in)                :: xval
      ! 補間する点の座標
      real(DP)                            :: Interpolate_s
      ! 補間した結果の値
      integer :: k
      ! DO 文変数
      Interpolate_s = 0.0D0
      do k=1,km
        Interpolate_s = Interpolate_s + s_data(k)*sin(PI*k/xl*xval-g_X(0))
      enddo
    end function Interpolate_s
    function a_Interpolate_as(as_data,xval)
      real(DP), dimension(:,:), intent(IN)   :: as_data
      !(in) 入力フーリエデータ
      real(DP), intent(in)                   :: xval
      ! 補間する点の座標
      real(DP), dimension(size(as_data,1))   :: a_Interpolate_as
      ! 補間した結果の値
      integer :: k
      ! DO 文変数
      a_Interpolate_as = 0.0D0
      do k=1,km
         a_Interpolate_as = a_Interpolate_as + as_data(:,k)*sin(PI*k/xl*(xval-g_X(0)))
      enddo
    end function a_Interpolate_as
    function Interpolate_c(c_data,xval)
      !
      real(DP), dimension(0:km), intent(IN) :: c_data
      !(in) 入力フーリエデータ
      real(DP), intent(in)                  :: xval
      ! 補間する点の座標
      real(DP)                              :: Interpolate_c
      ! 補間した結果の値
      integer :: k
      ! DO 文変数
      Interpolate_c = c_data(0)/2
      do k=1,km
         Interpolate_c = Interpolate_c + c_data(k)*cos(PI*k/xl*(xval-g_X(0)))
      enddo
    end function Interpolate_c
    function a_Interpolate_ac(ac_data,xval)
      real(DP), dimension(:,0:), intent(IN)  :: ac_data
      !(in) 入力フーリエデータ
      real(DP), intent(in)                   :: xval
      ! 補間する点の座標
      real(DP), dimension(size(ac_data,1))   :: a_Interpolate_ac
      ! 補間した結果の値
      integer :: k
      ! DO 文変数
      a_Interpolate_ac = ac_data(:,0)/2
      do k=1,km
         a_Interpolate_ac = a_Interpolate_ac + ac_data(:,k)*cos(PI*k/xl*(xval-g_X(0)))
      enddo
    end function a_Interpolate_ac
end module asc_module