!----------------------------------------------------------------------
!     Copyright (c) 2002 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  CTPACK/ISPACK-F90
!      Chebyshev 展開
!
!履歴  2002/01/24  竹広真一
!
module ctpack

  use lumatrix
  implicit none
  private 
  public g_x
  public ctinitial
  public g_t, t_g, t_dx_t
  public ctboundaries_DD, ctboundaries_DN, ctboundaries_ND, ctboundaries_NN

  interface g_t                        ! チェビシェフ逆変換
     module procedure g_t_2d, g_t_1d
  end interface

  interface t_g                        ! チェビシェフ正変換
     module procedure t_g_2d, t_g_1d
  end interface

  interface t_dx_t                     ! チェビシェフ微分
     module procedure t_dx_t_2d, t_dx_t_1d
  end interface

  interface ctboundaries_DD            ! 両端ディリクレ条件
     module procedure ctboundaries_DD_1d, ctboundaries_DD_2d
  end interface

  interface ctboundaries_DN            ! 始点ディリクレ終点ノイマン
     module procedure ctboundaries_DN_1d, ctboundaries_DN_2d
  end interface

  interface ctboundaries_ND            ! 始点ノイマン終点ディリクレ
     module procedure ctboundaries_ND_1d, ctboundaries_ND_2d
  end interface

  interface ctboundaries_NN            ! 両端ノイマン条件
     module procedure ctboundaries_NN_1d, ctboundaries_NN_2d
  end interface

  integer, dimension(5)              :: it
  real(8), dimension(:), allocatable :: t
  real(8), parameter                 :: pi=3.1415926535897932385D0

  integer :: im, km                 ! 格子点数, 波数
  real(8) :: xl                     ! 領域の大きさ

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

  save :: im, km, xl, it, t, g_x

contains

! ---- 初期化 ---- 
  subroutine ctinitial(i,k,xmin,xmax)
    integer, intent(in) :: i, k
    real(8), intent(in) :: xmin,xmax

    integer :: ii

    im=i ; km=k 
    xl = xmax-xmin

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

    allocate(t(3*im))
    call fttcti(im,it,t)

    allocate(g_x(0:im))
    do ii=0,im
       g_x(ii) = (xmax+xmin)/2 + xl/2 * cos(pi*ii/im)
    enddo

  end subroutine ctinitial

! ---- 逆変換 ---- 
  function g_t_2d(t_data)  ! スペクトル -> 台形格子

    double precision, dimension(:,:), intent(in)      :: t_data
    double precision, dimension(size(t_data,1),0:im)  :: g_t_2d

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

    m = size(t_data,1)
    if ( size(t_data,2)-1 < km ) then
       call msgdmp('E','g_t', &
            'The Chebyshev dimension of input data too small.')
    elseif ( size(t_data,2)-1 > km ) then
       call msgdmp('W','g_t', &
            'The Chebyshev dimension of input data too large.')
    endif

    g_t_2d = 0.0D0
    g_t_2d(:,0:km)=t_data
    call fttctb(m,im,g_t_2d,y,it,t)

  end function g_t_2d

  function g_t_1d(t_data)  ! スペクトル -> 台形格子
    double precision, dimension(:), intent(in)  :: t_data
    double precision, dimension(0:im)           :: g_t_1d

    double precision, dimension(1,size(t_data)) :: t_work
    double precision, dimension(1,0:im)         :: g_work

    t_work(1,:) = t_data  
    g_work = g_t_2d(t_work)
    g_t_1d = g_work(1,:)

  end function g_t_1d


! ---- 正変換 ---- 
  function t_g_2d(g_data)  ! 台形格子 -> スペクトル
    double precision, dimension(:,:), intent(in)     :: g_data
    double precision, dimension(size(g_data,1),0:km) :: t_g_2d

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

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

    call fttctf(m,im,g_work,y,it,t)
    t_g_2d = g_work(:,0:km)

  end function t_g_2d

  function t_g_1d(g_data)  ! 台形格子 -> スペクトル
    double precision, dimension(:), intent(in)     :: g_data
    double precision, dimension(0:km)              :: t_g_1d

    real(8), dimension(1,size(g_data)) :: g_work
    real(8), dimension(1,0:km)         :: t_work

    g_work(1,:) = g_data
    t_work = t_g_2d(g_work)
    t_g_1d = t_work(1,:)
    
  end function t_g_1d

