!--
!----------------------------------------------------------------------
! Copyright(c) 2008-2013 SPMDODEL Development Group. All rights reserved.
!----------------------------------------------------------------------
!
!表題  aq_module
!
!      spml/aq_module モジュールは 1 次元有限領域の下での流体運動を
!      Matsushima and Marcus (1994) で提唱された多項式を用いた
!      スペクトル数値計算ための Fortran90 関数を提供する.
!      このルーチンは離散化にチェビシェフ--ガウス--ラダウ格子点を適用
!      しており, 主に 2 次元極座標,  円筒座標, 球座標の原点の
!      特異性を回避しながらスペクトル計算を行うために用いることを
!      念頭においている.
!
!      2 次元データの 1 次元に関して同時にスペクトル計算を実行するための
!      関数も提供しており, 2, 3 次元領域での計算のベースも提供する.
!
!      Matsushima and Marcus (1994) の多項式に関する説明は
!      doc/spectral_radial.tex を参照のこと.
!
!++
module aq_module
  !
  != aq_module
  !
  ! Authors:: Shin-ichi Takehiro, Youhei SASAKI
  ! Version:: $Id$
  ! Copyright&License:: See COPYRIGHT[link:../COPYRIGHT]
  !
  !== 概要
  !
  ! spml/aq_module モジュールは 1 次元有限領域の下での流体運動を
  ! Matsushima and Marcus (1994) で提唱された多項式を用いたスペクトル法によって
  ! 数値計算するための Fortran90 関数を提供する.
  !
  ! このルーチンは離散化にチェビシェフ--ガウス--ラダウ格子点を適用して
  ! おり, 主に 2 次元極座標, 円筒座標, 球座標の原点の特異性を回避しな
  ! がらスペクトル計算を行うために用いることを念頭においている.
  !
  ! 2 次元データの 1 次元に関して同時にスペクトル計算を実行するための
  ! 関数も提供しており, 2, 3 次元領域での計算のベースも提供する.
  !
  ! Matsushima and Marcus (1994) の多項式に関する説明は
  ! 動径座標のスペクトル法(spectral_radial.pdf[link:spectral_radial.pdf])
  ! を参照のこと.
  !
  !== 関数・変数の名前と型について
  !
  !=== 命名法
  !
  ! * 関数名の先頭 (u_, g_, aq_, ag_) は, 返す値の形を示している.
  !   複数並んだ 2 次元データの第 1 次元の大きさは aq_Initial で
  !   設定する重みの指数の配列 nd の大きさと同じでなければならない.
  !   q_  :: スペクトルデータ
  !   g_  :: 1 次元格子点データ
  !   aq_ :: 1 次元スペクトルデータが複数並んだ 2 次元データ
  !   ag_ :: 1 次元格子点データが複数並んだ 2 次元データ.
  !
  ! * 関数名の間の文字列(Dr)は, その関数の作用を表している.
  !
  ! * 関数名の最後 (_e, _aq, _g, _ag) は, 入力変数の形がスペクトルデータ
  !   および格子点データであることを示している.
  !   _q  :: スペクトルデータ
  !   _g  :: 1 次元格子点データ
  !   _aq :: 1 次元スペクトルデータが複数並んだ 2 次元データ
  !   _ag :: 1 次元格子点データが複数並んだ 2 次元データ
  !
  !=== 各データの種類の説明
  !
  ! * g : 1 次元格子点データ.
  !   * 変数の種類と次元は real(8), dimension(im).
  !   * im は R 座標の格子点数であり, サブルーチン aq_Initial にて
  !     あらかじめ設定しておく.
  !
  ! * q : スペクトルデータ.
  !   * 変数の種類と次元は real(8), dimension(0:km).
  !   * km は R 方向の最大波数であり, サブルーチン aq_Initial にて
  !     あらかじめ設定しておく. スペクトルデータの格納のされ方については...
  !
  ! * ag : 1 次元(R)格子点データの並んだ 2 次元データ.
  !   * 変数の種類と次元は real(8), dimension(size(nd),im).
  !     第 2 次元が R 方向を表す.
  !
  ! * aq : 1 次元スペクトルデータの並んだ 2 次元データ.
  !   * 変数の種類と次元は real(8), dimension(size(nd),0:km).
  !     第 2 次元がスペクトルを表す.
  !
  ! * g_ で始まる関数が返す値は 1 次元格子点データに同じ.
  !
  ! * q_ で始まる関数が返す値はスペクトルデータに同じ.
  !
  ! * ag_ で始まる関数が返す値は 1 次元格子点データの並んだ
  !   2 次元データに同じ.
  !
  ! * aq_ で始まる関数が返す値は 1 次元スペクトルデータの並んだ
  !   2 次元データに同じ.
  !
  ! * スペクトルデータに対する微分等の作用とは, 対応する格子点データに
  !   微分などを作用させたデータをスペクトル変換したもののことである.
  !
  !== 変数・手続き群の要約
  !
  !==== 初期化
  !
  ! aq_Initial  :: スペクトル変換の格子点数, 波数, 領域の大きさの設定
  !
  !==== 座標変数
  !
  ! g_R        :: 格子点座標(R)を格納した 1 次元配列
  ! g_R_Weight :: 重み座標を格納した 1 次元配列
  !
  !==== 基本変換
  !
  ! g_q, ag_aq :: スペクトルデータから格子データへの変換
  ! q_g, aq_ag :: 格子データからスペクトルデータへの変換
  !
  !==== 微分
  !
  ! q_rDr_q, aq_rDr_aq :: スペクトルデータに r(d/dR) 微分を作用させる
  !
  !==== 積
  !
  ! q_r2_q,    aq_r2_aq    :: スペクトルデータに r^2 をかける
  ! q_r2Inv_q, aq_r2Inv_aq :: スペクトルデータに r^-2 をかける
  !
  !==== 積分・平均
  !
  ! a_Int_ag, a_Avr_ag :: 1 次元格子点データの並んだ 2 次元配列の積分および平均
  ! Int_g, Avr_g       :: 1 次元格子点データの積分および平均
  !
  !==== 境界値問題
  !
  ! aq_Boundary_D, aq_Boundary_N         :: 外側ディリクレ条件, ノイマン条件
  ! aq_BoundaryTau_D, aq_BoundaryTau_N   :: 外側ディリクレ条件,
  !                                         ノイマン条件(タウ法)
  ! ag_BoundaryGrid_D, ag_BoundaryGrid_N :: 外側ディリクレ条件,
  !                                         ノイマン条件(選点法)
  !
  use dc_message
  use lumatrix
  implicit none

  private

  public  aq_initial                           ! 初期化
  public  g_R, g_R_Weight                      ! 座標変数

  public  ag_aq, aq_ag, g_q, q_g               ! 基本変換

  public  aq_rDr_aq, q_rDr_q                   ! 微分
  public  aq_r2_aq, q_r2_q                     ! 積
  public  aq_r2Inv_aq, q_r2Inv_q               ! 積

  public  a_Int_ag, a_Avr_ag                   ! 積分・平均
  public  Int_g, Avr_g                         ! 積分・平均

  public  aq_Boundary_D, aq_Boundary_N         ! 境界条件
  public  aq_BoundaryTau_D, aq_BoundaryTau_N   ! 境界条件
  public  ag_BoundaryGrid_D, ag_BoundaryGrid_N ! 境界条件

  interface aq_Boundary_D
     !
     ! 外側境界ディリクレ型境界条件の適用(タウ法).
     !
     ! * 境界条件を適用する配列の次元によって内部でサブルーチンを
     !   使い分けている. ユーザーインターフェースは共通であるので
     !   下部ルーチンを呼ぶ必要はない.
     !
     ! 引数と結果の型
     !
     ! * 1 次元スペクトルデータの並んだ 2 次元配列の場合
     !
     !   real(8), dimension(size(md),0:km),intent(inout)   :: aq_data
     !   !(inout) 境界条件を適用するスペクトルデータ
     !
     !   real(8), dimension(size(md)), intent(in), optional :: value
     !   !(in) 適用する境界値
     !
     ! * 1 次元チェビシェフデータの場合
     !
     !   real(8), dimension(0:km),intent(inout)       :: q_data
     !   !(inout) 境界条件を適用するスペクトルデータ
     !
     !   real(8), intent(in), optional                :: value
     !   !(in) 適用する境界値
     !
     module procedure aq_BoundaryTau_D_1d, aq_BoundaryTau_D_2d
  end interface

  interface aq_Boundary_N
     !
     ! 外側境界ノイマン型境界条件の適用(タウ法).
     !
     ! * 境界条件を適用する配列の次元によって内部でサブルーチンを
     !   使い分けている. ユーザーインターフェースは共通であるので
     !   下部ルーチンを呼ぶ必要はない.
     !
     ! 引数と結果の型
     !
     ! * 1 次元スペクトルデータの並んだ 2 次元配列の場合
     !
     !   real(8), dimension(size(md),0:km),intent(inout)    :: aq_data
     !   !(inout) 境界条件を適用するスペクトルデータ
     !
     !   real(8), dimension(size(md)), intent(in), optional :: value
     !   !(in) 適用する境界値
     !
     ! * 1 次元チェビシェフデータの場合
     !
     !   real(8), dimension(0:km),intent(inout)       :: q_data
     !   !(inout) 境界条件を適用するスペクトルデータ
     !
     !   real(8), intent(in), optional                :: value
     !   !(in) 適用する境界値
     !
     module procedure aq_BoundaryTau_N_1d, aq_BoundaryTau_N_2d
  end interface

  interface aq_BoundaryTau_D
     !
     ! 外側境界ディリクレ型境界条件の適用(タウ法).
     !
     ! * 境界条件を適用する配列の次元によって内部でサブルーチンを
     !   使い分けている. ユーザーインターフェースは共通であるので
     !   下部ルーチンを呼ぶ必要はない.
     !
     ! 引数と結果の型
     !
     ! * 1 次元スペクトルデータの並んだ 2 次元配列の場合
     !
     !   real(8), dimension(size(md),0:km),intent(inout)    :: aq_data
     !   !(inout) 境界条件を適用するチェビシェフデータ
     !
     !   real(8), dimension(size(md)), intent(in), optional :: value
     !   !(in) 適用する境界値
     !
     ! * 1 次元スペクトルデータの場合
     !
     !   real(8), dimension(0:km),intent(inout)       :: q_data
     !   !(inout) 境界条件を適用するチェビシェフデータ
     !
     !   real(8), intent(in), optional                :: value
     !   !(in) 適用する境界値
     !
     module procedure aq_BoundaryTau_D_1d, aq_BoundaryTau_D_2d
  end interface

  interface aq_BoundaryTau_N
     !
     ! 外側境界ノイマン型境界条件の適用(タウ法).
     !
     ! * 境界条件を適用する配列の次元によって内部でサブルーチンを
     !   使い分けている. ユーザーインターフェースは共通であるので
     !   下部ルーチンを呼ぶ必要はない.
     !
     ! 引数と結果の型
     !
     ! * 1 次元スペクトルデータの並んだ 2 次元配列の場合
     !
     !   real(8), dimension(size(md),0:km),intent(inout)    :: aq_data
     !   !(inout) 境界条件を適用するチェビシェフデータ
     !
     !   real(8), dimension(size(md)), intent(in), optional :: value
     !   !(in) 適用する境界値
     !
     ! * 1 次元スペクトルデータの場合
     !
     !   real(8), dimension(0:km),intent(inout)       :: q_data
     !   !(inout) 境界条件を適用するチェビシェフデータ
     !
     !   real(8), intent(in), optional                :: value
     !   !(in) 適用する境界値
     !
     module procedure aq_BoundaryTau_N_1d, aq_BoundaryTau_N_2d
  end interface

  interface ag_BoundaryGrid_D
     !
     ! 外側境界ディリクレ型境界条件の適用(選点法).
     !
     ! * 境界条件を適用する配列の次元によって内部でサブルーチンを
     !   使い分けている. ユーザーインターフェースは共通であるので
     !   下部ルーチンを呼ぶ必要はない.
     !
     ! 引数と結果の型
     !
     ! * 1 次元スペクトルデータの並んだ 2 次元配列の場合
     !
     !   real(8), dimension(size(md),0:km),intent(inout)    :: aq_data
     !   !(inout) 境界条件を適用するチェビシェフデータ
     !
     !   real(8), dimension(size(md)), intent(in), optional :: value
     !   !(in) 適用する境界値
     !
     ! * 1 次元スペクトルデータの場合
     !
     !   real(8), dimension(0:km),intent(inout)       :: q_data
     !   !(inout) 境界条件を適用するチェビシェフデータ
     !
     !   real(8), intent(in), optional                :: value
     !   !(in) 適用する境界値
     !
     module procedure ag_BoundaryGrid_D_1d, ag_BoundaryGrid_D_2d
  end interface

  interface ag_BoundaryGrid_N
     !
     ! 外側境界ノイマン型境界条件の適用(選点法).
     !
     ! * 境界条件を適用する配列の次元によって内部でサブルーチンを
     !   使い分けている. ユーザーインターフェースは共通であるので
     !   下部ルーチンを呼ぶ必要はない.
     !
     ! 引数と結果の型
     !
     ! * 1 次元スペクトルデータの並んだ 2 次元配列の場合
     !
     !   real(8), dimension(size(md),0:km),intent(inout)    :: aq_data
     !   !(inout) 境界条件を適用するチェビシェフデータ
     !
     !   real(8), dimension(size(md)), intent(in), optional :: value
     !   !(in) 適用する境界値
     !
     ! * 1 次元スペクトルデータの場合
     !
     !   real(8), dimension(0:km),intent(inout)       :: q_data
     !   !(inout) 境界条件を適用するチェビシェフデータ
     !
     !   real(8), intent(in), optional                :: value
     !   !(in) 適用する境界値
     !
     module procedure ag_BoundaryGrid_N_1d, ag_BoundaryGrid_N_2d
  end interface

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

  real(8) :: alpha                        ! 展開多項式パラメター  0 < α <= 1
  real(8) :: beta                         ! 展開多項式パラメター  0 < β
  real(8) :: gamma                        ! 展開多項式パラメター  γ=2α+β
  integer, allocatable :: md(:)           ! 展開多項式上付次数のならび

  integer :: jmax                         ! a 座標(データ第 1 次元)の大きさ

  real(8), allocatable :: g_R(:)          ! ガウス--ラダウ格子点
  real(8), allocatable :: g_R_Weight(:)   ! ガウス重み

  real(8), allocatable :: CF(:,:,:)       ! 正変換用行列
  real(8), allocatable :: CB(:,:,:)       ! 逆変換用行列


  logical :: first_r2inv  =.true.
  logical :: first_Tau_D  =.true.
  logical :: first_Tau_N  =.true.
  logical :: first_Grid_N =.true.

  save im, km, ra, alpha, beta, gamma, CF, CB, g_R, g_R_Weight, md, jmax
  save first_r2inv, first_Tau_D, first_Tau_N, first_Grid_N

