!-- !---------------------------------------------------------------------- ! Copyright(c) 2008-2025 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_Initialize で ! 設定する重みの指数の配列 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(DP), dimension(im). ! * im は R 座標の格子点数であり, サブルーチン aq_Initial にて ! あらかじめ設定しておく. ! ! * q : スペクトルデータ. ! * 変数の種類と次元は real(DP), dimension(0:km). ! * km は R 方向の最大波数であり, サブルーチン aq_Initial にて ! あらかじめ設定しておく. スペクトルデータのq(l) には ! nd+2l 次の多項式成分が格納されている. ! ! * ag : 1 次元(R)格子点データの並んだ 2 次元データ. ! * 変数の種類と次元は real(DP), dimension(size(nd),im). ! 第 2 次元が R 方向を表す. ! ! * aq : 1 次元スペクトルデータの並んだ 2 次元データ. ! * 変数の種類と次元は real(DP), 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_types, only: DP use dc_message use lumatrix implicit none private public aq_Initial public aq_Initialize ! 初期化 public aq_Finalize ! 終了処理 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_Initialize module procedure aq_Initialize_I4 module procedure aq_Initialize_I8 end interface aq_Initialize interface aq_Boundary_D ! ! 外側境界ディリクレ型境界条件の適用(タウ法). ! ! * 境界条件を適用する配列の次元によって内部でサブルーチンを ! 使い分けている. ユーザーインターフェースは共通であるので ! 下部ルーチンを呼ぶ必要はない. ! ! 引数と結果の型 ! ! * 1 次元スペクトルデータの並んだ 2 次元配列の場合 ! ! real(DP), dimension(size(md),0:km),intent(inout) :: aq_data ! !(inout) 境界条件を適用するスペクトルデータ ! ! real(DP), dimension(size(md)), intent(in), optional :: value ! !(in) 適用する境界値 ! ! * 1 次元チェビシェフデータの場合 ! ! real(DP), dimension(0:km),intent(inout) :: q_data ! !(inout) 境界条件を適用するスペクトルデータ ! ! real(DP), intent(in), optional :: value ! !(in) 適用する境界値 ! module procedure aq_BoundaryTau_D_1d, aq_BoundaryTau_D_2d end interface aq_Boundary_D interface aq_Boundary_N ! ! 外側境界ノイマン型境界条件の適用(タウ法). ! ! * 境界条件を適用する配列の次元によって内部でサブルーチンを ! 使い分けている. ユーザーインターフェースは共通であるので ! 下部ルーチンを呼ぶ必要はない. ! ! 引数と結果の型 ! ! * 1 次元スペクトルデータの並んだ 2 次元配列の場合 ! ! real(DP), dimension(size(md),0:km),intent(inout) :: aq_data ! !(inout) 境界条件を適用するスペクトルデータ ! ! real(DP), dimension(size(md)), intent(in), optional :: value ! !(in) 適用する境界値 ! ! * 1 次元チェビシェフデータの場合 ! ! real(DP), dimension(0:km),intent(inout) :: q_data ! !(inout) 境界条件を適用するスペクトルデータ ! ! real(DP), intent(in), optional :: value ! !(in) 適用する境界値 ! module procedure aq_BoundaryTau_N_1d, aq_BoundaryTau_N_2d end interface aq_Boundary_N interface aq_BoundaryTau_D ! ! 外側境界ディリクレ型境界条件の適用(タウ法). ! ! * 境界条件を適用する配列の次元によって内部でサブルーチンを ! 使い分けている. ユーザーインターフェースは共通であるので ! 下部ルーチンを呼ぶ必要はない. ! ! 引数と結果の型 ! ! * 1 次元スペクトルデータの並んだ 2 次元配列の場合 ! ! real(DP), dimension(size(md),0:km),intent(inout) :: aq_data ! !(inout) 境界条件を適用するチェビシェフデータ ! ! real(DP), dimension(size(md)), intent(in), optional :: value ! !(in) 適用する境界値 ! ! * 1 次元スペクトルデータの場合 ! ! real(DP), dimension(0:km),intent(inout) :: q_data ! !(inout) 境界条件を適用するチェビシェフデータ ! ! real(DP), intent(in), optional :: value ! !(in) 適用する境界値 ! module procedure aq_BoundaryTau_D_1d, aq_BoundaryTau_D_2d end interface aq_BoundaryTau_D interface aq_BoundaryTau_N ! ! 外側境界ノイマン型境界条件の適用(タウ法). ! ! * 境界条件を適用する配列の次元によって内部でサブルーチンを ! 使い分けている. ユーザーインターフェースは共通であるので ! 下部ルーチンを呼ぶ必要はない. ! ! 引数と結果の型 ! ! * 1 次元スペクトルデータの並んだ 2 次元配列の場合 ! ! real(DP), dimension(size(md),0:km),intent(inout) :: aq_data ! !(inout) 境界条件を適用するチェビシェフデータ ! ! real(DP), dimension(size(md)), intent(in), optional :: value ! !(in) 適用する境界値 ! ! * 1 次元スペクトルデータの場合 ! ! real(DP), dimension(0:km),intent(inout) :: q_data ! !(inout) 境界条件を適用するチェビシェフデータ ! ! real(DP), intent(in), optional :: value ! !(in) 適用する境界値 ! module procedure aq_BoundaryTau_N_1d, aq_BoundaryTau_N_2d end interface aq_BoundaryTau_N interface ag_BoundaryGrid_D ! ! 外側境界ディリクレ型境界条件の適用(選点法). ! ! * 境界条件を適用する配列の次元によって内部でサブルーチンを ! 使い分けている. ユーザーインターフェースは共通であるので ! 下部ルーチンを呼ぶ必要はない. ! ! 引数と結果の型 ! ! * 1 次元スペクトルデータの並んだ 2 次元配列の場合 ! ! real(DP), dimension(size(md),0:km),intent(inout) :: aq_data ! !(inout) 境界条件を適用するチェビシェフデータ ! ! real(DP), dimension(size(md)), intent(in), optional :: value ! !(in) 適用する境界値 ! ! * 1 次元スペクトルデータの場合 ! ! real(DP), dimension(0:km),intent(inout) :: q_data ! !(inout) 境界条件を適用するチェビシェフデータ ! ! real(DP), intent(in), optional :: value ! !(in) 適用する境界値 ! module procedure ag_BoundaryGrid_D_1d, ag_BoundaryGrid_D_2d end interface ag_BoundaryGrid_D interface ag_BoundaryGrid_N ! ! 外側境界ノイマン型境界条件の適用(選点法). ! ! * 境界条件を適用する配列の次元によって内部でサブルーチンを ! 使い分けている. ユーザーインターフェースは共通であるので ! 下部ルーチンを呼ぶ必要はない. ! ! 引数と結果の型 ! ! * 1 次元スペクトルデータの並んだ 2 次元配列の場合 ! ! real(DP), dimension(size(md),0:km),intent(inout) :: aq_data ! !(inout) 境界条件を適用するチェビシェフデータ ! ! real(DP), dimension(size(md)), intent(in), optional :: value ! !(in) 適用する境界値 ! ! * 1 次元スペクトルデータの場合 ! ! real(DP), dimension(0:km),intent(inout) :: q_data ! !(inout) 境界条件を適用するチェビシェフデータ ! ! real(DP), intent(in), optional :: value ! !(in) 適用する境界値 ! module procedure ag_BoundaryGrid_N_1d, ag_BoundaryGrid_N_2d end interface ag_BoundaryGrid_N integer(8) :: im, km ! 格子点数, 切断波数 real(DP) :: ra ! 領域の大きさ real(DP) :: alpha ! 展開多項式パラメター 0 < α <= 1 real(DP) :: beta ! 展開多項式パラメター 0 < β real(DP) :: gamma ! 展開多項式パラメター γ=2α+β integer(8), allocatable :: md(:) ! 展開多項式上付次数のならび integer(8) :: jmax ! a 座標(データ第 1 次元)の大きさ real(DP), allocatable :: g_R(:) ! ガウス--ラダウ格子点 real(DP), allocatable :: g_R_Weight(:) ! ガウス重み real(DP), allocatable :: CF(:, :, :) ! 正変換用行列 real(DP), 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, k, r, a, b, m) integer, intent(in) :: i, k integer, dimension(:), intent(in) :: m real(DP), intent(in) :: a, b, r call MessageNotify('W', 'aq_Initial', & & 'aq_Initial is deprecated; will be removed; use aq_Initialize.') call aq_Initialize_I8(int(i, kind(im)), int(k, kind(km)), r, a, b, int(m, kind(md))) end subroutine aq_Initial ! -------------------------------------------------------------------------- subroutine aq_Initialize_I4(i, k, r, a, b, m) integer, intent(in) :: i, k integer, dimension(:), intent(in) :: m real(DP), intent(in) :: a, b, r call MessageNotify('W', 'aq_Initialize', & & 'INTEGER(4) im/km/md deprecated; will be removed; use INTEGER(8).') call aq_Initialize_I8(int(i, kind(im)), int(k, kind(km)), r, a, b, int(m, kind(md))) end subroutine aq_Initialize_I4 !--- 初期化 ---------------------------------------------------------------- subroutine aq_Initialize_I8(i_in, k_in, r_in, alpha_in, beta_in, md_in) ! ! スペクトル変換の格子点数, 波数, 領域の大きさ, 重みを設定する. ! ! 他の関数や変数を呼ぶ前に, 最初にこのサブルーチンを呼んで ! 初期設定をしなければならない. ! integer(8), intent(in) :: i_in !(in) 格子点数 integer(8), intent(in) :: k_in !(in) 切断波数 real(DP), intent(in) :: r_in !(in) 外側境界の座標(半径) real(DP), intent(in) :: alpha_in !(in) 展開多項式パラメター real(DP), intent(in) :: beta_in !(in) 展開多項式パラメター integer(8), intent(in) :: md_in(:) !(in) 展開多項式上付次数のならび integer(8) i, j, n, nn !--- パラメターのチェック・記憶 im = i_in km = k_in ra = r_in alpha = alpha_in beta = beta_in if (km .ge. 2.0_DP * im) then call MessageNotify('E', 'aq_initial', 'KM shoud be less than 2*IM') end if if (alpha .le. 0.0_DP) then call MessageNotify('E', 'aq_initial', 'alpha must be larger than 0') end if if (alpha .gt. 1.0_DP) then call MessageNotify('E', 'aq_initial', 'alpha must be smaller equal to 1') end if if (beta .le. 0.0D0) then call MessageNotify('E', 'aq_initial', 'beta must be larger than 0') end if gamma = 2.0_DP * 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 = 0, km nn = n * 2 + md(j) do i = 1, im CF(j, n, i) = Phi(g_R(i) / ra, nn, md(j)) * g_R_Weight(i) / ra**(gamma - 1) end do end do end do CB = 0.0D0 do j = 1, jmax do i = 1, im do n = 0, km nn = n * 2 + md(j) CB(j, i, n) = Phi(g_R(i) / ra, nn, md(j)) end do end do end do !--- 各ルーチン変換行列の初期化スイッチ first_r2inv = .true. first_Tau_D = .true. first_Tau_N = .true. first_Grid_N = .true. call MessageNotify( & & 'M', 'aq_Initialize', 'aq_module (2026/02/12) is initialized') end subroutine aq_Initialize_I8 subroutine aq_Finalize first_r2inv = .false. first_Tau_D = .false. first_Tau_N = .false. first_Grid_N = .false. end subroutine aq_Finalize !---- 逆変換 ---- function ag_aq(aq_data) ! ! スペクトルデータから格子データへ変換する(2 次元配列用). ! real(DP), dimension(:, 0:), intent(in) :: aq_data !(in) スペクトルデータ real(DP), dimension(size(aq_data, 1), im) :: ag_aq !(out) 格子点データ integer(8) :: 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.') end if ag_aq = 0.0D0 do i = 1, im !$omp parallel do do j = 1, jmax ag_aq(j, i) = sum(CB(j, i, :) * aq_data(j, :)) end do !$omp end parallel do end do end function ag_aq function g_q(q_data) ! ! スペクトルデータから格子データへ変換する(1 次元配列用). ! real(DP), dimension(:), intent(in) :: q_data !(in) スペクトルデータ real(DP), dimension(im) :: g_q !(out) 格子点データ real(DP), dimension(1, size(q_data)) :: q_work ! 作業用配列 real(DP), 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(DP), dimension(:, :), intent(in) :: ag_data !(in) 格子点データ real(DP), dimension(size(ag_data, 1), 0:km) :: aq_ag !(out) スペクトルデータ integer(8) :: 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.') end if aq_ag = 0.0D0 !$omp parallel do private(n) do j = 1, jmax do n = 0, km aq_ag(j, n) = sum(CF(j, n, :) * ag_data(j, :)) end do end do !$omp end parallel do end function aq_ag function q_g(g_data) ! ! 格子データからスペクトルデータへ変換する(1 次元配列用). ! real(DP), dimension(:), intent(in) :: g_data !(in) 格子点データ real(DP), dimension(0:km) :: q_g !(out) スペクトルデータ real(DP), dimension(1, size(g_data)) :: ag_work real(DP), 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(DP), dimension(:, 0:), intent(in) :: aq_data !(in) 入力スペクトルデータ real(DP), dimension(size(aq_data, 1), 0:size(aq_data, 2) - 1) :: aq_rDr_aq !(out) 出力スペクトルデータ integer(8) :: j, n, nn real(DP) :: 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.') end if aq_rDr_aq = 0.0D0 !$omp parallel do private(sqInm,sqInm2,n,nn) do j = 1, jmax n = md(j) + km * 2 sqInm = sqrt(Inm(n, md(j))) aq_rDr_aq(j, km) = n * aq_data(j, km) do nn = km - 1, 0, -1 n = md(j) + nn * 2 sqInm2 = sqInm sqInm = sqrt(Inm(n, md(j))) aq_rDr_aq(j, nn) = sqInm / sqInm2 * ( & (2 * n + gamma - 1) / (2 * n + gamma + 3) * aq_rDr_aq(j, nn + 1) & + (2 * n + gamma - 1) * (n + gamma + 1) / (2 * n + gamma + 3) * aq_data(j, nn + 1)) & + n * aq_data(j, nn) end do end do !$omp end parallel do end function aq_rDr_aq function q_rDr_q(q_data) ! ! 入力スペクトルデータに r(d/dR) 微分を作用する(1 次元配列用). ! ! スペクトルデータの r(d/dR) 微分とは, 対応する格子点データに R 微分を ! 作用させたデータのスペクトル変換のことである. ! ! real(DP), dimension(:), intent(in) :: q_data !(in) 入力チェビシェフデータ real(DP), dimension(0:km) :: q_rDr_q !(out) チェビシェフデータの R 微分 real(DP), 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(DP), dimension(:, 0:), intent(in) :: aq_data !(in) 入力スペクトルデータ real(DP), dimension(size(aq_data, 1), 0:size(aq_data, 2) - 1) :: aq_r2_aq !(out) 出力スペクトルデータ integer(8) :: j, n, m, nn real(DP) :: 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.') end if aq_r2_aq = 0.0D0 !$omp parallel do private(sqrInm,sqrInm2m,sqrInp2m,nn,n,m) do j = 1, jmax m = abs(md(j)) n = md(j) sqrInm = sqrt(Inm(n, m)) sqrInp2m = sqrt(Inm(n + 2, m)) aq_r2_aq(j, 0) = & (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, 0) & + (n - m + gamma - beta) * (n + m + gamma - 1) & / ((2 * n + gamma + 3) * (2 * n + gamma + 1)) & * aq_data(j, 1) * (sqrInm / sqrInp2m) do nn = 1, km - 1 n = md(j) + nn * 2 sqrInm2m = sqrInm sqrInm = sqrInp2m sqrInp2m = sqrt(Inm(n + 2, m)) aq_r2_aq(j, nn) = & (n - m) * (n + m + beta - 1) / ((2 * n + gamma - 5) * (2 * n + gamma - 3)) & * aq_data(j, nn - 1) * (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, nn) & + (n - m + gamma - beta) * (n + m + gamma - 1) & / ((2 * n + gamma + 3) * (2 * n + gamma + 1)) & * aq_data(j, nn + 1) * (sqrInm / sqrInp2m) end do n = md(j) + km * 2 sqrInm2m = sqrInm sqrInm = sqrInp2m aq_r2_aq(j, km) = & (n - m) * (n + m + beta - 1) / ((2 * n + gamma - 5) * (2 * n + gamma - 3)) & * aq_data(j, km - 1) * (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, km) end do !$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(DP), dimension(:), intent(in) :: q_data !(in) 入力チェビシェフデータ real(DP), dimension(0:km) :: q_r2_q !(out) チェビシェフデータの R 微分 real(DP), 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(DP), dimension(:, 0:), intent(in) :: aq_data !(in) 入力スペクトルデータ real(DP), dimension(size(aq_data, 1), 0:size(aq_data, 2) - 1) :: aq_r2Inv_aq !(out) 出力スペクトルデータ real(DP), dimension(:, :, :), allocatable :: R2MTX integer, dimension(:, :), allocatable :: kp integer(8) :: j, n, m, nn real(DP) :: 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.') end if 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 end do !$omp parallel do private(sqrInm,sqrInm2m,sqrInp2m,n,m,nn) do j = 1, jmax m = abs(md(j)) n = md(j) sqrInm = sqrt(Inm(n, m)) sqrInp2m = sqrt(Inm(n + 2, m)) R2MTX(j, 0, 0) = & (2 * n * (n + gamma - 1) + 2 * m * (m + beta - 1) + (gamma - 3) * (beta + 1)) & / ((2 * n + gamma + 1) * (2 * n + gamma - 3)) R2MTX(j, 0, 1) = & (n - m + gamma - beta) * (n + m + gamma - 1) & / ((2 * n + gamma + 3) * (2 * n + gamma + 1)) & * (sqrInm / sqrInp2m) do nn = 1, km - 1 n = md(j) + nn * 2 sqrInm2m = sqrInm sqrInm = sqrInp2m sqrInp2m = sqrt(Inm(n + 2, m)) R2MTX(j, nn, nn - 1) = & (n - m) * (n + m + beta - 1) / ((2 * n + gamma - 5) * (2 * n + gamma - 3)) & * (sqrInm / sqrInm2m) R2MTX(j, nn, nn) = & (2 * n * (n + gamma - 1) + 2 * m * (m + beta - 1) + (gamma - 3) * (beta + 1)) & / ((2 * n + gamma + 1) * (2 * n + gamma - 3)) R2MTX(j, nn, nn + 1) = & (n - m + gamma - beta) * (n + m + gamma - 1) & / ((2 * n + gamma + 3) * (2 * n + gamma + 1)) & * (sqrInm / sqrInp2m) end do n = md(j) + km * 2 sqrInm2m = sqrInm sqrInm = sqrInp2m R2MTX(j, km, km - 1) = & (n - m) * (n + m + beta - 1) / ((2 * n + gamma - 5) * (2 * n + gamma - 3)) & * (sqrInm / sqrInm2m) R2MTX(j, km, km) = & (2 * n * (n + gamma - 1) + 2 * m * (m + beta - 1) + (gamma - 3) * (beta + 1)) & / ((2 * n + gamma + 1) * (2 * n + gamma - 3)) end do !$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(DP), dimension(:), intent(in) :: q_data !(in) 入力チェビシェフデータ real(DP), dimension(0:km) :: q_r2Inv_q !(out) チェビシェフデータの R 微分 real(DP), 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(DP), dimension(:, :), intent(in) :: ag !(in)入力格子点データ real(DP), dimension(size(ag, 1)) :: a_Int_ag !(out) 積分したデータ integer(8) :: 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.') end if a_Int_ag = 0.0d0 do i = 1, im a_Int_ag(:) = a_Int_ag(:) + ag(:, i) * g_R_Weight(i) end do 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(DP), dimension(im), intent(in) :: g !(in) 格子点データ real(DP) :: 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(DP), dimension(:, 0:), intent(in) :: ag !(in)入力格子点データ real(DP), 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(DP), dimension(im), intent(in) :: g !(in) 格子点データ real(DP) :: 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(DP), dimension(:, 0:), intent(inout) :: aq_data !(inout) 境界条件を適用するチェビシェフデータ(jmax,0:km) real(DP), dimension(:), intent(in), optional :: value !(in) 境界値(jmax) real(DP), dimension(:, :, :), allocatable :: alu integer, dimension(:, :), allocatable :: kp real(DP), dimension(size(aq_data, 1), 0:km) :: aq_work real(DP), dimension(size(aq_data, 1), im) :: ag_work real(DP), dimension(size(aq_data, 1)) :: value0 ! 境界値 integer(8) :: k, j 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.') end if if (.not. present(value)) then value0 = 0.0D0 else value0 = value end if 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 end do do k = 0, km aq_work = 0.0D0 aq_work(:, k) = 1.0D0 ag_work = ag_aq(aq_work) alu(:, km, k) = ag_work(:, im) end do call LUDECOMP(alu, kp) end if do j = 1, jmax aq_data(j, km) = value0(j) end do aq_data = LUSOLVE(alu, kp, aq_data) end subroutine aq_BoundaryTau_D_2d subroutine aq_BoundaryTau_D_1d(q_data, value) ! ! Dirichlet 型境界条件の適用(タウ法, 1 次元配列用) ! * 両境界での値を与える. ! real(DP), dimension(0:km), intent(inout) :: q_data !(inout) 境界条件を適用するチェビシェフデータ(0:km) real(DP), intent(in), optional :: value !(in) 境界値 real(DP), dimension(1, 0:km) :: aq_work real(DP), dimension(1) :: vwork ! 境界値 if (.not. present(value)) then vwork(1) = 0.0D0 else vwork(1) = value end if 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(DP), dimension(:, 0:), intent(inout) :: aq_data !(inout) 境界条件を適用するチェビシェフデータ(m,0:km) real(DP), dimension(:), intent(in), optional :: value !(in) 境界値(m) real(DP), dimension(:, :, :), allocatable :: alu integer, dimension(:, :), allocatable :: kp real(DP), dimension(size(aq_data, 1), 0:km) :: aq_work real(DP), dimension(size(aq_data, 1), im) :: ag_work real(DP), dimension(size(aq_data, 1)) :: value0 ! 境界値 integer(8) :: j, k 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.') end if if (.not. present(value)) then value0 = 0.0D0 else value0 = value end if 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 end do do k = 0, km aq_work = 0.0D0 aq_work(:, k) = 1.0D0 ag_work = ag_aq(aq_rDr_aq(aq_work)) / spread(g_R, 1, jmax) alu(:, km, k) = ag_work(:, im) end do call LUDECOMP(alu, kp) end if do j = 1, jmax aq_data(j, km) = value0(j) end do 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(DP), dimension(0:km), intent(inout) :: q_data !(inout) 境界条件を適用するチェビシェフデータ(0:km) real(DP), intent(in), optional :: value !(in) 境界値 real(DP), dimension(1, 0:km) :: aq_work real(DP), dimension(1) :: vwork ! 境界値 if (.not. present(value)) then vwork(1) = 0.0D0 else vwork(1) = value end if 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(DP), dimension(:, :), intent(inout) :: ag_data !(inout) 境界条件を適用するチェビシェフデータ(jmax,im) real(DP), dimension(:), intent(in), optional :: value !(in) 境界値(jmax) real(DP), dimension(size(ag_data, 1)) :: value0 if (.not. present(value)) then value0 = 0.0d0 else value0 = value end if ag_data(:, im) = value0 end subroutine ag_BoundaryGrid_D_2d subroutine ag_BoundaryGrid_D_1d(g_data, value) ! ! Dirichlet 型境界条件の適用(タウ法, 1 次元配列用) ! * 両境界での値を与える. ! real(DP), dimension(im), intent(inout) :: g_data !(inout) 境界条件を適用するチェビシェフデータ(im) real(DP), intent(in), optional :: value !(in) 境界値 real(DP), dimension(1, im) :: ag_work real(DP), dimension(1) :: vwork if (.not. present(value)) then vwork(1) = 0.0d0 else vwork(1) = value end if 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(DP), dimension(:, :), intent(inout) :: ag_data !(inout) 境界条件を適用するチェビシェフデータ(m,im) real(DP), dimension(:), intent(in), optional :: value !(in) 境界値(m) real(DP), dimension(:, :, :), allocatable :: alu integer, dimension(:, :), allocatable :: kp real(DP), dimension(size(ag_data, 1), im) :: ag_work real(DP), dimension(size(ag_data, 1)) :: value0 ! 境界値 integer(8) :: 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.') end if if (.not. present(value)) then value0 = 0.0D0 else value0 = value end if 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 end do 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) end do call LUDECOMP(alu, kp) end if 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(DP), dimension(im), intent(inout) :: g_data !(inout) 境界条件を適用するチェビシェフデータ(0:km) real(DP), intent(in), optional :: value !(in) 境界値 real(DP), dimension(1, im) :: ag_work real(DP), dimension(1) :: vwork ! 境界値 if (.not. present(value)) then vwork(1) = 0.0D0 else vwork(1) = value end if 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(8), intent(in) :: M ! (in) 多項式打切次数(0 < M, M even) real(DP), intent(out) :: g_R(M / 2) ! ガウス-ラダウ格子点 real(DP), intent(out) :: g_W(M / 2) ! ガウス重み real(DP) :: r1, r2, dr real(DP) :: pi1, pi2 integer(8) :: i integer(8) :: ir character(len=4) :: CNUM real(DP), parameter :: eps = 1.0D-15 ! 打切誤差 integer, parameter :: nrfact = 100 ! 格子点走査間隔を定める数 !--- 入力チェック if (mod(M, 2_8) .ne. 0) then call MessageNotify('E', 'gauss_radau', 'M must be even') end if !--- ガウス--ラダウ格子点 dr = 1.0_DP / (real(nrfact, DP) * real(M, DP)) 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 end if end if end do 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_8)**2) end do g_W(M / 2) = Inm(0_8, 0_8) - 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(DP) :: r1, r2 integer(8), intent(in) :: n real(DP), intent(in) :: eps real(DP) :: bisec_pir real(DP) :: rm, pim real(DP) :: 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 end if goto 10 end function bisec_pir function pir(r, n) ! ! ガウス-ラダウ格子点を定義するための関数 ! real(DP), intent(in) :: r ! 動径座標 ( 0 < r < ra ) integer(8), intent(in) :: n ! 多項式次数(0 <= n) real(DP) :: pir ! 多項式の値 pir = qnm(r, n, 0_8) - (n + gamma - beta - 2) * (n + gamma - 3) * qnm(r, n - 2, 0_8) / (n * (n + beta - 1)) end function pir function Phi(r, n_ind, m_ind) ! ! 正規化された多項式計算 ! real(DP), intent(in) :: r ! 動径座標 ( 0 < r < ra ) integer(8), intent(in) :: n_ind ! 多項式次数(0 <= |m| <= n) integer(8), intent(in) :: m_ind ! 多項式次数(0 <= |m| <= n) real(DP) :: 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(DP), intent(in) :: r ! 動径座標 ( 0 < r < ra ) integer(8), intent(in) :: n_ind ! 多項式次数(0 <= |m| <= n) integer(8), intent(in) :: m_ind ! 多項式次数(0 <= |m| <= n) real(DP) :: qnm ! 多項式の値 real(DP) :: qnp2m, qnm2m integer(8) :: m, n m = abs(m_ind) !--- 入力チェック if (n_ind .lt. 0) then call MessageNotify('E', 'qnm', 'n must be larger equal to 0') end if if (m .gt. n_ind) then call MessageNotify('E', 'qnm', 'abs(m) must be smaller equal to n') end if if (mod(n_ind + m_ind, 2_8) .ne. 0) then call MessageNotify('E', 'qnm', 'n+m must be even') end if !--- 初期値 if (m .eq. n_ind) then qnm = r**m return end if 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 end if 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 end do return end function qnm function Inm(n_ind, m_ind) ! ! 多項式の規格化定数を漸化式により求める ! integer(8), intent(in) :: n_ind ! 多項式次数(0 <= |m| <= n) integer(8), intent(in) :: m_ind ! 多項式次数(0 <= |m| <= n) real(DP) :: Inm ! 規格化定数 integer(8) :: n, m real(DP) :: Inm2m real(DP) :: Inmln ! 規格化定数のlog real(DP) :: dlgamma external dlgamma m = abs(m_ind) !--- 入力チェック if (n_ind .lt. 0) then call MessageNotify('E', 'Inm', 'n must be larger equal to 0') end if if (m .gt. n_ind) then call MessageNotify('E', 'Inm', 'abs(m) must be smaller equal to n') end if if (mod(n_ind + m_ind, 2_8) .ne. 0) then call MessageNotify('E', 'Inm', 'n+m must be even') end if !--- 初期値 Inmln = dlgamma((gamma - beta) / 2) + dlgamma(m + (beta + 1) / 2) & & - dlgamma(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 end do return end function Inm end module aq_module