!     
!     チェビシェフ多項式展開法(collocation法)を用いた
!     平面ポワズイユ流の成長率の計算.
!
program SolveEigenvalueProblemOfPoiseuilleFlow

  implicit none

  integer, parameter :: m = 257
  integer, parameter :: l = 100

  real(8) re, al
  real(8) y(1:m), ty(1:5, 1:m, 1:m), u(1:m), d2u(1:m)
  complex(8) ca(1:m, 1:m), cb(1:m, 1:m), eigenvector(1:m)
  integer i,j
  real(8) re_ary(1:l), al_ary(1:l)
  complex(8) lambda

  re_ary(1) = 1.0d3
  do i=2,l
     re_ary(i) = re_ary(i-1)*10.0d0**(6.0d0/dble(l))
  end do
  al_ary(1) = 0.0d0
  do i=2,l
     al_ary(i) = dble(i)*1.2d0/dble(l)
  end do

  do i=1,l
     do j=1,l
        re = re_ary(i)
        al = al_ary(j)
        call calculate_matrices(m, ty, u, d2u, re, al, ca, cb)
        call solve_eigenvalue_problem(m, ca, cb, eigenvector, lambda)
        print *, re, al, dble(lambda), dimag(lambda)
     end do
     print *, ""
  end do
  stop
end program SolveEigenvalueProblemOfPoiseuilleFlow

!
! y方向の級数展開の基底関数とy座標値のテーブルを作る
! 
subroutine create_tables(m, y, ty)

  implicit none

  integer, intent(in)  :: m
  real(8), intent(out) :: y(1:m), &           ! 格子点のy座標
       &                  ty(0:5, 1:m, 0:m-1) ! 基底関数

  real(8) tn(0:4, 1:m, 0:m-1)

  call set_axis(m, y)
  call create_chebychev_polynomial_table(m, y, tn)
  call create_basis_function_table(m, y, tn, ty)

end subroutine create_tables


!
! 格子点のy座標を設定する
! 
! y(i): i番目のインデックスの格子点のy座標
!
! ここで, y座標として Gauss-Lobatto grid を取る.
!
subroutine set_axis(m, y)

  implicit none

  integer, intent(in)  :: m
  real(8), intent(out) :: y(1:m)  ! 格子点のy座標

  integer j
  real(8) pi, dtheta

  pi = acos(-1.0d0)
  dtheta = pi/dble(m+1)
  do j=1,m
     y(j) = cos( dble(j)*dtheta )
  end do

end subroutine set_axis


!
! Chebychev多項式のテーブルを作る.
!
subroutine create_chebychev_polynomial_table(m, y, tn)

  implicit none

  integer, intent(in) :: m
  real(8), intent(in) :: y(1:m)
  real(8), intent(out) :: tn(0:4, 1:m, 0:m-1)
  ! tn(k,j,s)に y(j) における s 次のChebychev多項式の 
  ! k 階微分 d^k T_s/dy^k (y(j)) の値を格納する

  integer j, s
  
  ! d^k T_s/dy^k (y) の s=0,1 については予め代入しておく
  tn(0,1:m,0) = 1.0d0
  tn(0,1:m,1) = y(1:m)
  tn(1,1:m,0) = 0.0d0
  tn(1,1:m,1) = 1.0d0
  tn(2:4,1:m,0:1) = 0.0d0
  do s=2,m-1
     do j=1,m
        ! T_{n+1}(y) - 2yT_n(y) + T_{n-1}(y) = 0を用いた.
        tn(0,j,s) = 2.0d0*y(j)*tn(0,j,s-1) - tn(0,j,s-2)
        ! 昇降関係式 T_{n \pm 1}(y) = ( y \mp (1-y^2)/n d/dy ) T_n(y) を用いた.
        tn(1,j,s) = dble(s)*( tn(0,j,s-1) - y(j)*tn(0,j,s) ) / (1.0d0-y(j)**2)
        ! チェビシェフの微分方程式 (1-y^2) T_n'' - yT_n' + n^2T_n = 0を用いた.
        tn(2,j,s) = ( y(j)*tn(1,j,s) - dble(s)**2*tn(0,j,s) ) / (1.0d0-y(j)**2)
        ! チェビシェフの微分方程式をy微分したものを用いた.
        tn(3,j,s) = ( 3.0d0*y(j)*tn(2,j,s) - (dble(s)**2 - 1.0d0)*tn(1,j,s) ) &
             &                 / (1.0d0-y(j)**2)
        ! チェビシェフの微分方程式をyで2階微分したものを用いた.
        tn(4,j,s) = ( 5.0d0*y(j)*tn(3,j,s) - (dble(s)**2 - 4.0d0)*tn(2,j,s) ) &
             &                 / (1.0d0-y(j)**2)
     end do
  end do

end subroutine create_chebychev_polynomial_table