contains

! --- 初期化
  subroutine aq_Initial(i_in,k_in,r_in,alpha_in,beta_in,md_in)
    !
    ! スペクトル変換の格子点数, 波数, 領域の大きさ, 重みを設定する.
    !
    ! 他の関数や変数を呼ぶ前に, 最初にこのサブルーチンを呼んで
    ! 初期設定をしなければならない.
    !
    integer,intent(in) :: i_in              !(in) 格子点数
    integer,intent(in) :: k_in              !(in) 切断波数
    real(8),intent(in) :: r_in              !(in) 外側境界の座標(半径)
    real(8),intent(in) :: alpha_in          !(in) 展開多項式パラメター
    real(8),intent(in) :: beta_in           !(in) 展開多項式パラメター
    integer,intent(in) :: md_in(:)          !(in) 展開多項式上付次数のならび

    integer i, j, n

   !--- パラメターのチェック・記憶
    im=i_in ; km=k_in ; ra = r_in; alpha=alpha_in; beta=beta_in

    if ( km .ge. 2*im ) then
       call MessageNotify('E','aq_initial','KM shoud be less than 2*IM')
    endif

    if ( alpha .le. 0.0D0 ) then
       call MessageNotify('E','aq_initial','alpha must be larger than 0')
    endif
    if ( alpha .gt. 1.0D0 ) then
       call MessageNotify('E','aq_initial','alpha must be smaller equal to 1')
    endif
    if ( beta .le. 0.0D0 ) then
       call MessageNotify('E','aq_initial','beta must be larger than 0')
    endif

    gamma = 2.0D0 * alpha + beta

   !--- 次数情報の計算
    jmax = size(md_in)
    if ( allocated(md) ) deallocate(md)
    allocate(md(jmax))

    md = md_in

   !--- 格子点, 重みの設定
    if ( allocated(g_R) ) deallocate(g_R)
    if ( allocated(g_R_Weight) ) deallocate(g_R_Weight)
    allocate(g_R(im),g_R_Weight(im))

    call gauss_radau(2*im, g_R, g_R_Weight)
    g_R = ra * g_R
    g_R_Weight = ra**(gamma-1) * g_R_Weight

   !--- 変換行列の設定
    if ( allocated(CF) ) deallocate(CF)
    if ( allocated(CB) ) deallocate(CB)
    allocate(CF(jmax,0:km,im),CB(jmax,im,0:km))

    CF = 0.0D0
    do j=1,jmax
       do n=md(j),km,2
          do i=1,im
             CF(j,n,i) = Phi(g_R(i)/ra,n,md(j)) * g_R_Weight(i)/ra**(gamma-1)
          enddo
       enddo
    enddo

    CB = 0.0D0
    do j=1,jmax
       do i=1,im
          do n=md(j),km,2
             CB(j,i,n) = Phi(g_R(i)/ra,n,md(j))
          enddo
       enddo
    enddo

   !--- 各ルーチン変換行列の初期化スイッチ

    first_r2inv  = .true.
    first_Tau_D  = .true.
    first_Tau_N  = .true.
    first_Grid_N = .true.

    call MessageNotify(&
      'M','aq_initial','aq_module (2009/07/31) is initialized')

  end subroutine aq_Initial

