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

#
# = エネルギー方程式の弾性エネルギーの項の領域積分の時間微分を計算するスクリプト
#
# * 履歴 (新しい日付を上に書く)
#   * 2011/10/25  新規作成
#
# * めも
#   * 出力結果の単位が正しくないのは, 重力加速度を定数で安直に与えているから. 
#     * なんとかしたいけど方法がわからないのでひとまず放置. 
#   * 手順
#     * 各時刻での Term の領域積分の値を求める
#       * とある時刻での Term を計算
#       * GPhys のコマンド (.integrate) を使って領域積分
#     * 各時刻での Term の領域積分の時間微分を求める 
#     * 結果を result.dat に出力
#   # これで正しいのかはわからない
#
##########################################################################

require "numru/ggraph"
require "narray"
include NumRu

###--- 定数とか用意 ---###

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

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

Grav = 9.800000000000001    # 重力加速度 (arare5 で使っている値)
CpDry    = 1039.642343614574  # 乾燥定圧比熱 (arare5 で使っている値)l


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

gphys10 = GPhys::IO.open('../../netCDF/thermal-moist_restart.nc','DensBZ')
gphys20 = GPhys::IO.open('../../netCDF/thermal-moist_restart.nc','PTempBZ')
gphys30 = GPhys::IO.open('../../netCDF/thermal-moist_restart.nc','VelSoundBZ')
gphys4  = GPhys::IO.open('../../netCDF/thermal-moist_ExnerAll.nc','ExnerAll')




###--- Term の計算とその領域積分 ---###

  ###--- 変数の配列を用意 ---###
  gphysIntElasE = [0]    # Term の領域積分値

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

    time = tn * dt
    p time

    gphysDensBZ   = gphys10.cut( 'x'=>x1..x2 ).cut( 'y'=>100 ).cut( 'z'=>z1..z2 )
    gphysPTempBZ  = gphys20.cut( 'x'=>x1..x2 ).cut( 'y'=>100 ).cut( 'z'=>z1..z2 )
    gphysVelSoundBZ  = gphys30.cut( 'x'=>x1..x2 ).cut( 'y'=>100 ).cut( 'z'=>z1..z2 )
    gphysExnerAll = gphys4.cut( 'x'=>x1..x2 ).cut( 'y'=>100 ).cut( 'z'=>z1..z2 ).cut( 't'=>time )

    ###--- 弾性エネルギーとその領域積分を計算 ---###
    gphysElasE = 0.5 * gphysDensBZ * ( CpDry * gphysPTempBZ * gphysExnerAll / gphysVelSoundBZ ) ** 2    # 弾性エネルギーを計算
    gphysIntElasE[tn] = gphysElasE.integrate( 'x' ).integrate( 'z' )                               # 弾性エネルギーの領域積分を計算

    ###--- テスト出力 ---###
#    p gphysIntElasE
    
  end


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

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

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

    time = tn * dt

    ###--- 各時刻での Term の領域積分の時間微分の値を計算 ---###
    gphysDtIntElasE[tn] = ( gphysIntElasE[tn+1] - gphysIntElasE[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(gphysDtIntElasE[tn])        # 弾性エネルギーの領域積分の時間微分を出力
    foo.print("\n")                       # 改行
    foo.close                             # 出力先のファイルを閉じる

  end