! ---- 微分計算 ---- 
  function t_dx_t_2d(t_data)
    real(8), dimension(:,0:), intent(in)                  :: t_data
    real(8), dimension(size(t_data,1),0:size(t_data,2)-1) :: t_dx_t_2d

    integer :: m, k
    integer :: nm, kmax

    nm=size(t_data,1)
    kmax=size(t_data,2)-1
    if ( kmax  < km ) then
       call msgdmp('W','t_dx_t', &
            'The Chebyshev dimension of input data too small.')
    elseif ( kmax > km ) then
       call msgdmp('E','t_dx_t', &
            'The Chebyshev dimension of input data too large.')
    endif

    if ( kmax == im ) then
       do m=1,nm
          t_dx_t_2d(m,kmax)   = 0. 
          t_dx_t_2d(m,kmax-1) = 2 * km * t_data(m,kmax) /2
       enddo
    else
       do m=1,nm
          t_dx_t_2d(m,kmax)   = 0. 
          t_dx_t_2d(m,kmax-1) = 2 * km * t_data(m,kmax) 
          ! スタートはグリッド対応最大波数未満. Factor 1/2 不要
       enddo
    endif

    do k=kmax-2,0,-1
       do m=1,nm
          t_dx_t_2d(m,k) = t_dx_t_2d(m,k+2) + 2*(k+1)*t_data(m,k+1)
       enddo
    enddo

    do k=0,kmax
       do m=1,nm
          t_dx_t_2d(m,k) = 2/xl * t_dx_t_2d(m,k)
       enddo
    enddo

  end function t_dx_t_2d

  function t_dx_t_1d(t_data)
    real(8), dimension(:), intent(in)   :: t_data
    real(8), dimension(size(t_data))    :: t_dx_t_1d

    real(8), dimension(1,size(t_data))  :: t_work
    
    t_work(1,:) = t_data
    t_work = t_dx_t_2d(t_work)
    t_dx_t_1d = t_work(1,:)

  end function t_dx_t_1d

!---- Dirichlet 型境界条件 ----

  subroutine ctboundaries_DD_2d(t_data,values)
   ! Dirichlet 型境界条件の適用

    real(8), dimension(:,0:),intent(inout)         :: t_data  ! データ(m,0:km)
    real(8), dimension(:,:), intent(in), optional  :: values  ! 境界値(m,2)

    real(8), dimension(:,:), allocatable  :: alu
    integer, dimension(:), allocatable    :: kp
    real(8), dimension(0:km,0:km)         :: cheb
    real(8), dimension(0:km,0:im)         :: g
    real(8), dimension(size(t_data,1))    :: value1, value2           ! 境界値

    logical :: first = .true.
    integer :: k
    save    :: alu, kp, first

    if ( size(t_data,2)-1 < km ) then
       call msgdmp('E','ctboundaries_DD', &
            'The Chebyshev dimension of input data too small.')
    elseif ( size(t_data,2)-1 > km ) then
       call msgdmp('W','ctboundaries_DD', &
            'The Chebyshev dimension of input data too large.')
    endif

    if (.not. present(values)) then
       value1=0 ; value2=0
    else
       value1 = values(:,1) ; value2 = values(:,2)
    endif

    if ( first ) then
       first = .false.

       allocate(alu(0:km,0:km),kp(0:km))

       cheb=0
       do k=0,km
          cheb(k,k)=1
       enddo
       alu = cheb

       g=g_t(cheb)
       alu(km-1,:) = g(:,0)
       alu(km,:)   = g(:,im)

       call ludecomp(alu,kp)
    endif

    t_data(:,km-1) = value1
    t_data(:,km)   = value2
    t_data = lusolve(alu,kp,t_data)

  end subroutine ctboundaries_DD_2d

  subroutine ctboundaries_DD_1d(t_data,values)
   ! Dirichlet 型境界条件の適用

    real(8), dimension(0:km),intent(inout)       :: t_data
    real(8), dimension(2), intent(in), optional  :: values          ! 境界値

    real(8), dimension(1,0:km)                   :: t_work
    real(8), dimension(1,2)                      :: vwork           ! 境界値

    if (.not. present(values)) then
       vwork(1,1)=0 ; vwork(1,2)=0
    else
       vwork(1,:) = values
    endif

    t_work(1,:)=t_data
    call ctboundaries_DD_2d(t_work,vwork)
    t_data=t_work(1,:)

  end subroutine ctboundaries_DD_1d


