!----------------------------------------------------------------------
!     Copyright (c) 2002 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  FTPACK/ISPACK-F90
!
!履歴  2002/01/25 竹広真一 
!      
!
module ftrpack
  implicit none

  private
  public ftrinitial
  public g_e, e_g, e_dx_e
  public g_x

  interface g_e                        ! Fourier 逆変換
     module procedure g_e_2d, g_e_1d
  end interface

  interface e_g                        ! Fourier 正変換
     module procedure e_g_2d, e_g_1d
  end interface

  interface e_dx_e                     ! x 微分
     module procedure e_dx_e_2d, e_dx_e_1d
  end interface

  integer            :: im=32             ! 格子点の設定
  integer            :: km=10             ! 切断波数の設定
  double precision   :: xl=1.0            ! 領域の大きさ

  integer,dimension(5)              :: iti
  real(8),dimension(:),allocatable  :: ti
  real(8), parameter                :: pi=3.1415926535897932385D0

  real(8), allocatable :: g_x(:)          ! 格子点座標

  save im, km, iti, ti, xl, g_x

  contains
  !--------------- 初期化 -----------------
    subroutine ftrinitial(i,k,xlength)

      integer,intent(in) :: i           ! 格子点の設定
      integer,intent(in) :: k           ! 切断波数の設定

      real(8),intent(in)  :: xlength    ! 領域の大きさ

      integer :: ii

      im = i
      km = k
      xl = xlength

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

      allocate(ti(im*2))

      call fttrui(im,iti,ti)

      allocate(g_x(0:im-1))
      do ii=0,im
         g_x(ii) = xlength/im*ii
      enddo

    end subroutine ftrinitial

  !--------------- 基本変換 -----------------

    function g_e_2d(e) ! スペクトル -> 格子
      real(8), dimension(:,-km:), intent(in) :: e
      real(8), dimension(size(e,1),0:im-1)   :: g_e_2d

      real(8), dimension(size(e,1)*im)       :: y
      integer :: m, k

      m=size(e,1)
      if ( size(e,2) < 2*km+1 ) then
         call msgdmp('E','g_e', &
              'The Fourier dimension of input data too small.')
      elseif ( size(e,2) > 2*km+1 ) then
         call msgdmp('W','g_e', &
              'The Fourier dimension of input data too large.')
      endif

      g_e_2d = 0.0D0
      g_e_2d(:,0)=e(:,0)
      g_e_2d(:,1)=0
      do k=1,km
         g_e_2d(:,2*k)=e(:,k)
         g_e_2d(:,2*k+1)=e(:,-k)
      enddo
      g_e_2d(:,2*km+2:im-1)=0

      call fttrub(m,im,g_e_2d,y,iti,ti)
    end function g_e_2d

    function g_e_1d(e) ! スペクトル -> 格子
      real(8), dimension(0:im-1)             :: g_e_1d
      real(8), dimension(-km:km), intent(in) :: e

      real(8), dimension(1,size(e))  :: e_work
      real(8), dimension(1,0:im-1)   :: g_work

      e_work(1,:) = e
      g_work = g_e_2d(e_work)
      g_e_1d = g_work(1,:)

    end function g_e_1d

    function e_g_2d(g)  ! 格子 -> スペクトル
      real(8), dimension(:,:), intent(in)     :: g
      real(8), dimension(size(g,1),-km:km)    :: e_g_2d

      real(8), dimension(size(g,1)*im)     :: y
      real(8), dimension(size(g,1),0:im-1) :: g_work
      integer :: m, k

      m = size(g,1)
      if ( size(g,2) < im ) then
         call msgdmp('E','e_g', &
              'The Grid points of input data too small.')
      elseif ( size(g,2) > im ) then
         call msgdmp('W','e_g', &
              'The Grid points of input data too large.')
      endif
      g_work = g

      call fttruf(m,im,g_work,y,iti,ti)

      do k=1,km
         e_g_2d(:,k) = g_work(:,2*k)
         e_g_2d(:,-k) = g_work(:,2*k+1)
      enddo
      e_g_2d(:,0) = g_work(:,0)

    end function e_g_2d

    function e_g_1d(g)  ! 格子 -> スペクトル
      real(8), dimension(-km:km)              :: e_g_1d
      real(8), dimension(0:im-1), intent(in)  :: g

      real(8), dimension(1,size(g))        :: g_work
      real(8), dimension(1,-km:km)         :: e_work

      g_work(1,:) = g
      e_work = e_g_2d(g_work)
      e_g_1d = e_work(1,:)

    end function e_g_1d

  !--------------- 微分計算 -----------------
    function e_dx_e_2d(e)   ! スペクトルに作用する x 微分演算子

      real(8), dimension(:,-km:), intent(in)     :: e
      real(8), dimension(size(e,1),-km:km)       :: e_dx_e_2d

      integer k

      if ( size(e,2) < 2*km+1 ) then
         call msgdmp('W','e_dx_e', &
              'The Chebyshev dimension of input data too small.')
      elseif ( size(e,2) > 2*km+1 ) then
         call msgdmp('W','e_dx_e', &
              'The Chebyshev dimension of input data too large.')
      endif

      do k=-km,km
         e_dx_e_2d(:,k) = -(2*pi*k/xl)*e(:,-k)
      enddo
    end function e_dx_e_2d

    function e_dx_e_1d(e)   ! スペクトルに作用する x 微分演算子

      real(8), dimension(-km:km), intent(in)     :: e
      real(8), dimension(-km:km)                 :: e_dx_e_1d

      real(8), dimension(1,-km:km)               :: e_work

      e_work(1,:) = e
      e_work = e_dx_e_2d(e_work)
      e_dx_e_1d = e_work(1,:)

    end function e_dx_e_1d

  end module ftrpack
