#!/usr/bin/env ruby

=begin

$  for file in *.xwd ; do xwdtopnm ${file} |  pnmcut 2 2 904 654 > test.pnm; ppmtogif test.pnm > ${file}.gif; done
$ gifsicle --delay=50 --no-loop dcl_*.gif > movie.gif

=end

# ----------------------------------------------
# プログラムメイン

END{

  num = 6
  set_filename(num)           # ファイル名等基本設定
  mradl_anim_main             # gif 作成
  print "gif convert start #{Time.now}\n"
  mkgifanim                   # gif アニメ作成シェルスクリプト

}

# ----------------------------------------------
# 必要なライブラリ, モジュールの読み込み

$local_path = '/home/yukiko/work/ape/yukiko/lib'
$: << $local_path

require "getopts"
require "numru/netcdf"
module NumRu
  class NetCDFVar
    alias put scaled_put
    alias get scaled_get
  end
end

require "numru/ggraph"
require "colorbar"
include NumRu
include NMath


#----------------------------------------------
# 設定

def set_filename(num)

  numnum = 0
  num = ARGV[numnum+1].to_i  if numnum = ARGV.index("-num")

  $exp = "T159L48_eml"
  $timestep = "2hour"

  $latitude = 0.4
  $lat_label = 0 # if $latitude == 0.4

  $fig_path = "/home/yukiko/work/ape/GrADS/T159L48_eml_expID01/hourly/fig/"
  $dir_path = "/home/yukiko/work/ape/GrADS/T159L48_eml_expID01/hourly/netcdf/"

  varlist = ["rain","q","t"]
  set_varinfo
  $var_array.push($rain_hash) 
  $var_array.push($q_hash) 
  $var_array.push($t_hash) 

  $gifname = 
    "#{$exp}_control_#{varlist[0]}-#{varlist[1]}-#{varlist[2]}_y#{$lat_label}"

  print "#{$gifname} name start #{Time.now}\n"

end

def set_varinfo

  rmiss = DCL.glpget('rmiss')
  $var_array = []
  
  $ncfilename =  [
    "#{$dir_path}#{$exp}_expID01_0001.nc"
#    "#{$dir_path}#{$exp}_expID01_0002.nc"
#    "#{$dir_path}#{$exp}_expID01_0003.nc",
#    "#{$dir_path}#{$exp}_expID01_0004.nc",
#    "#{$dir_path}#{$exp}_expID01_0005.nc",
#    "#{$dir_path}#{$exp}_expID01_0006.nc"
  ]

  # こうすい
  $rain_hash = {
    "variable" => "tr_tppn", 
    "diff"     => true, 
    "levels" => NArray[10.0**(-6), 10.0**(-5.5), 10.0**(-5), 10.0**(-4.5), 10.0**(-4), 10.0**(-3.5), DCL.glpget('rmiss')], 
#    "levels" => NArray[0.0001, 0.0002, 0.0003, 0.0004, 0.0005, 0.0006, DCL.glpget('rmiss')], 
    "patterns" => NArray[35999, 45999, 55999, 70999, 75999, 85999], 
    "lost_axis" => "(diff) from (mean) zonal",
    "title" => "precipitation flux (kg m-2 s-1)"
  }

  # 温度偏差
  $t_hash = {
    "variable" => "T", 
    "diff"     => true, 
    "levels" => NArray[rmiss,-1,-0.5, 0, 0.5, 1, 1.5, rmiss], 
    "lost_axis" => "(diff) from (mean) zonal",
    "title" => "temperature (K)"
  }
  
  # 比湿偏差
  $q_hash = {
    "variable" => "Q", 
    "diff"     => true, 
    "levels" => NArray[rmiss,-2e-3,-1e-3, 0, 1e-3, 2e-3, 3e-3,rmiss], 
    "lost_axis" => "(diff) from (mean) zonal", 
    "title" => "specific humidity (1)"
  }

  # 鉛直流
  $sigdot_hash = {
    "variable" => "OMG", 
    "diff"     => false, 
    "levels" => NArray[rmiss,-10e-7,-5e-7, 0, 5e-7, 10e-7, 15e-7, rmiss], 
    "lost_axis" => " ",
    "title" => "sigma-vertical velocity (s-1)"
  }

  # ジオポテンシャル
  $gph_hash = {
    "variable" => "Z", 
    "diff"     => true, 
    "levels" => NArray[rmiss,-10,-5, 0, 5, 10, 15, rmiss], 
    "lost_axis" => "(diff) from (mean) zonal",
    "title" => "geopotential height (m)"
  }

  # 東西風
  $u_hash = {
    "variable" => "U", 
    "diff"     => true, 
    "levels" => NArray[rmiss,-5,-2.5, 0, 2.5, 5, 7.5, rmiss], 
    "lost_axis" => "(diff) from (mean) zonal",
    "title" => "eastward wind (m s-1)"
  }

  # 南北風
  $v_hash = {
    "variable" => "V", 
    "diff"     => true, 
    "levels" => NArray[rmiss,-4,-2, 0, 2, 4, 6, rmiss], 
    "lost_axis" => "(diff) from (mean) zonal",
    "title" => "northward wind (m s-1)"
  }