!---- Dirichlet/Neumann 型境界条件 ----

  subroutine ctboundaries_DN_2d(t_data,values)
   ! Dirichlet/Neumann 型境界条件の適用

    real(8), dimension(:,0:),intent(inout)         :: t_data  ! データ(m,0:km)
    real(8), dimension(:,:), intent(in), optional  :: values  ! 境界値(m,2)

    real(8), dimension(:,:), allocatable  :: alu
    integer, dimension(:), allocatable    :: kp
    real(8), dimension(0:km,0:km)         :: cheb
    real(8), dimension(0:km,0:im)         :: g
    real(8), dimension(size(t_data,1))    :: value1, value2           ! 境界値

    logical :: first = .true.
    integer :: k
    save    :: alu, kp, first

    if ( size(t_data,2)-1 < km ) then
       call msgdmp('E','ctboundaries_DN', &
            'The Chebyshev dimension of input data too small.')
    elseif ( size(t_data,2)-1 > km ) then
       call msgdmp('W','ctboundaries_DN', &
            'The Chebyshev dimension of input data too large.')
    endif

    if (.not. present(values)) then
       value1=0 ; value2=0
    else
       value1 = values(:,1) ; value2 = values(:,2)
    endif

    if ( first ) then
       first = .false.

       allocate(alu(0:km,0:km),kp(0:km))

       cheb=0
       do k=0,km
          cheb(k,k)=1
       enddo
       alu = cheb

       g=g_t(cheb)
       alu(km-1,:) = g(:,0)
       g=g_t(t_dx_t(cheb))
       alu(km,:)   = g(:,im)

       call ludecomp(alu,kp)
    endif

    t_data(:,km-1) = value1
    t_data(:,km)   = value2
    t_data = lusolve(alu,kp,t_data)

  end subroutine ctboundaries_DN_2d

  subroutine ctboundaries_DN_1d(t_data,values)
   ! Dirichlet/Neumann 型境界条件の適用

    real(8), dimension(0:km),intent(inout)       :: t_data
    real(8), dimension(2), intent(in), optional  :: values          ! 境界値

    real(8), dimension(1,0:km)                   :: t_work
    real(8), dimension(1,2)                      :: vwork           ! 境界値

    if (.not. present(values)) then
       vwork(1,1)=0 ; vwork(1,2)=0
    else
       vwork(1,:) = values
    endif

    t_work(1,:)=t_data
    call ctboundaries_DN_2d(t_work,vwork)
    t_data=t_work(1,:)

  end subroutine ctboundaries_DN_1d

