# -*- coding: utf-8 -*-

#
# = エネルギー方程式の運動エネルギーの項の領域積分の時間微分を計算するスクリプト
#
# * 履歴 (新しい日付を上に書く)
#   * 2011/10/18  結果を dat ファイルに出力するように編集
#   * 2011/10/13  新規作成
#
# * めも
#   * 手順
#     * 各時刻での Term の領域積分の値を求める
#       * とある時刻での Term を計算
#       * GPhys のコマンド (.integrate) を使って領域積分
#     * 各時刻での Term の領域積分の時間微分を求める 
#     * 結果を result.dat に出力
#   # これで正しいのかはわからない
#
##########################################################################

require "numru/ggraph"
include NumRu

dt = 10.0
tt = 20000
#tt = 200
tnmin = 0
tnmax = tt/dt
tn = 0

x1 = 0
x2 = 10000
z1 = 0
z2 = 10000

###--- 必要な netCDF ファイルを open する ---###

gphys10 = GPhys::IO.open('../../netCDF/thermal-moist_restart.nc','DensBZ')
gphys1  = GPhys::IO.open('../../netCDF/thermal-moist_VelX.nc','VelX')
gphys2  = GPhys::IO.open('../../netCDF/thermal-moist_VelZ.nc','VelZ')


###--- 運動エネルギーの計算とその領域積分 ---###

  ###--- 変数の配列を用意 ---###
  gphysIntKE = [0]

  ###--- 時間について for を回す ---###
  for tn in tnmin..tnmax

    time = tn * dt

    gphysDensBZ = gphys10.cut( 'x'=>x1..x2 ).cut( 'y'=>100 ).cut( 'z'=>z1..z2 )
    gphysVelX   = gphys1.cut( 'x'=>x1..x2 ).cut( 'y'=>100 ).cut( 'z'=>z1..z2 ).cut( 't'=>time )
    gphysVelZ   = gphys2.cut( 'x'=>x1..x2 ).cut( 'y'=>100 ).cut( 'z'=>z1..z2 ).cut( 't'=>time )
#    gphysdisp0 = gphys10.cut( 'x'=>x1..x2 ).cut( 'y'=>100 ).cut( 'z'=>z1..z2 ) #+ gphys2.cut( 't'=>time )

    ###--- 運動エネルギーとその領域積分を計算 ---###
    gphysKE = gphysDensBZ * ( gphysVelX ** 2 + gphysVelZ ** 2 ) * 0.50   # 運動エネルギーを計算
    gphysIntKE[tn] = gphysKE.integrate( 'x' ).integrate( 'z' )           # 運動エネルギーの領域積分を計算

    ###--- テスト出力 ---###
#    p gphysIntKE[tn]

end


###--- 時間微分の計算 ---###

  ###--- 配列を用意 ---###
  gphysDtIntKE = [0]

  ###--- 時間について for を回す ---###
  for tn in tnmin..tnmax-1

    time = tn * dt

    ###--- 各時刻での運動エネルギーの領域積分の時間微分の値を計算 ---###
    gphysDtIntKE[tn] = ( gphysIntKE[tn+1] - gphysIntKE[tn] ) / dt

  end


###--- 結果を出力してみる ---###

  ###--- 時間について for を回す ---###
  for tn in tnmin..tnmax-1

    foo = File.open("result.dat", 'a')    # 出力先のファイルを開く (a -> add)
#    foo.puts var                          # 出力 (改行あり)
    foo.print(tn*dt)                      # 時刻を出力 (改行なし)
    foo.print(",    " )                   # これがないとスペースがあかず, 前の値と次の値の境目がわからない
    foo.print(gphysDtIntKE[tn])           # 運動エネルギーの領域積分の時間微分を出力
    foo.print("\n")                       # 改行
    foo.close                             # 出力先のファイルを閉じる

  end

### 描画してみる ###

#  DCL.sgpset('lcntl', false)
#  GGraph.contour( gphysDtIntKE ) 
#  DCL.grcls