fftについて

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