#!/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{

#   $variable = "t"
#  $variable = "q"
#  $variable = "qcnd"

#  $vector = true
#  $latitude = 1.4

  setopt                 # 引数解析
  set_filename           # ファイル名等基本設定
  mradl_anim_main        # gif 作成
  print "gif convert start\n"
  mkgifanim              # gif アニメ作成シェルスクリプト

}


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

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

#require "getopts"
require "optparse"
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

$exp_hash = {
  "con" => "con", 
  "a" => "mradlAa", 
  "b" => "mradlAb", 
  "c" => "mradlAc", 
  "d" => "mradlAd"
}


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

def set_filename

  $lat_label = 0  if $latitude == 1.4

#  if $vector
# $gifname = "#{$cum}_#{$exp}_#{$variable}-vect_lat#{$latitude}-#{$timestep}"
#  else
#    $gifname = "#{$cum}_#{$exp}_#{$variable}_lat#{$latitude}"
#  end
  if $timestep == "1hour"
    $gifname = "#{$cum}_#{$exp}_#{$variable}_y#{$lat_label}-#{$timestep}"
  else
    $gifname = "#{$cum}_#{$exp}_#{$variable}_y#{$lat_label}"
  end

  $fig_path = "/home/yukiko/work/aqua3d/yukiko/mradl/20051213/fig/anim/"
  $dir_path = "/home/yukiko/work/aqua3d/yukiko/mradl/20051213/#{$cum}700#{$exp_hash[$exp]}/"

  rmiss = DCL.glpget('rmiss')
  if $variable == "t"       # 温度偏差
    $filename = ["#{$dir_path}t_out_10.nc", "#{$dir_path}t_out_20.nc"]
    $levels = NArray[rmiss,-1,-0.5, 0, 0.5, 1, 1.5, rmiss]
    $lost_axis = "(diff) from (mean) zonal"
    if $vector
      $title = "temperature, wind (K, (m s-1, s-1))"
    else
      $title = "temperature (K)"
    end
    $comment = "lat = #{$latitude}, (diff) zonal"
  elsif $variable == "qcnd" # 凝結加熱
    $filename = ["#{$dir_path}qcnd_sum_10.nc", "#{$dir_path}qcnd_sum_20.nc"]
    $levels = NArray[0, 0.4e-4, 0.8e-4, 1.2e-4, 1.6e-4, 2.0e-4, 2.4e-4, rmiss]
    $lost_axis = ""
    if $vector
      $title = "condensation heating, wind (K s-1, (m s-1, s-1))"
    else
      $title = "condensation heating (K s-1)"
    end
    $comment = "lat = #{$latitude}"
  elsif $variable == "q"    # 比湿偏差
    $filename = ["#{$dir_path}q_out_10.nc"  , "#{$dir_path}q_out_20.nc"] 
    $levels = NArray[rmiss,-2e-3,-1e-3, 0, 1e-3, 2e-3, 3e-3,rmiss]
    $lost_axis = "(diff) from (mean) zonal"
    if $vector
      $title = "specific humidity, wind (1, (m s-1, s-1))"
    else
      $title = "specific humidity (1)"
    end
    $comment = "lat = #{$latitude}"
#  elsif $variable == "q_sigdot"    # 比湿偏差
#    $filename = ["#{$dir_path}q_sigdot_out_10.nc"  , "#{$dir_path}q_sigdot_out_20.nc"] 
#    $levels = NArray[rmiss,-5e-9,-2.5e-9, 0, 2.5e-9, 5e-9, 7.5e-9, rmiss] 
#    $lost_axis = " "
#    if $vector
#      $title = "q*sigdot, wind (s-1, (m s-1, s-1))"
#    else
#      $title = "q*sigdot (s-1)"
#    end
#    $comment = "lat = #{$latitude}"
  end


end


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