! ---- 逆変換 ----
  function ag_aq(aq_data)
    !
    ! スペクトルデータから格子データへ変換する(2 次元配列用).
    !
    real(8), dimension(:,0:), intent(in)    :: aq_data
    !(in) スペクトルデータ

    real(8), dimension(size(aq_data,1),im)  :: ag_aq
    !(out) 格子点データ

    integer :: i, j

    if ( size(aq_data,1) /= jmax ) then
       call MessageNotify('E','ag_aq', &
            '1st dim. of the spectral data should be same as dim. of MD.')
    end if

    if ( size(aq_data,2)-1 < km ) then
       call MessageNotify('E','ag_aq', &
            'The spectral dimension of input data too small.')
    elseif ( size(aq_data,2)-1 > km ) then
       call MessageNotify('W','ag_aq', &
            'The spectral dimension of input data too large.')
    endif

    ag_aq = 0.0D0
    do i=1,im
!$omp parallel do
       do j=1,jmax
             ag_aq(j,i) = sum(CB(j,i,md(j):km:2)*aq_data(j,md(j):km:2))
       enddo
!$omp end parallel do
    enddo

  end function ag_aq

  function g_q(q_data)
    !
    ! スペクトルデータから格子データへ変換する(1 次元配列用).
    !
    real(8), dimension(:), intent(in)  :: q_data
    !(in) スペクトルデータ

    real(8), dimension(im)             :: g_q
    !(out) 格子点データ

    real(8), dimension(1,size(q_data)) :: q_work
    ! 作業用配列
    real(8), dimension(1,im)         :: g_work
    ! 作業用配列

    q_work(1,:) = q_data
    g_work = ag_aq(q_work)
    g_q = g_work(1,:)

  end function g_q

! ---- 正変換 ----
  function aq_ag(ag_data)
    !
    ! 格子データからスペクトルデータへ変換する(2 次元配列用).
    !
    real(8), dimension(:,:), intent(in)      :: ag_data
    !(in) 格子点データ

    real(8), dimension(size(ag_data,1),0:km) :: aq_ag
    !(out) スペクトルデータ

    integer :: j, n

    if ( size(ag_data,1) /= jmax ) then
       call MessageNotify('E','aq_ag', &
            '1st dim. of the grid data should be same as dim. of MD.')
    end if

    if ( size(ag_data,2) < im ) then
       call MessageNotify('E','aq_ag', &
            'The Grid points of input data too small.')
    elseif ( size(ag_data,2) > im ) then
       call MessageNotify('W','aq_ag', &
            'The Grid points of input data too large.')
    endif

    aq_ag = 0.0D0
!$omp parallel do private(n)
    do j=1,jmax
       do n=md(j),km,2
          aq_ag(j,n) = sum(CF(j,n,:)*ag_data(j,:))
       enddo
    enddo
