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  !作業領域番号

    call prcopn('DclInitRealFFT')
    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) 
    call prccls('DclInitRealFFT')
  end subroutine
!---------------------------------------------------------
  subroutine DclDeallocRealFFT(index)
    integer, intent(in), optional :: index  !作業領域番号

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

    call prcopn('DclRealFFT_F')
    if(present(index)) then; idx = index
                       else; idx = 1     ; end if

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

    call prcopn('DclRealFFT_B')
    if(present(index)) then; idx = index
                       else; idx = 1     ; end if

    n = size(x)
    if(n.ne.wk(idx)%n) call msgdmp('E', 'DclrealFFT_B', &
                   & 'Wrong working area.')
    DclRealFFT_B = x
    call rfftb(n, DclRealFFT_B, wk(idx)%array) 
    call prccls('DclRealFFT_B')
  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  !作業領域番号

    call prcopn('DclInitEasyFFT')
    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) 
    call prccls('DclInitEasyFFT')
  end subroutine
!---------------------------------------------------------
  subroutine DclDeallocEasyFFT(index)
    integer, intent(in), optional :: index  !作業領域番号

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

    call prcopn('DclEasyFFT_F')
    if(present(index)) then; idx = index
                       else; idx = 1     ; end if

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

    call prcopn('DclEasyFFT_B')
    if(present(index)) then; idx = index
                       else; idx = 1     ; end if

    n = size(x)
    if(n.ne.wk(idx)%n) call msgdmp('E', 'DclEasyFFT_B', &
                   & 'Wrong working area.')
    call ezfftb(n,x,a0,a,b,wk(idx)%array) 
    call prccls('DclEasyFFT_B')
  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  !作業領域番号

    call prcopn('DclInitSinFFT')
    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)
    call prccls('DclInitSinFFT')
  end subroutine
!---------------------------------------------------------
  subroutine DclDeallocSinFFT(index)
    integer, intent(in), optional :: index  !作業領域番号

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

    call prcopn('DclSinFFT')
    if(present(index)) then; idx = index
                       else; idx = 1     ; end if

    n = size(x)
    if(n.ne.wk(idx)%n) call msgdmp('E', 'DclSinFFT', &
                   & 'Wrong working area.')
    DclSinFFT = x
    call sint(n, DclSinFFT, wk(idx)%array)
    call prccls('DclSinFFT')
  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  !作業領域番号

    call prcopn('DclInitCosFFT')
    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)
    call prccls('DclInitCosFFT')
  end subroutine
!---------------------------------------------------------
  subroutine DclDeallocCosFFT(index)
    integer, intent(in), optional :: index  !作業領域番号

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

    call prcopn('DclCosFFT')
    if(present(index)) then; idx = index
                       else; idx = 1     ; end if

    n = size(x)
    if(n.ne.wk(idx)%n) call msgdmp('E', 'DclCosFFT', &
                   & 'Wrong working area.')
    DclCosFFT = x
    call cost(n,DclCosFFT,wk(idx)%array)
    call prccls('DclCosFFT')
  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  !作業領域番号

    call prcopn('DclInitSinQFT')
    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)
    call prccls('DclInitSinQFT')
  end subroutine
!---------------------------------------------------------
  subroutine DclDeallocSinQFT(index)
    integer, intent(in), optional :: index  !作業領域番号

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

    call prcopn('DclSinQFT_F')
    if(present(index)) then; idx = index
                       else; idx = 1     ; end if

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

    call prcopn('DclSinQFT_B')
    if(present(index)) then; idx = index
                       else; idx = 1     ; end if

    n = size(x)
    if(n.ne.wk(idx)%n) call msgdmp('E', 'DclSinQFT_B', &
                   & 'Wrong working area.')
    DclSinQFT_B = x
    call sinqb(n, DclSinQFT_B, wk(idx)%array)
    call prccls('DclSinQFT_B')
  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  !作業領域番号

    call prcopn('DclInitCosQFT')
    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)
    call prccls('DclInitCosQFT')
  end subroutine
!---------------------------------------------------------
  subroutine DclDeallocCosQFT(index)
    integer, intent(in), optional :: index  !作業領域番号

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

    call prcopn('DclCosQFT_F')
    if(present(index)) then; idx = index
                       else; idx = 1     ; end if

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

    call prcopn('DclCosQFT_B')
    if(present(index)) then; idx = index
                       else; idx = 1     ; end if

    n = size(x)
    if(n.ne.wk(idx)%n) call msgdmp('E', 'DclCosQFT_B', &
                   & 'Wrong working area.')
    DclCosQFT_B = x
    call cosqb(size(x), DclCosQFT_B, wk(idx)%array)
    call prccls('DclCosQFT_B')
  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  !作業領域番号

    call prcopn('DclInitComplexFFT')
    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)
    call prccls('DclInitComplexFFT')
  end subroutine
!---------------------------------------------------------
  subroutine DclDeallocComplexFFT(index)
    integer, intent(in), optional :: index  !作業領域番号

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

    call prcopn('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)
    call prccls('DclComplexFFT_F')
  end function
!---------------------------------------------------------
  function DclComplexFFT_B(c,index)
    complex, intent(in), dimension(:) :: c
    integer, intent(in), optional     :: index
    complex, dimension(size(c))       :: DclComplexFFT_B

    call prcopn('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)
    call prccls('DclComplexFFT_B')
  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