def mradl_anim_main

  # ----------------------------------------------------------
  gphys = 
    GPhys::NetCDF_IO.open(["#{$dir_path}rain_sum_10.nc","#{$dir_path}rain_sum_20.nc"], 
                          "rain") * 0.4*10**(-6)
  gphys_t = GPhys::NetCDF_IO.open($filename, $variable)
  gphys_u = GPhys::NetCDF_IO.open(["#{$dir_path}u_out_10.nc","#{$dir_path}u_out_20.nc"], "u")
  gphys_sig = 
    GPhys::NetCDF_IO.open(["#{$dir_path}sigdot_out_10.nc","#{$dir_path}sigdot_out_20.nc"], "sigdot")

  # ----------------------------------------------------------
  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 軸単位を書く位置

  $colorbar_hash = {
    'vlength'=>0.13,
    "vwidth"=>0.01,
    'labelintv'=>1, 
    'voff'=>0.02, 
    "constwidth" => true
  }


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

  if $timestep == "1hour"
    timeall = 15*24
    timestep = 1
  elsif $timestep == "2hour"
    timeall = 15*12
    timestep = 2
  end

#  (25*24).times{ |num|
#  (15*12).times{ |num|
#  (5).times{ |num|

  timeall.times{ |num|

#    numnum= num * 10/ 24 * timestep
#    $file_label = "day=#{numnum*1.0/10+5}"
#    $file_label = "day=#{num*timestep/24}"
#    $file_label = "day=#{(num*1.0*timestep/24).to_.scanf("%f2.1")}"
    $file_label = "day=#{format("%.1f",(num*1.0*timestep/24) )}"
    print "  #{$file_label} #{Time.now}\n"

# 1 枚目

#    GGraph.next_fig('viewport'=>[0.13,0.78,0.45,0.60]) 
    GGraph.next_fig('viewport'=>[0.10,0.80,0.50,0.65]) 
    GGraph.next_axes('xunits'=>'','yunits'=>'','xtitle'=>'') 
    DCL.uzpset('labelxb',false)
#    g = gphys.cut(true,-30..30,true)[true,true,num*timestep+119].
    g = gphys.cut(true,-30..30,true)[true,true,num*timestep].
      set_lost_axes($lost_axis). 
      set_att("ape_name","#{gphys.data.get_att("long_name")}").
      set_att("long_name","").set_att("unit","")

    patterns = NArray[35999, 45999, 55999, 70999, 75999, 85999]
    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')]
    tone_hash = { 'patterns'=> patterns, 'levels'=>levels }
    GGraph.tone( g, true, tone_hash ) 
    DCL.uzpset('labelxb',true)
    DCL.uulin([0,360],[$latitude,$latitude])

    # タイトル とカラーバー
    DCL::uxsttl("t", "#{$file_label}" , 1.0 )
    DCL::uzinit 
    DCL::uxsttl("t", "precipitation (kg m-2 s-1)", -1.0 ) 
    GGraph.color_bar($colorbar_hash)

# 2 枚目

    patterns = NArray[30999, 35999, 45999, 55999, 70999, 75999, 85999]
    tone_hash = { 'patterns'=> patterns, 'levels'=>$levels }

#    g = gphys_t.cut(true,$latitude,true,true)[true,true,num*timestep+119].
    g = gphys_t.cut(true,$latitude,true,true)[true,true,num*timestep].
      set_lost_axes($lost_axis). 
      set_att("long_name","").set_att("unit","")
#    g_tmp = gphys_t.cut(true,$latitude,true,true).mean(0)[true,num*timestep+119].
    g_tmp = gphys_t.cut(true,$latitude,true,true).mean(0)[true,num*timestep].
      set_lost_axes($lost_axis).
      set_att("long_name","").set_att("unit","")

#    GGraph.next_fig('viewport'=>[0.13,0.78,0.15,0.40], 'new_frame'=>false )
    GGraph.next_fig('viewport'=>[0.10,0.80,0.20,0.45], 'new_frame'=>false )

    if $variable == "qcnd"
      GGraph.tone( g, true , tone_hash ) 
    elsif $variable == "q_sigdot"
      GGraph.tone( -g, true , tone_hash ) 
    else
      GGraph.tone( g-g_tmp, true , tone_hash ) 
    end

    vect_plot(gphys_u,gphys_sig,num,timestep) if $vector 
    DCL::uzinit 
    DCL::uxsttl("t", $title, -1.0 ) 
    DCL::uzinit 
    DCL::uxsttl("t", "lat=#{$latitude}" , 1.0 )
#    $colorbar_hash.store("vcent", 0.30)
    $colorbar_hash.store("vcent", 0.35)
    GGraph.color_bar($colorbar_hash )
    $colorbar_hash.delete("vcent")    

  }
  
  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 `
#  `gifsicle --delay=50 --loopcount=100 dcl_*.gif > #{$gifname}_anim.gif`
  `gifsicle --loopcount=100 dcl_*.gif > #{$fig_path}#{$gifname}_anim.gif`
  `mv dcl_001.xwd.gif #{$fig_path}#{$gifname}_t0.gif`
  `rm dcl_* tmp.pnm`