!$omp end parallel do

  end function aq_ag

  function q_g(g_data)
    !
    ! 格子データからスペクトルデータへ変換する(1 次元配列用).
    !
    real(8), dimension(:), intent(in)     :: g_data
    !(in) 格子点データ

    real(8), dimension(0:km)              :: q_g
    !(out) スペクトルデータ

    real(8), dimension(1,size(g_data)) :: ag_work
    real(8), dimension(1,0:km)         :: aq_work

    ag_work(1,:) = g_data
    aq_work = aq_ag(ag_work)
    q_g = aq_work(1,:)

  end function q_g

! ---- 微分計算 ----
  function aq_rDr_aq(aq_data)
    !
    ! 入力スペクトルデータに対して微分 r(d/dr) のスペクトル係数
    ! を計算する(2 次元配列用).
    !
    !   a_n = aq_data/sqrt(Inm),  b_n = aq_rDr_aq/sqrt(Inm)
    !
    !   b_n = (2n+gamma-1)/(2n+gamma+3)b_n+2
    !       + (2n+gamma-1)(n+gamma+1)/(2n+gamma+3)a_n+2 + n a_n
    !
    real(8), dimension(:,0:), intent(in)     :: aq_data
    !(in) 入力スペクトルデータ

    real(8), dimension(size(aq_data,1),0:size(aq_data,2)-1) :: aq_rDr_aq
    !(out) 出力スペクトルデータ

    integer :: j, n
    integer :: nstr, nend
    real(8) :: sqInm, sqInm2

    if ( size(aq_data,1) /= jmax ) then
       call MessageNotify('E','aq_rDr_aq', &
            '1st dim. of the spectral data should be same as dim. of MD.')
    end if

    if ( size(aq_data,2)-1 < km ) then
       call MessageNotify('E','aq_rDr_aq', &
            'The spectral dimension of input data too small.')
    elseif ( size(aq_data,2)-1 > km ) then
       call MessageNotify('W','aq_rDr_aq', &
            'The spectral dimension of input data too large.')
    endif

    aq_rDr_aq = 0.0D0
!$omp parallel do private(nstr,nend,sqInm,sqInm2,n)
    do j=1,jmax
       if ( mod(md(j),2) .eq. mod(km,2) ) then
          nstr=md(j) ; nend=km
       else
          nstr=md(j) ; nend=km-1
       endif

       sqInm = sqrt(Inm(nend,md(j)))
       aq_rDr_aq(j,nend) = nend * aq_data(j,nend)

       do n=nend-2,nstr,-2
          sqInm2 = sqInm
          sqInm = sqrt(Inm(n,md(j)))

          aq_rDr_aq(j,n) = sqInm/sqInm2 * (&
                 (2*n+gamma-1)/(2*n+gamma+3) * aq_rDr_aq(j,n+2) &
               + (2*n+gamma-1)*(n+gamma+1)/(2*n+gamma+3) * aq_data(j,n+2) )&
               + n * aq_data(j,n)
       enddo
    enddo
!$omp end parallel do

  end function aq_rDr_aq

  function q_rDr_q(q_data)
    !
    ! 入力スペクトルデータに r(d/dR) 微分を作用する(1 次元配列用).
    !
    ! スペクトルデータの r(d/dR) 微分とは, 対応する格子点データに R 微分を
    ! 作用させたデータのスペクトル変換のことである.
    !
    !
    real(8), dimension(:), intent(in)   :: q_data
    !(in) 入力チェビシェフデータ

    real(8), dimension(0:km)            :: q_rDr_q
    !(out) チェビシェフデータの R 微分

    real(8), dimension(1,size(q_data))  :: aq_work
    ! 作業用配列

    aq_work(1,:) = q_data
    aq_work = aq_rDr_aq(aq_work)
    q_rDr_q = aq_work(1,:)

  end function q_rDr_q

