module fft_work
  use dcl_common

  type work
    integer                     :: n
    real, dimension(:), pointer :: array
  end type work

end module
!-------------------------------------------------
!  FFT  周期実数値データのフーリエ変換 (1)
!-------------------------------------------------
module fftreal
  use fft_work
  type(work), dimension(10), private :: wk

  contains
!---------------------------------------------------------
  subroutine DclInitRealFFT(n, index)       !初期化をおこなう
    integer, intent(in)           :: n      !処理するデータの長さ
    integer, intent(in), optional :: index  !作業領域番号

    if(present(index)) then; idx = index
                       else; idx = 1     ; end if

    if(associated(wk(idx)%array)) call msgdmp('E', 'DclInitRealFFT', &
                   & 'Working area has been allocated already.')
    wk(idx)%n = n
    len = 2*n+15
    allocate(wk(idx)%array(len))

    call rffti(n, wk(idx)%array) 
  end subroutine
!---------------------------------------------------------
  subroutine DclDeallocRealFFT(index)
    integer, intent(in), optional :: index  !作業領域番号

    if(present(index)) then; idx = index
                       else; idx = 1     ; end if
    deallocate(wk(idx)%array)
  end subroutine
!---------------------------------------------------------
  function DclRealFFT_F(r, index)
    real,    intent(in), dimension(:) :: r
    integer, intent(in), optional     :: index
    real,    dimension(size(r))       :: DclRealFFT_F

    if(present(index)) then; idx = index
                       else; idx = 1     ; end if

    n = size(r)
    if(n.ne.wk(idx)%n) call msgdmp('E', 'DclRealFFT_F', &
                   & 'Wrong working area.')
    DclRealFFT_F = r
    call rfftf(n, DclRealFFT_F, wk(idx)%array) 
  end function
!---------------------------------------------------------
  function DclRealFFT_B(r, index)
    real,    intent(in), dimension(:) :: r
    integer, intent(in), optional     :: index
    real,    dimension(size(r))       :: DclRealFFT_B

    if(present(index)) then; idx = index
                       else; idx = 1     ; end if

    n = size(r)
    if(n.ne.wk(idx)%n) call msgdmp('E', 'DclRealFFT_B', &
                   & 'Wrong working area.')
    DclRealFFT_B = r
    call rfftb(n, DclRealFFT_B, wk(idx)%array) 
  end function
end module
!-------------------------------------------------
! FFT: rffti, rfftf, rfftbの簡易型サブルーチン (2)
!-------------------------------------------------
module ffteasy
  use fft_work
  type(work), dimension(10), private :: wk

  contains
!---------------------------------------------------------
  subroutine DclInitEasyFFT(n, index)       !初期化をおこなう
    integer, intent(in)           :: n      !処理するデータの長さ
    integer, intent(in), optional :: index  !作業領域番号

    if(present(index)) then; idx = index
                       else; idx = 1     ; end if

    if(associated(wk(idx)%array)) call msgdmp('E', 'DclInitEasyFFT', &
                   & 'Working area has been allocated already.')
    wk(idx)%n = n
    len = 3*n+15
    allocate(wk(idx)%array(len))

    call ezffti(n, wk(idx)%array) 
  end subroutine
!---------------------------------------------------------
  subroutine DclDeallocEasyFFT(index)
    integer, intent(in), optional :: index  !作業領域番号

    if(present(index)) then; idx = index
                       else; idx = 1     ; end if
    deallocate(wk(idx)%array)
  end subroutine
!---------------------------------------------------------
  subroutine DclEasyFFT_F(r,a0,a,b,index)
    real,    intent(in),  dimension(:) :: r
    real,    intent(out)               :: a0 
    real,    intent(out), dimension(:) :: a
    real,    intent(out), dimension(:) :: b
    integer, intent(in),  optional     :: index

    if(present(index)) then; idx = index
                       else; idx = 1     ; end if

    n = size(r)
    if(n.ne.wk(idx)%n) call msgdmp('E', 'DclEasyFFT_F', &
                   & 'Wrong working area.')
   call ezfftf(n,r,a0,a,b,wk(idx)%array) 
  end subroutine
!---------------------------------------------------------
  subroutine DclEasyFFT_B(r,a0,a,b,index)
    real,    intent(out), dimension(:) :: r
    real,    intent(in)                :: a0
    real,    intent(in),  dimension(:) :: a
    real,    intent(in),  dimension(:) :: b 
    integer, intent(in),  optional     :: index

    if(present(index)) then; idx = index
                       else; idx = 1     ; end if

    n = size(r)
    if(n.ne.wk(idx)%n) call msgdmp('E', 'DclEasyFFT_B', &
                   & 'Wrong working area.')
    call ezfftb(n,r,a0,a,b,wk(idx)%array) 
  end subroutine
end module
!-------------------------------------------------
! 奇の周期データのsine変換をおこなう． (3)
!-------------------------------------------------
module fftsin
  use fft_work
  type(work), dimension(10), private :: wk

  contains