end

# ----------------------------------------------
# 引数解析
def setopt

  num = 0

  $vector = true
  $cum = "adj"
  $exp = "con"
  $timestep = "2hour"

  $variable = ARGV[num+1]      if num = ARGV.index("-var")
  $latitude = ARGV[num+1].to_f if num = ARGV.index("-lat")
  $timestep = "1hour"          if num = ARGV.index("-1hour")
  $cum = "kuo"                 if num = ARGV.index("-kuo")
  $exp = ARGV[num+1]           if num = ARGV.index("-exp")
#  $vector   = true             if num = ARGV.index("-vect")
#  $timestep = "2hour"          if num = ARGV.index("-2hour")

  print "$cum = #{$cum}, $exp = #{$exp}, $variable = #{$variable}, $vector = #{$vector}, $latitude = #{$latitude}, $timestep = #{$timestep}\n"

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.08 - charsize/2
  vy = 0.13 - 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 vect_plot(gphys_u,gphys_sig,num,timestep)

#    g_u =  gphys_u.cut(true,$latitude,true,true)[true,true,num*timestep+119].
    g_u =  gphys_u.cut(true,$latitude,true,true)[true,true,num*timestep].
      set_lost_axes($lost_axis)
#    g_sig = gphys_sig.cut(true,$latitude,true,true)[true,true,num*timestep+119].
    g_sig = gphys_sig.cut(true,$latitude,true,true)[true,true,num*timestep].
      set_lost_axes($lost_axis)

    # ベクトルの間引
    u_data  = NArrayMiss.sfloat( g_u.axis(0).pos.val.size , 
			   g_u.axis(1).pos.val.size ).fill!(0.0)
 
    sig_data = NArrayMiss.sfloat( g_sig.axis(0).pos.val.size , 
			    g_sig.axis(1).pos.val.size ).fill!(0.0)

    num_l = 3 
    num_z = 1

    g_u.axis(0).pos.val.size.times{ |num|
      if (num % num_l) == 0
	u_data[num,true] =  g_u.data.val[num,true]
	sig_data[num,true] =  g_sig.data.val[num,true]
      elsif
	u_data[num,true]  = 0.0
	sig_data[num,true]  =  0.0
      end
    }
    
    g_u.axis(1).pos.val.size.times{ |num|
      if (num % num_z) == 0
	u_data[true,num] = u_data[true,num]
	sig_data[true,num] = sig_data[true,num]
      elsif
	u_data[true,num] = 0.0
	sig_data[true,num] = 0.0
      end
    }

    #    (u,-sigdot) 
    DCL::uglset('LNRMAL', false)
    DCL::ugrset('XFACT1', 10E-6*100)
    DCL::ugrset('YFACT1', 15.0*1000)
    DCL::uglset('LUNIT', true)
    DCL::ugrset('VXUNIT', 0.04)
    DCL::ugrset('VYUNIT', 0.04)
    DCL::ugrset('VXUOFF', 0.02)

    DCL::ugvect( u_data, -sig_data[true,0..-2] ) 

end