!---- Neumann/Dirichlet 型境界条件 ----

  subroutine ctboundaries_ND_2d(t_data,values)
   ! Neumann/Dirichlet 型境界条件の適用

    real(8), dimension(:,0:),intent(inout)         :: t_data  ! データ(m,0:km)
    real(8), dimension(:,:), intent(in), optional  :: values  ! 境界値(m,2)

    real(8), dimension(:,:), allocatable  :: alu
    integer, dimension(:), allocatable    :: kp
    real(8), dimension(0:km,0:km)         :: cheb
    real(8), dimension(0:km,0:im)         :: g
    real(8), dimension(size(t_data,1))    :: value1, value2           ! 境界値

    logical :: first = .true.
    integer :: k
    save    :: alu, kp, first

    if ( size(t_data,2)-1 < km ) then
       call msgdmp('E','ctboundaries_ND', &
            'The Chebyshev dimension of input data too small.')
    elseif ( size(t_data,2)-1 > km ) then
       call msgdmp('W','ctboundaries_ND', &
            'The Chebyshev dimension of input data too large.')
    endif

    if (.not. present(values)) then
       value1=0 ; value2=0
    else
       value1 = values(:,1) ; value2 = values(:,2)
    endif

    if ( first ) then
       first = .false.

       allocate(alu(0:km,0:km),kp(0:km))

       cheb=0
       do k=0,km
          cheb(k,k)=1
       enddo
       alu = cheb

       g=g_t(t_dx_t(cheb))
       alu(km-1,:) = g(:,0)
       g=g_t(cheb)
       alu(km,:)   = g(:,im)

       call ludecomp(alu,kp)
    endif

    t_data(:,km-1) = value1
    t_data(:,km)   = value2
    t_data = lusolve(alu,kp,t_data)

  end subroutine ctboundaries_ND_2d

  subroutine ctboundaries_ND_1d(t_data,values)
   ! Neumann/Dirichlet 型境界条件の適用

    real(8), dimension(0:km),intent(inout)       :: t_data
    real(8), dimension(2), intent(in), optional  :: values          ! 境界値

    real(8), dimension(1,0:km)                   :: t_work
    real(8), dimension(1,2)                      :: vwork           ! 境界値

    if (.not. present(values)) then
       vwork(1,1)=0 ; vwork(1,2)=0
    else
       vwork(1,:) = values
    endif

    t_work(1,:)=t_data
    call ctboundaries_ND_2d(t_work,vwork)
    t_data=t_work(1,:)

  end subroutine ctboundaries_ND_1d


!---- Neumann 型境界条件 ----

  subroutine ctboundaries_NN_2d(t_data,values)
   ! Neumann/Dirichlet 型境界条件の適用

    real(8), dimension(:,0:),intent(inout)         :: t_data  ! データ(m,0:km)
    real(8), dimension(:,:), intent(in), optional  :: values  ! 境界値(m,2)

    real(8), dimension(:,:), allocatable  :: alu
    integer, dimension(:), allocatable    :: kp
    real(8), dimension(0:km,0:km)         :: cheb
    real(8), dimension(0:km,0:im)         :: g
    real(8), dimension(size(t_data,1))    :: value1, value2           ! 境界値

    logical :: first = .true.
    integer :: k
    save    :: alu, kp, first

    if ( size(t_data,2)-1 < km ) then
       call msgdmp('E','ctboundaries_NN', &
            'The Chebyshev dimension of input data too small.')
    elseif ( size(t_data,2)-1 > km ) then
       call msgdmp('W','ctboundaries_NN', &
            'The Chebyshev dimension of input data too large.')
    endif

    if (.not. present(values)) then
       value1=0 ; value2=0
    else
       value1 = values(:,1) ; value2 = values(:,2)
    endif

    if ( first ) then
       first = .false.

       allocate(alu(0:km,0:km),kp(0:km))

       cheb=0
       do k=0,km
          cheb(k,k)=1
       enddo
       alu = cheb

       g=g_t(t_dx_t(cheb))
       alu(km-1,:) = g(:,0)
       alu(km,:)   = g(:,im)

       call ludecomp(alu,kp)
    endif

    t_data(:,km-1) = value1
    t_data(:,km)   = value2
    t_data = lusolve(alu,kp,t_data)

  end subroutine ctboundaries_NN_2d

  subroutine ctboundaries_NN_1d(t_data,values)
   ! Neumann/Dirichlet 型境界条件の適用

    real(8), dimension(0:km),intent(inout)       :: t_data
    real(8), dimension(2), intent(in), optional  :: values          ! 境界値

    real(8), dimension(1,0:km)                   :: t_work
    real(8), dimension(1,2)                      :: vwork           ! 境界値

    if (.not. present(values)) then
       vwork(1,1)=0 ; vwork(1,2)=0
    else
       vwork(1,:) = values
    endif

    t_work(1,:)=t_data
    call ctboundaries_NN_2d(t_work,vwork)
    t_data=t_work(1,:)

  end subroutine ctboundaries_NN_1d

end module ctpack