end


# ----------------------------------------------
# お絵描きメイン

def mradl_anim_main

  # ------------------------------------------------------------

  # こうすい
  gphys_1 = GPhys::NetCDF_IO.open("#{$dir_path}#{$exp}_expID01_0001.nc", 
                                  "tr_tppn")

  gphys_2 = GPhys::NetCDF_IO.open("#{$dir_path}#{$exp}_expID01_0001.nc", 
                                  "Q")

  gphys_3 = GPhys::NetCDF_IO.open("#{$dir_path}#{$exp}_expID01_0001.nc", 
                                  "T")
 
  # ------------------------------------------------------------
  DCL.swpset('LDUMP', true) 
  DCL.swpset('LWAIT',false) ; DCL.swpset('LWAIT0',false) 
  DCL.swpset('LWAIT1',false)
  DCL.swpset('IPOSX', 50) ; DCL.swpset('IPOSY',50)  

  DCL.gropn(1)

  DCL.uzfact(0.6)                # 文字の大きさ倍
  DCL.sgpset('lcorner',false)    # コーナーマークはつけない
  DCL.slmgn(0.0, 0.0, 0.0, 0.0)  # 描画マージン指定
  DCL.sgpset('lfprop',true)      # プロポーショナルフォントを使う
  DCL.sgpset('lfull',true)       # 全画面表示
  DCL.sgpset('lcntl', false)     # アンダーバーは普通に表示
  DCL.uscset('cyspos', 'B' )      # y 軸単位を書く位置

  ipara = NumRu::DCL.uzpget('RSIZEC1')
  DCL.uzpset('RSIZEC0',0.012)    # 小タイトル 文字大きさ
  DCL.uzpset('INDEXL0',1)       # 小タイトル 文字ラインインデックス
  rmiss = DCL.glpget('rmiss')
  colorbar_hash = {
    'vlength'=>0.13,
    "vwidth"=>0.01,
    'labelintv'=>1, 
    'voff'=>0.02,
    "constwidth" => true
  }
  patterns = NArray[30999, 35999, 45999, 55999, 70999, 75999, 85999]

    
#    $file_label = "lat=#{$latitude}, day=#{format("%.1f",(num*1.0*timestep/24) )}"

# 1 枚目: 凝結加熱
  $file_label = "lat=#{$latitude}"
  DCL.uzpset('labelxb',false)
  GGraph.set_fig('viewport'=>[0.10,0.80,0.52,0.65],'new_frame'=>true) 
  GGraph.next_axes('xunits'=>'','yunits'=>'','xtitle'=>'') 

  tone_hash = { 'patterns'=> $var_array[0]["patterns"], 
    'levels'=>$var_array[0]["levels"] }


  g = gphys_make(gphys_1,$latitude)
  g = g.set_att("long_name","").set_att("unit","").set_lost_axes("")

  GGraph.tone( g, true , tone_hash ) 

  DCL::uxsttl("t", $var_array[0]["title"], -1.0 ) 
  DCL::uzinit 
#  DCL.uxpttl("t", 0, $var_array[0]["lost_axis"], 1.0)
  DCL::uxsttl("t", "#{$file_label}" , 1.0 )
  GGraph.color_bar( colorbar_hash )
  DCL::uzinit 

# 2 枚目: 温度偏差
  sigma = 0.83
  $file_label = "lat=#{$latitude}, sigma=#{sigma}"
  DCL.uzpset('labelxb',false)
  GGraph.next_axes('xunits'=>'','yunits'=>'','xtitle'=>nil) 
  GGraph.set_fig('viewport'=>[0.10,0.80,0.34,0.47], 'new_frame'=>false )
  
  tone_hash = { 'patterns'=> patterns, 'levels'=>$var_array[1]["levels"] }


  g = gphys_make(gphys_2,sigma)
  g_tmp = g.
    mean(0).set_lost_axes("")
  g = g.set_att("long_name","").set_att("unit","").set_lost_axes("")

  GGraph.tone( g-g_tmp, true , tone_hash ) 

  DCL::uxsttl("t", $var_array[1]["title"], -1.0 ) 
  DCL::uzinit 
  DCL::uxsttl("t", "#{$file_label}" , 1.0 )
#  DCL.uxpttl("t", 0, $var_array[1]["lost_axis"], 1.0)
  GGraph.color_bar( colorbar_hash )
  DCL::uzinit 
  