!
! Chebychev多項式を用いたy方向の展開に用いる基底関数の表を作る.
! 
! 基底関数は F_n(y) := (1-y^2)^2 T_n(y) とする.
! d^k F_n/dy^k (y)の値をk=0から4まで求める.
!
! Note: この基底関数は, Boyd(2001)のp.144にHeinrichs' basisとして載っている.
!
subroutine create_basis_function_table(m, y, tn, ty)

  implicit none

  integer, intent(in) :: m
  real(8), intent(in) :: y(1:m), tn(0:4, 1:m, 0:m-1)
  real(8), intent(out) :: ty(0:4, 1:m, 0:m-1) ! 基底関数

  integer j, s
  
  ! s 次の基底関数の k 階微分の y(j) での値 d^k F_s/dy^k (y(j)) の表 ty(k,j,s) を作る
  ! 以下の関係は基底関数を単にyでk階微分すると導出できる. 
  do s=0,m-1
     do j=1,m
        ty(0,j,s) = (1.0d0-y(j)**2)**2*tn(0,j,s)

        ty(1,j,s) = -4.0d0*y(j)*(1.0d0-y(j)**2)*tn(0,j,s) &
             &      +(1.0d0-y(j)**2)**2*tn(1,j,s)
        
        ty(2,j,s) = (-4.0d0+12.0d0*y(j)**2)*tn(0,j,s) &
             &      -8.0d0*y(j)*(1.0d0-y(j)**2)*tn(1,j,s) &
             &      +(1.0d0-y(j)**2)**2*tn(2,j,s)

        ty(3,j,s) = 24.0d0*y(j)*tn(0,j,s) &
             &      +(-12.0d0+36.0d0*y(j)**2)*tn(1,j,s) &
             &      -12.0d0*y(j)*(1.0d0-y(j)**2)*tn(2,j,s) &
             &      +(1.0d0-y(j)**2)**2*tn(3,j,s)

        ty(4,j,s) = 24.0d0*tn(0,j,s) &
             &      +96.0d0*y(j)*tn(1,j,s) &
             &      +(-24.0d0+72.0d0*y(j)**2)*tn(2,j,s) &
             &      -16.0d0*y(j)*(1.0d0-y(j)**2)*tn(3,j,s) &
                    +(1.0d0-y(j)**2)**2*tn(4,j,s)
     end do
  end do
  
end subroutine create_basis_function_table


!
! 平面ポワズイユ流の基本場を設定する
!
subroutine set_meanflow(m, y, u, d2u)
  
  implicit none
  
  integer, intent(in)  :: m         ! y方向のインデックスの最大値
  real(8), intent(in)  :: y(1:m)    ! 格子点の y 座標
  real(8), intent(out) :: u(1:m), & ! 東向きの流速, $ u = 1-y**2 $
       &                  d2u(1:m)  ! uのyに関する2階微分, $ d^2 u/dy^2 = -2 $

  u(1:m) = 1.0d0 - y(1:m)**2
  d2u(1:m) = -2.0d0
  
end subroutine set_meanflow


! 
! Orr-Sommerfeld方程式の解を基底関数で級数展開したときの展開係数についての
! 一般化固有値問題 $A a = \sigma B a$ に現れる行列 A, Bを計算する
! 
! ここで a は 基底関数 F_n(y(iy))の係数 a_n (1行 m列), 
! A,Bのj,i要素をca(i,j),cb(i,j)に格納する
! 
! 
subroutine calculate_matrices(m, ty, u, d2u, re, al, ca, cb)

  implicit none

  integer, intent(in) :: m
  real(8), intent(in) :: ty(0:4, 1:m, 0:m-1), u(1:m), d2u(1:m), re, al
  complex(8), intent(out) :: ca(1:m, 1:m), cb(1:m, 1:m)

  integer i, j

  ca(1:m, 1:m) = 0.0d0
  cb(1:m, 1:m) = 0.0d0
  do j=0,m-1
     do i=1,m
        ca(i,j+1) = u(i)*(ty(2,i,j) - al**2*ty(0,i,j)) &
             &      -d2u(i)*ty(0,i,j) &
             &      +dcmplx(0.0d0, 1.0d0)*( ty(4,i,j) - 2.0d0*al**2*ty(2,i,j) + al**4*ty(0,i,j) )/(al*re)
        cb(i,j+1) = ty(2,i,j) - al**2*ty(0,i,j)
     end do
  end do

end subroutine calculate_matrices


! 
! 固有値問題を解き, 第一不安定モードまたは最小減衰モードと, それに対応する
! 固有関数を求める
! 
! Note:
!    非対称行列の一般化固有値問題を解くLapackのサブルーチン zgegv を用いる.
!
!    zgegvの使い方は, 
!    call zgegv( 'N', 'V', N, A, LDA, B, LDB, ALPHA, BETA,
!                VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO )
!
!    LDA, LDBは leading dimension of A, B
!    固有値 lambda は lambda = alpha/beta の形で返る.
!    VLは左固有ベクトル, VRは右固有ベクトル
!    詳しい使い方は http://www.netlib.org/lapack/complex16/zgegv.f を参照すること.
!
subroutine solve_eigenvalue_problem(m, ca, cb, eigenvector, eigenvalue)

  implicit none

  integer, intent(in) :: m
  complex(8), intent(inout) :: ca(1:m, 1:m), cb(1:m, 1:m)
  complex(8), intent(out) :: eigenvector(1:m), eigenvalue

  integer info
  complex(8) alpha(1:m), beta(1:m), dummy(1:1, 1:m), rev(1:m, 1:m), work(1:24*m)
  real(8) work2(1:8*m)
  
  complex(8) lambda(1:m)
  integer i_max

  call zgegv( 'N', 'V', m, ca, m, cb, m, alpha, beta, &
       &      dummy, 1, rev, m, work, 24*m, work2, info )
  if( info /= 0 )then
     eigenvalue = dcmplx(-999.9d0, -999.9d0)
     return
  end if
  lambda(1:m) = alpha(1:m)/beta(1:m)
  
  i_max = maxloc(dimag(lambda), 1)

  eigenvector(1:m) = rev(1:m, i_max)
  eigenvalue = lambda(i_max)

end subroutine solve_eigenvalue_problem