!---------------------------------------------------------
  subroutine DclInitSinFFT(n,index)
    integer, intent(in)           :: n      !処理するデータの長さ
    integer, intent(in), optional :: index  !作業領域番号

    if(present(index)) then; idx = index
                       else; idx = 1     ; end if

    if(associated(wk(idx)%array)) call msgdmp('E', 'DclInitSinFFT', &
                   & 'Working area has been allocated already.')
    wk(idx)%n = n
    len = 2.5*n+15
    allocate(wk(idx)%array(len))

    call sinti(n, wk(idx)%array)
  end subroutine
!---------------------------------------------------------
  subroutine DclDeallocSinFFT(index)
    integer, intent(in), optional :: index  !作業領域番号

    if(present(index)) then; idx = index
                       else; idx = 1     ; end if
    deallocate(wk(idx)%array)
  end subroutine
!---------------------------------------------------------
  function DclSinFFT(r,index)
    real,    intent(in), dimension(:) :: r
    integer, intent(in), optional     :: index
    real,    dimension(size(r))       :: DclSinFFT

    if(present(index)) then; idx = index
                       else; idx = 1     ; end if

    n = size(r)
    if(n.ne.wk(idx)%n) call msgdmp('E', 'DclSinFFT', &
                   & 'Wrong working area.')
    DclSinFFT = r
    call sint(n, DclSinFFT, wk(idx)%array)
  end function
end module
!-------------------------------------------------
! FFT: 偶の周期データのcosine変換をおこなう (4)
!-------------------------------------------------
module fftcos
  use fft_work
  type(work), dimension(10), private :: wk

  contains
!---------------------------------------------------------
  subroutine DclInitCosFFT(n,index) 
    integer, intent(in)           :: n      !処理するデータの長さ
    integer, intent(in), optional :: index  !作業領域番号

    if(present(index)) then; idx = index
                       else; idx = 1     ; end if

    if(associated(wk(idx)%array)) call msgdmp('E', 'DclInitCosFFT', &
                   & 'Working area has been allocated already.')
    wk(idx)%n = n
    len = 3*n+15
    allocate(wk(idx)%array(len))

    call costi(n, wk(idx)%array)
  end subroutine
!---------------------------------------------------------
  subroutine DclDeallocCosFFT(index)
    integer, intent(in), optional :: index  !作業領域番号

    if(present(index)) then; idx = index
                       else; idx = 1     ; end if
    deallocate(wk(idx)%array)
  end subroutine
!---------------------------------------------------------
  function DclCosFFT(r, index)
    real,    intent(in), dimension(:) :: r
    integer, intent(in), optional     :: index
    real,    dimension(size(r))       :: DclCosFFT

    if(present(index)) then; idx = index
                       else; idx = 1     ; end if

    n = size(r)
    if(n.ne.wk(idx)%n) call msgdmp('E', 'DclCosFFT', &
                   & 'Wrong working area.')
    DclCosFFT = r
    call cost(n,DclCosFFT,wk(idx)%array)
  end function
end module
!-------------------------------------------------
! FFT: 奇数波数成分のみのsin変換をおこなう．(5)
!-------------------------------------------------
module fftqsin
  use fft_work
  type(work), dimension(10), private :: wk

  contains
!---------------------------------------------------------
  subroutine DclInitSinQFT(n,index)
    integer, intent(in)           :: n      !処理するデータの長さ
    integer, intent(in), optional :: index  !作業領域番号

    if(present(index)) then; idx = index
                       else; idx = 1     ; end if

    if(associated(wk(idx)%array)) call msgdmp('E', 'DclInitSinQFT', &
                   & 'Working area has been allocated already.')
    wk(idx)%n = n
    len = 3*n+15
    allocate(wk(idx)%array(len))

    call sinqi(n, wk(idx)%array)
  end subroutine
!---------------------------------------------------------
  subroutine DclDeallocSinQFT(index)
    integer, intent(in), optional :: index  !作業領域番号

    if(present(index)) then; idx = index
                       else; idx = 1     ; end if
    deallocate(wk(idx)%array)
  end subroutine
!---------------------------------------------------------
  function DclSinQFT_F(r, index)
    real,    intent(in), dimension(:) :: r
    integer, intent(in), optional     :: index
    real,    dimension(size(r))       :: DclSinQFT_F

    if(present(index)) then; idx = index
                       else; idx = 1     ; end if

    n = size(r)
    if(n.ne.wk(idx)%n) call msgdmp('E', 'DclSinQFT_F', &
                   & 'Wrong working area.')
    DclSinQFT_F = r
    call sinqf(n, DclSinQFT_F, wk(idx)%array)
  end function
!---------------------------------------------------------
  function DclSinQFT_B(r, index)
    real,    intent(in), dimension(:) :: r
    integer, intent(in), optional     :: index
    real,    dimension(size(r))       :: DclSinQFT_B

    if(present(index)) then; idx = index
                       else; idx = 1     ; end if

    n = size(r)
    if(n.ne.wk(idx)%n) call msgdmp('E', 'DclSinQFT_B', &
                   & 'Wrong working area.')
    DclSinQFT_B = r
    call sinqb(n, DclSinQFT_B, wk(idx)%array)
  end function
