乱流のエネルギースペクトル関数を求めるためにfftを勉強しているのですが，fortranにおけるfftライブラリは区間0~1で積分したものが返ってくるのでしょうか？また一般の系の長さをLとしたときは区間0~Lで積分したものが欲しいので2Pi/Lを割ったりかけたりする必要があるのでしょうか？
ざっくりとした質問で申し訳ありませんがわかる範囲で何かアドバイスあればよろしくお願いします．
例）
Ruby on RailsでQiitaのようなWebアプリをつくっています。
記事を投稿する機能の実装中にエラーが発生しました。
解決方法を教えて下さい。
例）
Ruby on RailsでQiitaのようなWebアプリをつくっています。
記事を投稿する機能の実装中にエラーが発生しました。
解決方法を教えて下さい。
出ているエラーメッセージを入力
例）
NameError (uninitialized constant World)
または、問題・エラーが起きている画像をここにドラッグアンドドロップ
ソースコードを入力
例）
def greet
puts Hello World
end
ここに問題・エラーに対して試したことを記載してください。
Fortranに組み込みのFFTライブラリは無いと認識しているので、一番よく使われるFFTWと仮定してお答えしますが、正規化はされないので信号の長さで除算する必要があります。
https://www.fftw.org/fftw3_doc/The-1d-Real_002ddata-DFT.html
離散フーリエ変換で調べてみるのをおすすめします、Qiitaにもわかりやすい記事がたくさんあります。
ただし-5/3を見たいだけなら正規化は定数にしか効いてこないのでそこまで意識する必要はないかもしれません。
@hattorisena
Questioner@hattorisena
Questioner@hattorisena
Questioner
subroutine energy_sample(ux,uy,istep,sum_Ek)!uxはx方向速度,uyはy方向速度
real(8), intent(in) :: ux(0:NX+1, 0:NY+1), uy(0:NX+1, 0:NY+1)
integer, intent(in) :: istep
complex(8) ux_hat(1:NX/2+1, 1:NY), uy_hat(1:NX/2+1, 1:NY)
real(8) E_tmp(1:NX/2+1, 1:NY)
real(8) K_abs(1:NX/2+1, 1:NY)
real(8) Energy(0:NX)
real(8) sum_Ek,im, Ek
integer(8) i, j, index
integer(8) plan
character(*), parameter :: str1='./128_0.85v4_0.60/enespe/'
character(*), parameter :: str2='.d'
character(8) str
character(64) filename
write(str, '(I8.8)') istep ! 数値を文字列に変換
filename = str1//str//str2 ! 文字列の結合
! 速度場をフーリエ変換
call dfftw_plan_dft_r2c_2d(plan, NX, NY, ux(1:NX,1:NY), ux_hat, FFTW_ESTIMATE) !uxのフーリエ変換を返す
call dfftw_execute(plan, ux(1:NX,1:NY), ux_hat)
call dfftw_destroy_plan(plan)
call dfftw_plan_dft_r2c_2d(plan, NX, NY, uy(1:NX,1:NY), uy_hat, FFTW_ESTIMATE)!uyのフーリエ変換を返す
call dfftw_execute(plan, uy(1:NX,1:NY), uy_hat)
call dfftw_destroy_plan(plan)
ux_hat(:,:) = ux_hat(:,:)/(NX*NY/2)
uy_hat(:,:) = uy_hat(:,:)/(NX*NY/2)
! 二乗を足す
E_tmp(:,:) = 0.0d0
do j = 1, NY
do i = 1, NX/2+1
E_tmp(i, j) = E_tmp(i, j) + real(ux_hat(i, j))**2 + imag(ux_hat(i, j))**2
E_tmp(i, j) = E_tmp(i, j) + real(uy_hat(i, j))**2 + imag(uy_hat(i, j))**2
enddo
enddo
! 配列要素に対する波数 そのまま代入するので0クリアいらない
do j = 1, NY
do i = 1, NX/2+1
if (i==1 .and. j==1) then ! sqrt(0)は計算不可能
K_abs(i, j) = 0.0d0
else
K_abs(i, j) = sqrt((i-1)**2.0d0 + (j-1)**2.0d0) !K_abs 波数の絶対値 2.0d0としないと実数にならずsqrtが処理してくれない
endif
enddo
enddo
! 波数を四捨五入し、対応する整数の波数にエネルギーを足し合わせる
Energy(:) = 0.0d0
do j = 1, NY/2+1 ! 波数は半分しか計算できない
do i = 1, NX/2+1
index = nint(K_abs(i, j)) !nint 実数から小数へ四捨五入 !小数で求まる
Energy(index) = Energy(index) + E_tmp(i, j)/2.0
enddo
enddo
sum_Ek = 0.0d0
do i = 0,91
sum_Ek = sum_Ek + Energy(i)
enddo
!Energy(:) = (Energy(:)*Xmax)/(2.0*2.0*Pi)
!delta_k = 2Pi/Lでエネルギースペクトルを割るとほしいものが求まりそう
open(10, file=filename)
do i = 0, NX
im = i*2*Pi/Xmax
Ek =Energy(i)*Xmax/(2.0d0*Pi)
write(10, '(I4, e12.4)') im, Ek
enddo
@hattorisena
Questioner