# 3 枚目: 比湿偏差
  sigma = 0.22
  $file_label = "lat=#{$latitude}, sigma=#{sigma}"
  DCL.uzpset('labelxb',true)
  GGraph.set_fig('viewport'=>[0.10,0.80,0.16,0.29], 'new_frame'=>false )
  tone_hash = { 'patterns'=> patterns, 'levels'=>$var_array[2]["levels"] }
  

  g = gphys_make(gphys_3,sigma)
  g_tmp = g.
    mean(0).set_lost_axes("")
  g = g.set_att("long_name","").set_att("unit","").set_lost_axes("")

  GGraph.tone( g-g_tmp, true , tone_hash ) 

  DCL::uzinit 
  DCL::uxsttl("t", $var_array[2]["title"], -1.0 ) 
  DCL::uzinit 
  DCL::uxsttl("t", "#{$file_label}" , 1.0 )
#  DCL.uxpttl("t", 0, $var_array[2]["lost_axis"], 1.0)
  GGraph.color_bar( colorbar_hash )
  
  DCL::grcls 

end


# gif アニメ作成シェルスクリプト
def mkgifanim

  `for file in *.xwd ; do xwdtopnm ${file} |  pnmcut 2 2 904 654 > tmp.pnm; ppmtogif tmp.pnm > ${file}.gif; done >& /dev/null `
  `mv dcl_001.xwd.gif #{$fig_path}#{$gifname}.gif`
  `rm dcl_* tmp.pnm`

end



# ----------------------------------------------
# GPhys メソッド追加定義

def GGraph.annotate(str_ary) # GGraph の annotate メソッド再定義

  charsize = 0.7 * DCL.uzpget('rsizec1')
  dvx = 0.01
  dvy = charsize*1.5
  raise TypeError,"Array expected" if ! str_ary.is_a?(Array)
  vxmin,vxmax,vymin,vymax = DCL.sgqvpt
  vx = 0.66 
#  vy = 0.12 - charsize/2
#  vy = 0.10 - charsize/2
  vy = 0.11 - charsize/2
  str_ary.each{|str|
    DCL::sgtxzv(vx,vy,"[ #{str} ]",charsize,0,-1,1)
    vy -= dvy
  }
  nil
end

class GPhys

  # 座標, データ操作ログ記録関数の定義 (Grid.lost_axes 参照)
  def add_lost_axes( lost )
    GPhys.new( @grid.copy.add_lost_axes( lost.to_a ), @data)
  end
  def set_lost_axes( lost )
    GPhys.new( @grid.copy.set_lost_axes( lost.to_a ), @data)
  end

  # rename 再定義 (gphys class の値を返す)
  def rename(name)
    GPhys.new(@grid,@data.copy.rename(name))
  end

  # attribute 変更再定義 (gphys class の値を返す)
  def set_att(name,val)
    GPhys.new(@grid,@data.copy.set_att(name,val))
  end

end

def gphys_make(g,dimcut)

  grid0size = g.grid_copy.axis(0).pos.val.size
  grid_0    = g.grid_copy.axis(0)
  timesize = 5*24*6
  grid_time =
      Axis.new().
      set_pos(VArray.new(NArray.sfloat(timesize).indgen!/24).rename("time").
    put_att("units","days"))

  data =  NArray.sfloat(grid0size,timesize).fill!(0.0)
  grid = Grid.new(grid_0, grid_time)
  
  data[true,0..119] = 
    GPhys::NetCDF_IO.open(
                          "#{$dir_path}#{$exp}_expID01_0001.nc", 
                          g.name).cut(true,dimcut,true).data.val.to_na
  data[true,120..239] = 
    GPhys::NetCDF_IO.open(
                          "#{$dir_path}#{$exp}_expID01_0002.nc", 
                          g.name).cut(true,dimcut,true).data.val.to_na
  data[true,240..359] = 
    GPhys::NetCDF_IO.open(
                          "#{$dir_path}#{$exp}_expID01_0003.nc", 
                          g.name).cut(true,dimcut,true).data.val.to_na
  data[true,360..479] = 
    GPhys::NetCDF_IO.open(
                          "#{$dir_path}#{$exp}_expID01_0004.nc", 
                          g.name).cut(true,dimcut,true).data.val.to_na
  data[true,480..599] = 
    GPhys::NetCDF_IO.open(
                          "#{$dir_path}#{$exp}_expID01_0005.nc", 
                          g.name).cut(true,dimcut,true).data.val.to_na
  data[true,600..719] = 
    GPhys::NetCDF_IO.open(
                          "#{$dir_path}#{$exp}_expID01_0006.nc", 
                          g.name).cut(true,dimcut,true).data.val.to_na

  data = VArray.new(data).rename("tr_tppn")
  gphys = GPhys.new(grid, data) 
  return gphys

end