end module
!-------------------------------------------------
! FFT:偶数波数成分のみのcossine変換をおこなう．(6)
!-------------------------------------------------
module fftqcos
  use fft_work
  type(work), dimension(10), private :: wk

  contains
!---------------------------------------------------------
  subroutine DclInitCosQFT(n,index)
    integer, intent(in)           :: n      !処理するデータの長さ
    integer, intent(in), optional :: index  !作業領域番号

    if(present(index)) then; idx = index
                       else; idx = 1     ; end if

    if(associated(wk(idx)%array)) call msgdmp('E', 'DclInitCosQFT', &
                   & 'Working area has been allocated already.')
    wk(idx)%n = n
    len = 3*n+15
    allocate(wk(idx)%array(len))

    call cosqi(n, wk(idx)%array)
  end subroutine
!---------------------------------------------------------
  subroutine DclDeallocCosQFT(index)
    integer, intent(in), optional :: index  !作業領域番号

    if(present(index)) then; idx = index
                       else; idx = 1     ; end if
    deallocate(wk(idx)%array)
  end subroutine
!---------------------------------------------------------
  function DclCosQFT_F(r, index)
    real,    intent(in), dimension(:) :: r
    integer, intent(in), optional     :: index
    real,    dimension(size(r))       :: DclCosQFT_F

    if(present(index)) then; idx = index
                       else; idx = 1     ; end if

    n = size(r)
    if(n.ne.wk(idx)%n) call msgdmp('E', 'DclCosQFT_F', &
                   & 'Wrong working area.')
    DclCosQFT_F = r
    call cosqf(n, DclCosQFT_F, wk(idx)%array)
  end function
!---------------------------------------------------------
  function DclCosQFT_B(r, index)
    real,    intent(in), dimension(:) :: r
    integer, intent(in), optional     :: index
    real,    dimension(size(r))       :: DclCosQFT_B

    if(present(index)) then; idx = index
                       else; idx = 1     ; end if

    n = size(r)
    if(n.ne.wk(idx)%n) call msgdmp('E', 'DclCosQFT_B', &
                   & 'Wrong working area.')
    DclCosQFT_B = r
    call cosqb(size(r), DclCosQFT_B, wk(idx)%array)
  end function
end module
!-------------------------------------------------
! 周期複素数データのフーリエ変換をおこなう  (7)
!-------------------------------------------------
module fftcmplx
  use fft_work
  type(work), dimension(10), private :: wk

  contains
!---------------------------------------------------------
  subroutine DclInitComplexFFT(n,index)
    integer, intent(in)           :: n      !処理するデータの長さ
    integer, intent(in), optional :: index  !作業領域番号

    if(present(index)) then; idx = index
                       else; idx = 1     ; end if

    if(associated(wk(idx)%array)) call msgdmp('E', 'DclInitComplexFFT', &
                   & 'Working area has been allocated already.')
    wk(idx)%n = n
    len = 4*n+15
    allocate(wk(idx)%array(len))

    call cffti(n, wk(idx)%array)
  end subroutine
!---------------------------------------------------------
  subroutine DclDeallocComplexFFT(index)
    integer, intent(in), optional :: index  !作業領域番号

    if(present(index)) then; idx = index
                       else; idx = 1     ; end if
    deallocate(wk(idx)%array)
  end subroutine
!---------------------------------------------------------
  function DclComplexFFT_F(c,index)
    complex, intent(in), dimension(:) :: c
    integer, intent(in), optional     :: index
    complex, dimension(size(c))       :: DclComplexFFT_F

    if(present(index)) then; idx = index
                       else; idx = 1     ; end if

    n = size(c)
    if(n.ne.wk(idx)%n) call msgdmp('E', 'DclComplexFFT_F', &
                   & 'Wrong working area.')
    DclComplexFFT_F = c
    call cfftf(n, DclComplexFFT_F, wk(idx)%array)
  end function
!---------------------------------------------------------
  function DclComplexFFT_B(c,index)
    complex, intent(in), dimension(:) :: c
    integer, intent(in), optional     :: index
    complex, dimension(size(c))       :: DclComplexFFT_B

    if(present(index)) then; idx = index
                       else; idx = 1     ; end if

    n = size(c)
    if(n.ne.wk(idx)%n) call msgdmp('E', 'DclComplexFFT_B', &
                   & 'Wrong working area.')
    DclComplexFFT_B = c
    call cfftb(n, DclComplexFFT_B, wk(idx)%array)
  end function
end module
!-------------------------------------------------
!  fftlib Module
!-------------------------------------------------
module fftlib
  use fftreal
  use ffteasy
  use fftsin
  use fftcos
  use fftqsin
  use fftqcos
  use fftcmplx
  private :: work
end module