! ---- r^2 積計算 ----
  function aq_r2_aq(aq_data)
    !
    ! 入力スペクトルデータに対して積 r^2 のスペクトル係数
    ! を計算する(2 次元配列用).
    !
    !   a_n^m = aq_data/sqrt(Inm),  b_n^m = aq_rDr_aq/sqrt(Inm)
    !
    !   b_n^m = (n-|m|)(n+|m|+beta-1)/((2n+gamma-5)(2n+gamma-3)) a_n-2^m
    !       + (2n(n+gamma-1) + 2|m|(|m|+beta-1)+(gamma-3)(beta+1)
    !            /((2n+gamma+1)(2n+gamma-3)) a_n^m
    !       + (n-|m|+gamma-beta)(n+|m|+gamma-1)
    !            /((2n+gamma+3)(2n+gamma+1)) a_n+2^m
    !
    real(8), dimension(:,0:), intent(in)     :: aq_data
    !(in) 入力スペクトルデータ

    real(8), dimension(size(aq_data,1),0:size(aq_data,2)-1) :: aq_r2_aq
    !(out) 出力スペクトルデータ

    integer :: j, n, m
    integer :: nstr, nend
    real(8) :: sqrInp2m, sqrInm, sqrInm2m

    if ( size(aq_data,1) /= jmax ) then
       call MessageNotify('E','aq_r2_aq', &
            '1st dim. of the spectral data should be same as dim. of MD.')
    end if

    if ( size(aq_data,2)-1 < km ) then
       call MessageNotify('E','aq_r2_aq', &
            'The spectral dimension of input data too small.')
    elseif ( size(aq_data,2)-1 > km ) then
       call MessageNotify('W','aq_r2_aq', &
            'The spectral dimension of input data too large.')
    endif

    aq_r2_aq = 0.0D0

!$omp parallel do private(nstr,nend,sqrInm,sqrInm2m,sqrInp2m,n,m)
    do j=1,jmax
       if ( mod(md(j),2) .eq. mod(km,2) ) then
          nstr=md(j) ; nend=km
       else
          nstr=md(j) ; nend=km-1
       endif

       m = abs(md(j))

       n = nstr
       sqrInm   = sqrt(Inm(n,m))
       sqrInp2m = sqrt(Inm(n+2,m))
       aq_r2_aq(j,n) =                                            &
            (2*n*(n+gamma-1) + 2*m*(m+beta-1)+(gamma-3)*(beta+1)) &
               /((2*n+gamma+1)*(2*n+gamma-3)) * aq_data(j,n)      &
          + (n-m+gamma-beta)*(n+m+gamma-1)                        &
               /((2*n+gamma+3)*(2*n+gamma+1))                     &
                * aq_data(j,n+2) * (sqrInm/sqrInp2m)

       do n=nstr+2,nend-2,2
          sqrInm2m = sqrInm
          sqrInm   = sqrInp2m
          sqrInp2m = sqrt(Inm(n+2,m))

          aq_r2_aq(j,n) = &
             (n-m)*(n+m+beta-1)/((2*n+gamma-5)*(2*n+gamma-3))      &
             * aq_data(j,n-2)* (sqrInm/sqrInm2m)                   &
           + (2*n*(n+gamma-1) + 2*m*(m+beta-1)+(gamma-3)*(beta+1)) &
                /((2*n+gamma+1)*(2*n+gamma-3)) * aq_data(j,n)      &
           + (n-m+gamma-beta)*(n+m+gamma-1)                        &
                /((2*n+gamma+3)*(2*n+gamma+1))                     &
                * aq_data(j,n+2) * (sqrInm/sqrInp2m)
       enddo

       n = nend
       sqrInm2m = sqrInm
       sqrInm   = sqrInp2m
       aq_r2_aq(j,n) = &
            (n-m)*(n+m+beta-1)/((2*n+gamma-5)*(2*n+gamma-3))       &
             * aq_data(j,n-2)* (sqrInm/sqrInm2m)                   &
           + (2*n*(n+gamma-1) + 2*m*(m+beta-1)+(gamma-3)*(beta+1)) &
                /((2*n+gamma+1)*(2*n+gamma-3)) * aq_data(j,n)
    enddo
!$omp end parallel do

    aq_r2_aq = aq_r2_aq * ra**2

  end function aq_r2_aq

  function q_r2_q(q_data)
    !
    ! 入力スペクトルデータに対して積 r^2 のスペクトル係数
    ! を計算する(1 次元配列用).
    !
    real(8), dimension(:), intent(in)   :: q_data
    !(in) 入力チェビシェフデータ

    real(8), dimension(0:km)            :: q_r2_q
    !(out) チェビシェフデータの R 微分

    real(8), dimension(1,size(q_data))  :: aq_work
    ! 作業用配列

    aq_work(1,:) = q_data
    aq_work = aq_r2_aq(aq_work)
    q_r2_q = aq_work(1,:)

  end function q_r2_q

! ---- r^-2 積計算 ----
  function aq_r2Inv_aq(aq_data)
    !
    ! 入力スペクトルデータに対して積 r^-2 のスペクトル係数
    ! を計算する(2 次元配列用).
    !
    !   a_n^m = aq_r2Inv_aq/sqrt(Inm),  b_n^m = aq_data/sqrt(Inm)
    !
    !   b_n^m = (n-|m|)(n+|m|+beta-1)/((2n+gamma-5)(2n+gamma-3)) a_n-2^m
    !       + (2n(n+gamma-1) + 2|m|(|m|+beta-1)+(gamma-3)(beta+1)
    !            /((2n+gamma+1)(2n+gamma-3)) a_n^m
    !       + (n-|m|+gamma-beta)(n+|m|+gamma-1)
    !            /((2n+gamma+3)(2n+gamma+1)) a_n+2^m
    !
    real(8), dimension(:,0:), intent(in)     :: aq_data
    !(in) 入力スペクトルデータ

    real(8), dimension(size(aq_data,1),0:size(aq_data,2)-1) :: aq_r2Inv_aq
    !(out) 出力スペクトルデータ

    real(8), dimension(:,:,:), allocatable    :: R2MTX
    integer, dimension(:,:), allocatable      :: kp

    integer :: j, n, m
    integer :: nstr, nend
    real(8) :: sqrInp2m, sqrInm, sqrInm2m

    save R2MTX, kp

    if ( size(aq_data,1) /= jmax ) then
       call MessageNotify('E','aq_r2Inv_aq', &
            '1st dim. of the spectral data should be same as dim. of MD.')
    end if

    if ( size(aq_data,2)-1 < km ) then
       call MessageNotify('E','aq_r2Inv_aq', &
            'The spectral dimension of input data too small.')
    elseif ( size(aq_data,2)-1 > km ) then
       call MessageNotify('W','aq_r2Inv_aq', &
            'The spectral dimension of input data too large.')
    endif

    if ( first_r2inv ) then
       first_r2inv = .false.
       if ( allocated(R2MTX) ) deallocate(R2MTX)
       if ( allocated(kp) ) deallocate(kp)
       allocate(R2MTX(jmax,0:km,0:km),kp(jmax,0:km))

       R2MTX(:,:,:) = 0.0D0
       do n=0,km
          R2MTX(:,n,n) = 1.0D0
       enddo

!$omp parallel do private(nstr,nend,sqrInm,sqrInm2m,sqrInp2m,n,m)
       do j=1,jmax
          if ( mod(md(j),2) .eq. mod(km,2) ) then
             nstr=md(j) ; nend=km
          else
             nstr=md(j) ; nend=km-1
          endif

          m = abs(md(j))

          n = nstr

          sqrInm   = sqrt(Inm(n,m))
          sqrInp2m = sqrt(Inm(n+2,m))

          R2MTX(j,n,n)   = &
               (2*n*(n+gamma-1) + 2*m*(m+beta-1)+(gamma-3)*(beta+1)) &
               /((2*n+gamma+1)*(2*n+gamma-3))

          R2MTX(j,n,n+2) = &
               (n-m+gamma-beta)*(n+m+gamma-1)                        &
               /((2*n+gamma+3)*(2*n+gamma+1))                        &
               * (sqrInm/sqrInp2m)

          do n=nstr+2,nend-2,2
             sqrInm2m = sqrInm
             sqrInm   = sqrInp2m
             sqrInp2m = sqrt(Inm(n+2,m))

             R2MTX(j,n,n-2) = &
                (n-m)*(n+m+beta-1)/((2*n+gamma-5)*(2*n+gamma-3))      &
                *(sqrInm/sqrInm2m)

             R2MTX(j,n,n)   = &
                (2*n*(n+gamma-1) + 2*m*(m+beta-1)+(gamma-3)*(beta+1)) &
                   /((2*n+gamma+1)*(2*n+gamma-3))

             R2MTX(j,n,n+2) = &
                (n-m+gamma-beta)*(n+m+gamma-1)                        &
                   /((2*n+gamma+3)*(2*n+gamma+1))                     &
                   * (sqrInm/sqrInp2m)
          enddo

          n = nend
          sqrInm2m = sqrInm
          sqrInm   = sqrInp2m
             R2MTX(j,n,n-2) = &
                (n-m)*(n+m+beta-1)/((2*n+gamma-5)*(2*n+gamma-3))      &
                *(sqrInm/sqrInm2m)
             R2MTX(j,n,n)   = &
                (2*n*(n+gamma-1) + 2*m*(m+beta-1)+(gamma-3)*(beta+1)) &
                   /((2*n+gamma+1)*(2*n+gamma-3))
       enddo
!$omp end parallel do

       R2MTX = R2MTX * ra**2

       call LuDecomp(R2MTX,kp)

    end if

    aq_r2Inv_aq = LuSolve(R2MTX,kp,aq_data)

  end function aq_r2Inv_aq

  function q_r2Inv_q(q_data)
    !
    ! 入力スペクトルデータに対して積 r^2 のスペクトル係数
    ! を計算する(1 次元配列用).
    !
    real(8), dimension(:), intent(in)   :: q_data
    !(in) 入力チェビシェフデータ

    real(8), dimension(0:km)            :: q_r2Inv_q
    !(out) チェビシェフデータの R 微分

    real(8), dimension(1,size(q_data))  :: aq_work
    ! 作業用配列

    aq_work(1,:) = q_data
    aq_work = aq_r2Inv_aq(aq_work)
    q_r2Inv_q = aq_work(1,:)

  end function q_r2Inv_q

!--------------- 積分計算 -----------------
  function a_Int_ag(ag)
    !
    ! 1 次元格子点データが並んだ 2 次元配列の積分
    !
    !      \int_0^a f(R) W(R) dR, W(R) = R^beta/(a^2-R^2)^(1-alpha)
    !
    ! を計算する.
    !
    real(8), dimension(:,:), intent(in)      :: ag
    !(in)入力格子点データ

    real(8), dimension(size(ag,1))           :: a_Int_ag
    !(out) 積分したデータ
    integer :: i

    if ( size(ag,2) < im ) then
       call MessageNotify('E','a_Int_ag', &
            'The Grid points of input data too small.')
    elseif ( size(ag,2) > im ) then
       call MessageNotify('W','a_Int_ag', &
            'The Grid points of input data too large.')
    endif

    a_Int_ag = 0.0d0
    do i=1,im
       a_Int_ag(:) = a_Int_ag(:) + ag(:,i)*g_R_Weight(i)
    enddo
  end function a_Int_ag

  function Int_g(g)
    !
    ! 1 次元格子点データの積分
    !
    !      \int_0^a f(R) W(R) dR, W(R) = R^beta/(a^2-R^2)^(1-alpha)
    !
    ! を行う
    !
    real(8), dimension(im), intent(in)   :: g
    !(in) 格子点データ

    real(8)                              :: Int_g
    !(out) 積分値

    Int_g = sum(g*g_R_Weight)

  end function Int_g

  function a_Avr_ag(ag)
    !
    ! 1 次元格子点データが並んだ 2 次元配列の平均
    !
    !      \int_0^a f(R) W(R) dR/\int_0^a W(R) dR,
    !          W(R) = R^beta/(a^2-R^2)^(1-alpha)
    !
    ! を計算する.
    !
    real(8), dimension(:,0:), intent(in)   :: ag
    !(in)入力格子点データ

    real(8), dimension(size(ag,1))         :: a_Avr_ag
    !(out) 平均したデータ

    a_Avr_ag = a_Int_ag(ag)/sum(g_R_Weight)

  end function a_Avr_ag

  function Avr_g(g)
    !
    ! 1 次元格子点データの平均
    !
    !      \int_0^a f(R) W(R) dR/\int_0^a W(R) dR,
    !          W(R) = R^beta/(a^2-R^2)^(1-alpha)
    !
    ! を計算する.
    !
    real(8), dimension(im), intent(in)   :: g
    !(in) 格子点データ

    real(8)                              :: Avr_g
    !(out) 積分値

    Avr_g = Int_g(g)/sum(g_R_Weight)
  end function Avr_g

!---- Dirichlet 型境界条件(タウ法) ----

  subroutine aq_BoundaryTau_D_2d(aq_data,value)
    !
    ! Dirichlet 型境界条件の適用(タウ法, 2 次元配列用)
    ! * 外側境界(i=im)での値を与える.
    !
    real(8), dimension(:,0:),intent(inout)         :: aq_data
    !(inout) 境界条件を適用するチェビシェフデータ(jmax,0:km)

    real(8), dimension(:), intent(in), optional  :: value
    !(in) 境界値(jmax)

    real(8), dimension(:,:,:), allocatable  :: alu
    integer, dimension(:,:), allocatable    :: kp
    real(8), dimension(size(aq_data,1),0:km):: aq_work
    real(8), dimension(size(aq_data,1),im)  :: ag_work
    real(8), dimension(size(aq_data,1))     :: value0           ! 境界値

    integer :: k, j, kstr, kend
    save    :: alu, kp

    if ( size(aq_data,2)-1 < km ) then
       call MessageNotify('E','aq_BoundaryTau_D', &
            'The spectral dimension of input data too small.')
    elseif ( size(aq_data,2)-1 > km ) then
       call MessageNotify('W','aq_BoundaryTau_D', &
            'The spectral dimension of input data too large.')
    endif

    if (.not. present(value)) then
       value0=0.0D0
    else
       value0 = value
    endif

    if ( first_Tau_D ) then
       first_Tau_D = .false.

       allocate(alu(size(aq_data,1),0:km,0:km),kp(size(aq_data,1),0:km))

       alu=0.0D0
       do k=0,km
          alu(:,k,k) = 1.0D0
       enddo

       do j=1,jmax
          if ( mod(md(j),2) .eq. mod(km,2) ) then
             kstr = md(j) ; kend = km
          else
             kstr = md(j) ; kend = km-1
          endif

          do k=kstr,kend,2
             aq_work = 0.0D0
             aq_work(j,k) = 1.0D0
             ag_work = ag_aq(aq_work)
             alu(j,kend,k) = ag_work(j,im)
          enddo
       enddo

       call ludecomp(alu,kp)
    endif

    do j=1,jmax
       if ( mod(md(j),2) .eq. mod(km,2) ) then
          aq_data(j,km)   = value0(j)
       else
          aq_data(j,km-1)   = value0(j)
       endif
    enddo

    aq_data = lusolve(alu,kp,aq_data)

  end subroutine aq_BoundaryTau_D_2d

  subroutine aq_BoundaryTau_D_1d(q_data,value)
    !
    ! Dirichlet 型境界条件の適用(タウ法, 1 次元配列用)
    ! * 両境界での値を与える.
    !
    real(8), dimension(0:km),intent(inout)       :: q_data
    !(inout) 境界条件を適用するチェビシェフデータ(0:km)

    real(8), intent(in), optional                :: value
    !(in) 境界値

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

    if (.not. present(value)) then
       vwork(1)=0.0D0
    else
       vwork(1) = value
    endif

    aq_work(1,:)=q_data
    call aq_BoundaryTau_D_2d(aq_work,vwork)
    q_data=aq_work(1,:)

  end subroutine aq_BoundaryTau_D_1d

!---- Neumann 型境界条件(タウ法) ----

  subroutine aq_BoundaryTau_N_2d(aq_data,value)
    !
    ! 外側境界 Neumann 型境界条件の適用(タウ法, 2 次元配列用)
    ! * i=im で勾配の値を与える.
    !
    real(8), dimension(:,0:),intent(inout)         :: aq_data
    !(inout) 境界条件を適用するチェビシェフデータ(m,0:km)

    real(8), dimension(:), intent(in), optional    :: value
    !(in) 境界値(m)

    real(8), dimension(:,:,:), allocatable    :: alu
    integer, dimension(:,:), allocatable      :: kp
    real(8), dimension(size(aq_data,1),0:km)  :: aq_work
    real(8), dimension(size(aq_data,1),im)    :: ag_work
    real(8), dimension(size(aq_data,1))       :: value0           ! 境界値

    integer :: j, k, kstr, kend

    save    :: alu, kp

    if ( size(aq_data,2)-1 < km ) then
       call MessageNotify('E','aq_BoundaryTau_N', &
            'The spectral dimension of input data too small.')
    elseif ( size(aq_data,2)-1 > km ) then
       call MessageNotify('W','aq_BoundaryTau_N', &
            'The spectral dimension of input data too large.')
    endif

    if (.not. present(value)) then
       value0=0.0D0
    else
       value0 = value
    endif

    if ( first_Tau_N ) then
       first_Tau_N = .false.
       allocate(alu(jmax,0:km,0:km),kp(jmax,0:km))

       alu=0.0D0
       do k=0,km
          alu(:,k,k) = 1.0D0
       enddo

       do j=1,jmax
          if ( mod(md(j),2) .eq. mod(km,2) ) then
             kstr=md(j) ; kend = km
          else
             kstr=md(j) ; kend = km-1
          endif

          do k=kstr,kend,2
             aq_work = 0.0D0
             aq_work(j,k) = 1.0D0
             ag_work = ag_aq(aq_rDr_aq(aq_work))/spread(g_R,1,jmax)
             alu(j,kend,k) = ag_work(j,im)
          enddo
       enddo

       call ludecomp(alu,kp)
    endif

    do j=1,jmax
       if ( mod(md(j),2) .eq. mod(km,2) ) then
          aq_data(j,km)   = value0(j)
       else
          aq_data(j,km-1)   = value0(j)
       endif
    enddo

    aq_data = lusolve(alu,kp,aq_data)

  end subroutine aq_BoundaryTau_N_2d

  subroutine aq_BoundaryTau_N_1d(q_data,value)
    !
    ! Dirichlet/Neumann 型境界条件の適用(タウ法, 1 次元配列用)
    ! * i=0 で勾配の値を与える.
    !
    real(8), dimension(0:km),intent(inout)       :: q_data
    !(inout) 境界条件を適用するチェビシェフデータ(0:km)

    real(8), intent(in), optional                :: value
    !(in) 境界値

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

    if (.not. present(value)) then
       vwork(1)=0.0D0
    else
       vwork(1) = value
    endif

    aq_work(1,:)=q_data
    call aq_BoundaryTau_N_2d(aq_work,vwork)
    q_data=aq_work(1,:)

  end subroutine aq_BoundaryTau_N_1d

!---- Dirichlet 型境界条件(選点法) ----

  subroutine ag_BoundaryGrid_D_2d(ag_data,value)
    !
    ! Dirichlet 型境界条件の適用(選点法, 2 次元配列用)
    ! * 外側境界(i=im)での値を与える.
    !
    real(8), dimension(:,:),intent(inout)         :: ag_data
    !(inout) 境界条件を適用するチェビシェフデータ(jmax,im)

    real(8), dimension(:), intent(in), optional  :: value
    !(in) 境界値(jmax)

    real(8), dimension(size(ag_data,1))       :: value0

    if (.not. present(value)) then
       value0=0.0d0
    else
       value0 = value
    endif
    ag_data(:,im) = value0

  end subroutine ag_BoundaryGrid_D_2d

  subroutine ag_BoundaryGrid_D_1d(g_data,value)
    !
    ! Dirichlet 型境界条件の適用(タウ法, 1 次元配列用)
    ! * 両境界での値を与える.
    !
    real(8), dimension(im),intent(inout)       :: g_data
    !(inout) 境界条件を適用するチェビシェフデータ(im)

    real(8), intent(in), optional              :: value
    !(in) 境界値

    real(8), dimension(1,im)                   :: ag_work
    real(8), dimension(1)                      :: vwork

    if (.not. present(value)) then
       vwork(1)=0.0d0
    else
       vwork(1) = value
    endif

    ag_work(1,:)=g_data
    call ag_BoundaryGrid_D_2d(ag_work,vwork)
    g_data=ag_work(1,:)

  end subroutine ag_BoundaryGrid_D_1d

!---- Neumann 型境界条件(選点法) ----

  subroutine ag_BoundaryGrid_N_2d(ag_data,value)
    !
    ! 外側境界 Neumann 型境界条件の適用(選点法, 2 次元配列用)
    ! * i=im で勾配の値を与える.
    !
    real(8), dimension(:,:),intent(inout)         :: ag_data
    !(inout) 境界条件を適用するチェビシェフデータ(m,im)

    real(8), dimension(:), intent(in), optional    :: value
    !(in) 境界値(m)

    real(8), dimension(:,:,:), allocatable    :: alu
    integer, dimension(:,:), allocatable      :: kp
    real(8), dimension(size(ag_data,1),im)    :: ag_work
    real(8), dimension(size(ag_data,1))       :: value0           ! 境界値

    integer :: i
    save    :: alu, kp

    if ( size(ag_data,2) < im ) then
       call MessageNotify('E','aq_BoundaryGrid_N', &
            'The dimension of input data too small.')
    elseif ( size(ag_data,2) > im ) then
       call MessageNotify('W','aq_BoundaryGrid_N', &
            'The dimension of input data too large.')
    endif

    if (.not. present(value)) then
       value0=0.0D0
    else
       value0 = value
    endif

    if ( first_Grid_N ) then
       first_Grid_N = .false.
       allocate(alu(jmax,im,im),kp(jmax,im))

       alu=0.0D0
       do i=1,im
          alu(:,i,i) = 1.0D0
       enddo

       do i=1,im
          ag_work = 0.0D0
          ag_work(:,i) = 1.0D0
          ag_work = ag_aq(aq_rDr_aq(aq_ag(ag_work)))/spread(g_R,1,jmax)
          alu(:,im,i) = ag_work(:,im)
       enddo

       call ludecomp(alu,kp)
    endif

    ag_data(:,im) = value0

    ag_data = lusolve(alu,kp,ag_data)

  end subroutine ag_BoundaryGrid_N_2d

  subroutine ag_BoundaryGrid_N_1d(g_data,value)
    !
    ! Dirichlet/Neumann 型境界条件の適用(選点法, 1 次元配列用)
    ! * i=0 で勾配の値を与える.
    !
    real(8), dimension(im),intent(inout)       :: g_data
    !(inout) 境界条件を適用するチェビシェフデータ(0:km)

    real(8), intent(in), optional                :: value
    !(in) 境界値

    real(8), dimension(1,im)                   :: ag_work
    real(8), dimension(1)                        :: vwork           ! 境界値

    if (.not. present(value)) then
       vwork(1)=0.0D0
    else
       vwork(1) = value
    endif

    ag_work(1,:)=g_data
    call ag_BoundaryGrid_N_2d(ag_work,vwork)
    g_data=ag_work(1,:)

  end subroutine ag_BoundaryGrid_N_1d

!---- 下部ルーチン
  subroutine gauss_radau(M,g_R,g_W)
    !
    ! M 次の多項式 π(r) の零点から
    ! ガウス-ラダウ格子点とガウス重みを計算する
    !
    integer, intent(in)  :: M
    ! (in) 多項式打切次数(0 < M, M even)
    real(8), intent(out) :: g_R(M/2)         ! ガウス-ラダウ格子点
    real(8), intent(out) :: g_W(M/2)         ! ガウス重み

    real(8) :: r1,r2,dr
    real(8) :: pi1, pi2
    integer(8) :: i
    integer :: ir
    character(len=4) :: CNUM
    real(8), parameter :: eps = 1.0D-15      ! 打切誤差
    integer, parameter :: nrfact = 100       ! 格子点走査間隔を定める数

    !--- 入力チェック
    if ( mod(M,2) .ne. 0 ) then
       call MessageNotify('E','gauss_radau','M must be even')
    endif

    !--- ガウス--ラダウ格子点
    dr = 1.0D0/(nrfact*M)

    i = 0
    do ir = 1, nrfact*M
       r1  = ir*dr     ;  r2  = (ir+1)*dr
       pi1 = pir(r1,M) ;  pi2 = pir(r2,M)

       if ( pi1*pi2 .le. 0.0D0 ) then
          i = i+1
          g_R(i) = bisec_pir( r1, r2, M , eps )

          if ( i .eq. M/2-1 ) then
             i = i+1
             g_R(M/2) = 1.0D0
             goto 10
          endif
       end if
    enddo

    write( CNUM,'(I4)') i
    call MessageNotify('E','gauss_radau',&
         'Only '//CNUM//' Gauss-Radau points are found')

    !--- ガウス重み
10  continue
    write( CNUM,'(I4)') i
    call MessageNotify('M','gauss_radau',CNUM//' Gauss-Radau points are set up.')

    do i=1,M/2-1
       g_W(i) = 2*(2*M+gamma-5)*g_R(i)**2     &
                /((M+gamma-beta-2)*(M+gamma-3)*Phi(g_R(i),M-2,0)**2)
    enddo

    g_W(M/2) = Inm(0,0) - sum(g_W(1:M/2-1))

  end subroutine gauss_radau

  function bisec_pir( r1, r2, n, eps )
    !
    !  2 分法で π(r,n) の零点を区間 [r1,r2] の間で誤差 eps で探す
    !
    real(8)             :: r1, r2
    integer, intent(in) :: n
    real(8), intent(in) :: eps
    real(8)             :: bisec_pir

    real(8) :: rm, pim
    real(8) :: pi1, pi2

    pi1 = pir(r1,n)
    if ( abs(pi1) .lt. eps ) then
       bisec_pir = r1
       return
    end if

    pi2 = pir(r2,n)
    if ( abs(pi2) .lt. eps ) then
       bisec_pir = r2
       return
    end if

10  rm = (r1+r2)/2
    pim = pir(rm,n)

    if ( abs(pim) .lt. eps ) then
       bisec_pir = rm
       return
    end if

    if ( abs(r1-r2) .lt. eps ) then
       bisec_pir = rm
       return
    end if

    if ( pim * pi1 .gt. 0.0D0 ) then
       r1 = rm
       pi1 = pim
    else
       r2 = rm
       pi2 = pim
    endif
    goto 10

  end function bisec_pir

  function pir(r,n)
    !
    ! ガウス-ラダウ格子点を定義するための関数
    !
    real(8), intent(in)  :: r                ! 動径座標 ( 0 < r < ra )
    integer, intent(in)  :: n                ! 多項式次数(0 <= n)
    real(8)              :: pir              ! 多項式の値

    pir = qnm(r,n,0)-(n+gamma-beta-2)*(n+gamma-3)*qnm(r,n-2,0)/(n*(n+beta-1))

  end function pir

  function Phi(r,n_ind,m_ind)
    !
    ! 正規化された多項式計算
    !
    real(8), intent(in)  :: r                ! 動径座標 ( 0 < r < ra )
    integer, intent(in)  :: n_ind            ! 多項式次数(0 <= |m| <= n)
    integer, intent(in)  :: m_ind            ! 多項式次数(0 <= |m| <= n)
    real(8)              :: Phi              ! 多項式の値

    Phi = qnm(r,n_ind,m_ind)/sqrt(Inm(n_ind,m_ind))
  end function Phi

  function qnm(r,n_ind,m_ind)
    !
    ! 漸化式による多項式計算
    !
    real(8), intent(in)  :: r                ! 動径座標 ( 0 < r < ra )
    integer, intent(in)  :: n_ind            ! 多項式次数(0 <= |m| <= n)
    integer, intent(in)  :: m_ind            ! 多項式次数(0 <= |m| <= n)
    real(8)              :: qnm              ! 多項式の値

    real(8) :: qnp2m, qnm2m
    integer :: m, n

    m = abs(m_ind)

    !--- 入力チェック
    if ( n_ind .lt. 0 ) then
       call MessageNotify('E','qnm','n must be larger equal to 0')
    endif

    if ( m .gt. n_ind ) then
       call MessageNotify('E','qnm','abs(m) must be smaller equal to n')
    endif

    if ( mod(n_ind+m_ind,2) .ne. 0 ) then
       call MessageNotify('E','qnm','n+m must be even')
    endif

    !--- 初期値
    if ( m .eq. n_ind ) then
       qnm = r**m
       return
    endif

    if ( n_ind .eq. m+2 ) then
       qnm = (2*m+gamma-1)/2*((2*m+gamma+1)/(2*m+beta+1)*r**2 -1) * r**m
       return
    endif

    qnm2m = r**m
    qnm = (2*m+gamma-1)/2*((2*m+gamma+1)/(2*m+beta+1)*r**2 -1)* r**m

    !--- 漸化式

    do n = m+2, n_ind-2, 2
       qnp2m =( ( 2*n + gamma - 1 ) &
                * (  ( 2*n + gamma - 3 )*( 2*n + gamma + 1) * r**2       &
                    - 2*n * (n + gamma -1) - 2*m * (m + beta -1)       &
                    -(gamma-3) * (beta+1)                            )*qnm &
               -(n-m+gamma-beta-2)*(n+m+gamma-3)*(2*n+gamma+1)*qnm2m )&
              /(( n - m + 2 )*( n + m + beta + 1 )*( 2*n + gamma - 3 ))
       qnm2m = qnm
       qnm = qnp2m
    enddo

    return
  end function qnm

  function Inm(n_ind,m_ind)
    !
    ! 多項式の規格化定数を漸化式により求める
    !
    integer, intent(in)  :: n_ind            ! 多項式次数(0 <= |m| <= n)
    integer, intent(in)  :: m_ind            ! 多項式次数(0 <= |m| <= n)
    real(8)              :: Inm              ! 規格化定数

    integer :: n,m
    real(8) :: Inm2m
    real(8) :: Inmln                         ! 規格化定数のlog
    real(8) :: gammaln
    external gammaln

    m=abs(m_ind)

    !--- 入力チェック
    if ( n_ind .lt. 0 ) then
       call MessageNotify('E','Inm','n must be larger equal to 0')
    endif

    if ( m .gt. n_ind ) then
       call MessageNotify('E','Inm','abs(m) must be smaller equal to n')
    endif

    if ( mod(n_ind+m_ind,2) .ne. 0 ) then
       call MessageNotify('E','Inm','n+m must be even')
    endif

    !--- 初期値
    Inmln = gammaln((gamma-beta)/2) + gammaln(m+(beta+1)/2) &
          - gammaln(m+(gamma+1)/2) - log(2.0D0)
    Inm = exp(Inmln)

    !--- 漸化式
    do n=m+2, n_ind, 2
       Inm2m = Inm
       Inm = (2*n+gamma-5)*(n-m+gamma-beta-2)*(n+m+gamma-3) &
               /((n-m)*(n+m+beta-1)*(2*n+gamma-1)) * Inm2m
    enddo

    return
  end function Inm

end module aq_module
