!----------------------------------------------------------------------
!   COPYRIGHT (c) 2006 Shin-ichi Takehiro. All rights reserved.
!----------------------------------------------------------------------
!
!表題  at_aq_galerkin_RRFF2
!
!      チェビシェフ−ガラーキン法
!      非圧縮流体の流線関数・流れポテンシャル用
!
!      case RR: 両端粘着条件(両端で値と 1 階微分が 0)
!               f(:,i=0)=f(:,i=im)=f'(:,i=0)=f'(:,i=im)=0 
!               [ f(:,x=xmax)=f(:,x=xmin)=f'(:,x=xmax)=0=f'(:,x=xmin)=0 ]
!
!      case FF: 両端自由すべり条件(両端で値と 2 階微分が 0)
!               f(:,i=0)=f(:,i=im)=f''(:,i=0)=f''(:,i=im)=0 
!               [ f(:,x=xmax)=f(:,x=xmin)=f''(:,x=xmax)=0=f''(:,x=xmin)=0 ]
!
!      case RF: 片端粘着条件, 片端自由すべり条件
!               (両端で値が 0, 片側で 1 階微分, もう一方で 2 階微分が 0)
!               f(:,i=0)=f(:,i=im)=f'(:,i=0)=f''(:,i=im)=0 
!               [ f(:,x=xmax)=f(:,x=xmin)=f''(:,x=xmax)=0=f'(:,x=xmin)=0 ]
!
!      case FR: 片端自由すべり条件, 片端粘着条件
!               (両端で値が 0, 片側で 2 階微分, もう一方で 1 階微分が 0)
!               f(:,i=0)=f(:,i=im)=f''(:,i=0)=f'(:,i=im)=0 
!               [ f(:,x=xmax)=f(:,x=xmin)=f'(:,x=xmax)=0=f''(:,x=xmin)=0 ]
!
!      k 次のガラーキン基底が
!        \phi_k(x)=T_k(x) + C_3 T_3(x) + C_2 T_2(x) 
!                         + C_1 T_1(x) + C_0 T_0(x) 
!      の型式(Type2)
!
!履歴  2006/01/24  竹広真一  新規作成
!      2006/01/25  竹広真一  モジュール初期化フラッグ導入
!      2010/09/21  竹広真一  at_module と境界条件スイッチを揃えた
!
module at_aq_galerkin_RRFF

  use lumatrix
  use at_module
  use dc_message

  implicit none

  private
  public :: at_aq_galerkin_RRFF_Initial  ! 初期化サブルーチン
  public :: aq_at, q_t                   ! チェビシェフ->ガラーキン変換
  public :: aq_ag, q_g                   ! 格子点->ガラーキンデータ変換
  public :: at_aq, t_q                   ! ガラーキン->チェビシェフ変換
  public :: ag_aq, g_q                   ! ガラーキン->格子点データ変換
  public :: aq_Dx_aq, q_Dx_q             ! X 微分

  real(8), allocatable :: TQ(:,:)   ! ガラーキン->チェビシェフ変換行列
  real(8), allocatable :: QT(:,:)   ! チェビシェフ->ガラーキン変換行列
  integer, allocatable :: kp(:)     ! LU 分解用ピボット格納配列

  real(8), allocatable :: alpha(:)  ! チェビシェフ<->ガラーキン行列積係数
  real(8), allocatable :: beta(:)   ! チェビシェフ<->ガラーキン行列和係数

  integer :: im                     ! 格子点数              
  integer :: km                     ! チェビシェフ切断波数  
  integer,parameter :: ks=4         ! ガラーキン基底最低次数
  character(LEN=2) :: bc            ! 境界条件(RR/FF/RF/FR)

  logical :: RRFF_Initialized=.false.    ! 初期化フラッグ

  save im, km, bc, TQ, QT, kp, alpha, beta, RRFF_Initialized
  public TQ, QT, alpha, beta

contains

  !
  ! チェビシェフ−ガラーキン法 
  ! 両端粘着条件用モジュール
  ! 初期化サブルーチン
  !
  subroutine at_aq_galerkin_RRFF_Initial(i_in,k_in,bc_in)

    integer, intent(IN)          :: i_in     ! (in)格子点数
    integer, intent(IN)          :: k_in     ! (in)チェビシェフ切断波数  
    character(LEN=2), intent(IN) :: bc_in    ! (in)境界条件(RR/FF/RF/FR)

    integer :: k, l, m, n

    im=i_in ; km=k_in ; bc=bc_in

    allocate(TQ(0:km,ks:km),QT(ks:km,ks:km),kp(ks:km))
    allocate(alpha(0:km),beta(0:km))


    select case (bc)

    !---- 両端粘着条件用変換行列設定
    case ('RR','rr')
       TQ = 0.0D0
       do k=ks,km
          TQ(0,k) =  ( (-1)**k + 1 ) * (k**2-4)/8.0D0
          TQ(1,k) = -( 1 - (-1)**k ) * (9-k**2)/16.0D0
          TQ(2,k) = -( (-1)**k + 1 ) * k**2    /8.0D0
          TQ(3,k) =  ( 1 - (-1)**k ) * (1-k**2)/16.0D0
          TQ(k,k) = 1.0
       enddo
       call MessageNotify('M','at_aq_galerkin_RRFF_Initial(Type2)',&
                          'Conversion matrices for RR-B.C. initialized.')

    !---- 両端自由すべり条件用変換行列設定
    case ('FF','ff','SS','ss')
       TQ = 0.0D0
       do k=ks,km
          TQ(0,k) =  ( (-1)**k + 1 ) * ( k**2*(k**2-1)-12 )/24.0D0
          TQ(1,k) =  ( 1 - (-1)**k ) * ( k**2*(k**2-1)-72 )/144.0D0
          TQ(2,k) = -( (-1)**k + 1 ) * k**2*(k**2-1)       /24.0D0
          TQ(3,k) = -( 1 - (-1)**k ) * k**2*(k**2-1)       /144.0D0
          TQ(k,k) = 1.0
       enddo
       call MessageNotify('M','at_aq_galerkin_RRFF_Initial(Type2)',&
                          'Conversion matrices for FF-B.C. initialized.')


    !---- 片端粘着条件, 片端自由すべり条件変換行列設定
    case ('RF','rf','RS','rs')
       TQ = 0.0D0
       do k=ks,km
          TQ(0,k) =      (-1)**k * k**2*(k**2-1)/48.0D0         &
                       + 3.0D0/16.0D0 * k**2                    &
                       - 3.0D0/32.0D0 * ( 1 - (-1)**k )         &
                       - 1.0D0/2.0D0  * ( 1 + (-1)**k )

          TQ(1,k) =  - (-1)**k * k**2*(k**2-1)/96.0D0           &
                     + 1.0D0/32.0D0 * k**2                      &
                     - 33.0D0/64.0D0 * ( 1 - (-1)**k )

          TQ(2,k) =  - (-1)**k * k**2*(k**2-1)/48.0D0           &
                     - 3.0D0/16.0D0 * k**2                      &
                     + 3.0D0/32.0D0 * ( 1 - (-1)**k )

          TQ(3,k) = 1.0D0/32.0D0 * (    (-1)**k * k**2*(k**2-1)/3.0D0 &
                                      - k**2                          &
                                      + ( 1 - (-1)**k )/2.0D0 )
          TQ(k,k) = 1.0
       enddo
       call MessageNotify('M','at_aq_galerkin_RRFF_Initial(Type2)',&
                          'Conversion matrices for RF-B.C. initialized.')

    !---- 片端自由すべり条件, 片端粘着条件変換行列設定
    case ('FR','fr','SR','sr')
       TQ = 0.0D0
       do k=ks,km
          TQ(0,k) =      k**2*(k**2-1)/48.0D0                   &
                       - 3.0D0/16.0D0 * (-1)**(k+1) * k**2      &
                       + 3.0D0/32.0D0 * ( 1 - (-1)**k )         &
                       - 1.0D0/2.0D0  * ( 1 + (-1)**k )

          TQ(1,k) =    k**2*(k**2-1)/96.0D0                     &
                     + 1.0D0/32.0D0 * (-1)**(k+1) * k**2        &
                     - 33.0D0/64.0D0 * ( 1 - (-1)**k )

          TQ(2,k) =  - k**2*(k**2-1)/48.0D0                     &
                     + 3.0D0/16.0D0 * (-1)**(k+1) * k**2        &
                     - 3.0D0/32.0D0 * ( 1 - (-1)**k )

          TQ(3,k) = -1.0D0/32.0D0 * (   k**2*(k**2-1)/3.0D0     &
                                      + (-1)**(k+1) * k**2      &
                                      - ( 1 - (-1)**k )/2.0D0 )
          TQ(k,k) = 1.0
       enddo
       call MessageNotify('M','at_aq_galerkin_RRFF_Initial(Type2)',&
                          'Conversion matrices for FR-B.C. initialized.')

    case default
       call MessageNotify('E','at_aq_galerkin_RRFF_Initial(Type2)',&
                'Argument for B.C. not valid. Should be RR/FF/RF/FR.')

    end select

    beta=1.0 ; beta(0)=0.5D0
    if (im .eq. km ) beta(km)=0.5D0

    ! 両端粘着条件用変換逆行列
    alpha=1.0 ; alpha(0)=2.0D0

    QT = 0.0D0
    do m=ks,km
       do n=ks,km
          do l=0,km
             QT(m,n) = QT(m,n) + alpha(l)*TQ(l,m)*TQ(l,n)
          enddo
       enddo
    enddo

    call LUDecomp(QT,kp)

    RRFF_Initialized=.true.

  end subroutine at_aq_galerkin_RRFF_Initial

  !
  ! 両端粘着条件
  ! チェビシェフ係数 -> ガラーキン係数変換(2次元データ)
  !
  function aq_at(at_data)
    real(8), intent(IN) :: at_data(:,0:)                !(in)  チェビシェフ係数
    real(8)             :: aq_at(size(at_data,1),ks:km) !(out) ガラーキン係数  

    real(8)             :: aq_work(size(at_data,1),ks:km)  ! 作業用配列
    
    integer :: k,m

    if ( .not. RRFF_Initialized ) &
         call MessageNotify('E','aq_at',&
                            'at_aq_galerkin_RRFF_module(Type2) not initialized')

    aq_work =0.0
    do m=ks,km
       do k=0,km
          aq_work(:,m) = aq_work(:,m) &
               + alpha(k) * beta(k) * at_data(:,k) * TQ(k,m)
       enddo
    enddo

    aq_at = LUSolve(QT,kp,aq_work)
  end function aq_at

  !
  ! 両端粘着条件
  ! チェビシェフ係数 -> ガラーキン係数変換(1次元データ)
  !
  function q_t(t_data)

    real(8), intent(IN) :: t_data(0:km)        !(in)  チェビシェフ係数
    real(8)             :: q_t(ks:km)          !(out) ガラーキン係数  

    real(8)             :: q_work(ks:km)       ! 作業用配列

    integer :: k,m

    if ( .not. RRFF_Initialized ) &
         call MessageNotify('E','q_t',&
                            'at_aq_galerkin_RRFF_module(Type2) not initialized')

    q_work =0.0
    do m=ks,km
       do k=0,km
          q_work(m) = q_work(m) &
               + alpha(k) * beta(k) * t_data(k) * TQ(k,m)
       enddo
    enddo

    q_t = LUSolve(QT,kp,q_work)
  end function q_t

  !
  ! 両端粘着条件
  ! ガラーキン係数 -> チェビシェフ係数変換(2次元データ)
  !
  function at_aq(aq_data)

    real(8), intent(IN)  :: aq_data(:,ks:)              !(in)  ガラーキン係数
    real(8)              :: at_aq(size(aq_data,1),0:km) !(out) チェビシェフ係数

    integer :: m, n

    if ( .not. RRFF_Initialized ) &
         call MessageNotify('E','at_aq',&
                            'at_aq_galerkin_RRFF_module(Type2) not initialized')
    at_aq = 0.0D0
    do m=0,km
       do n=ks,km
          at_aq(:,m) = at_aq(:,m) + TQ(m,n)*aq_data(:,n)/beta(m)
       enddo
    enddo
  end function at_aq

  !
  ! 両端粘着条件
  ! ガラーキン係数 -> チェビシェフ係数変換(1次元データ)
  !
  function t_q(q_data)

    real(8), intent(IN)  :: q_data(ks:km)       !(in)  ガラーキン係数  
    real(8)              :: t_q(0:km)           !(out) チェビシェフ係数

    integer :: m, n

    if ( .not. RRFF_Initialized ) &
         call MessageNotify('E','q_t',&
                            'at_aq_galerkin_RRFF_module(Type2) not initialized')

    t_q = 0.0D0
    do m=0,km
       do n=ks,km
          t_q(m) = t_q(m) + TQ(m,n)*q_data(n)/beta(m)
       enddo
    enddo
  end function t_q

  !
  ! 両端粘着条件
  ! 格子点データ -> ガラーキン係数変換(2次元データ)
  !
  function aq_ag(ag_data)
    real(8), intent(IN)  :: ag_data(:,0:)                !(in)  格子点データ
    real(8)              :: aq_ag(size(ag_data,1),ks:km) !(out) ガラーキン係数

    aq_ag = aq_at(at_ag(ag_data))
  end function aq_ag

  !
  ! 両端粘着条件
  ! 格子点データ -> ガラーキン係数変換(1次元データ)
  !
  function q_g(g_data)
    real(8), intent(IN)  :: g_data(0:im)        !(in)  格子点データ
    real(8)              :: q_g(ks:km)          !(out) ガラーキン係数

    q_g = q_t(t_g(g_data))
  end function q_g

  !
  ! 両端粘着条件
  ! ガラーキン係数 -> 格子点データ変換(2次元データ)
  !
  function ag_aq(aq_data)
    real(8), intent(IN) :: aq_data(:,ks:)              !(in)  ガラーキン係数  
    real(8)             :: ag_aq(size(aq_data,1),0:im) !(out) 格子点データ  
    
    ag_aq = ag_at(at_aq(aq_data))
  end function ag_aq

  !
  ! 両端粘着条件
  ! ガラーキン係数 -> 格子点データ変換(1次元データ)
  !
  function g_q(q_data)
    real(8), intent(IN) :: q_data(ks:km)        !(in)  ガラーキン係数  
    real(8)             :: g_q(0:im)            !(out) 格子点データ  
    
    g_q = g_t(t_q(q_data))
  end function g_q

  !
  ! 両端粘着条件
  ! X 微分計算(1 次元)
  !
  function aq_Dx_aq(aq_data)
    real(8), intent(IN) :: aq_data(:,ks:)                  !(in) ガラーキン係数
    real(8)             :: aq_Dx_aq(size(aq_data,1),ks:km) !(out) 微分ガラーキン
    aq_Dx_aq = aq_at(at_Dx_at(at_aq(aq_data)))
  end function aq_Dx_aq

  !
  ! 両端粘着条件
  ! X 微分計算(1 次元)
  !
  function q_Dx_q(q_data)
    real(8), intent(IN) :: q_data(ks:km)
    real(8)             :: q_Dx_q(ks:km)

    q_Dx_q = q_t(t_Dx_t(t_q(q_data)))

  end function q_Dx_q

end module at_aq_galerkin_RRFF
