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


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

def set_filename

  $lat_label = 0  if $latitude == - 1.5
  $lat_label = 0  if $latitude == - 0.4

  if $timestep == "1hour"
    $gifname = "#{$exp}_control_#{$variable}_y#{$lat_label}-#{$timestep}"
  else
    $gifname = "#{$exp}_control_#{$variable}_y#{$lat_label}"
  end


  $fig_path = "/home/yukiko/work/ape/GrADS/fig/"
  $dir_path = "/home/yukiko/work/ape/GrADS/#{$exp}_expID01/hourly/"


  rmiss = DCL.glpget('rmiss')
  if $variable == "T"       # 温度偏差
    $filename = ["#{$dir_path}/0001/T.ctl", "#{$dir_path}/0002/T.ctl"]
    $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 == "DTCUM" # 凝結加熱
    $filename = ["#{$dir_path}/0001/DTCUM.ctl", "#{$dir_path}/0002/DTCUM.ctl"]
#    $levels = NArray[0, 0.4e-4, 0.8e-4, 1.2e-4, 1.6e-4, 2.0e-4, 2.4e-4, rmiss]
    $levels = NArray[rmiss,-1.6e-4,-0.8e-4, 0, 0.8e-4, 1.6e-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}/0001/Q.ctl"  , "#{$dir_path}/0002/Q.ctl"] 
    $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}"
  end

  # -----------------------
  $day_hash = {
    "T39L48_eml"=> 40,
    "T39L48_ias"=> 45,
    "T39L48_ksc"=> 40,
    "T39L48_kuo"=> 15,
    "T39L48_mca"=> 40,
    "T39L48_non"=> 20, 
    "T159L48_non"=> 60,
    "T159L48_eml"=> 25
  }

end


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

def mradl_anim_main

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


  gphys_c1 = 
    GPhys::GrADS_IO.open("#{$dir_path}/0001/PRCPC.ctl", 
                         "PRCPC")
  gphys_l1 = 
    GPhys::GrADS_IO.open("#{$dir_path}/0001/PRCPL.ctl", 
                         "PRCPL")

  gphys_c2 = 
    GPhys::GrADS_IO.open("#{$dir_path}/0002/PRCPC.ctl", 
                         "PRCPC")
  gphys_l2 = 
    GPhys::GrADS_IO.open("#{$dir_path}/0002/PRCPL.ctl", 
                         "PRCPL")

  gphys_c3 = 
    GPhys::GrADS_IO.open("#{$dir_path}/0003/PRCPC.ctl", 
                         "PRCPC")
  gphys_l3 = 
    GPhys::GrADS_IO.open("#{$dir_path}/0003/PRCPL.ctl", 
                         "PRCPL")

  gphys_t1 = GPhys::GrADS_IO.open("#{$dir_path}/0001/#{$variable}.ctl", 
                                  $variable)
  gphys_t2 = GPhys::GrADS_IO.open("#{$dir_path}/0002/#{$variable}.ctl", 
                                  $variable)
  gphys_t3 = GPhys::GrADS_IO.open("#{$dir_path}/0003/#{$variable}.ctl", 
                                  $variable)

  if $variable == "DTCUM" 
    gphys_tt1 = 
      GPhys::GrADS_IO.open("#{$dir_path}/0001/DTLSC.ctl", 
                           "DTLSC")
    gphys_tt2 = 
      GPhys::GrADS_IO.open("#{$dir_path}/0002/DTLSC.ctl", 
                           "DTLSC")
    gphys_tt3 = 
      GPhys::GrADS_IO.open("#{$dir_path}/0003/DTLSC.ctl", 
                           "DTLSC")
  end

  gphys_u1 = GPhys::GrADS_IO.open("#{$dir_path}/0001/U.ctl", "U")
  gphys_u2 = GPhys::GrADS_IO.open("#{$dir_path}/0002/U.ctl", "U")
  gphys_u3 = GPhys::GrADS_IO.open("#{$dir_path}/0003/U.ctl", "U")
  
  gphys_sig1 = GPhys::GrADS_IO.open("#{$dir_path}/0001/OMG.ctl", "OMG")
  gphys_sig2 = GPhys::GrADS_IO.open("#{$dir_path}/0002/OMG.ctl", "OMG")
  gphys_sig3 = GPhys::GrADS_IO.open("#{$dir_path}/0003/OMG.ctl", "OMG")

  # ----------------------------------------------------------
  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

#  (5).times{ |num|

  timeall.times{ |num|

    $file_label = "day=#{format("%.1f",(num*1.0*timestep/24+ $day_hash[$exp]) )}"
    print "  #{$file_label} #{Time.now}\n"

# 1 枚目

    GGraph.next_fig('viewport'=>[0.10,0.80,0.50,0.65]) 
    GGraph.next_axes('xunits'=>'','yunits'=>'','xtitle'=>'') 
    DCL.uzpset('labelxb',false)

    if (120 > num*timestep)
      g = gphys_c1.cut(true,-30..30,true,true)[true,true,0,num*timestep] + 
        gphys_l1.cut(true,-30..30,true,true)[true,true,0,num*timestep]
    elsif (120 <= num*timestep) && (240 > num*timestep)
      g = gphys_c2.cut(true,-30..30,true,true)[true,true,0,num*timestep-120] + 
        gphys_l2.cut(true,-30..30,true,true)[true,true,0,num*timestep-120]
    else
      g = gphys_c3.cut(true,-30..30,true,true)[true,true,0,num*timestep-240] + 
        gphys_l3.cut(true,-30..30,true,true)[true,true,0,num*timestep-240]
    end
    g = g.lon_lotate.
      set_lost_axes($lost_axis). 
      set_att("ape_name","precipitation flux").
      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([-180,180],[$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 }


    if (120 > num*timestep)
      g = gphys_t1.cut(true,$latitude,true,true)[true,true,num*timestep]
      g_tmp = gphys_t1.cut(true,$latitude,true,true).mean(0)[true,num*timestep]

      if $variable == "DTCUM" 
        g = g + gphys_tt1.cut(true,$latitude,true,true)[true,true,num*timestep]
        g_tmp = g_tmp + gphys_tt1.cut(true,$latitude,true,true).
          mean(0)[true,num*timestep]
      end

    elsif (120 <= num*timestep) && (240 > num*timestep)
      g = gphys_t2.cut(true,$latitude,true,true)[true,true,num*timestep-120]
      g_tmp = gphys_t2.cut(true,$latitude,true,true).mean(0)[true,num*timestep-120]
      if $variable == "DTCUM" 
        g = g + gphys_tt2.cut(true,$latitude,true,true)[true,true,num*timestep-120]
        g_tmp = g_tmp + gphys_tt2.cut(true,$latitude,true,true).
          mean(0)[true,num*timestep-120]
      end
    else
      g = gphys_t3.cut(true,$latitude,true,true)[true,true,num*timestep-240]
      g_tmp = gphys_t3.cut(true,$latitude,true,true).mean(0)[true,num*timestep-240]
      if $variable == "DTCUM" 
        g = g + gphys_tt3.cut(true,$latitude,true,true)[true,true,num*timestep-240]
        g_tmp = g_tmp + gphys_tt3.cut(true,$latitude,true,true).
          mean(0)[true,num*timestep-240]
      end
    end

    if $variable == "DTCUM" 
      g = g.lon_lotate. 
        set_lost_axes("").
        set_att("long_name","").set_att("unit","")
    else
      g = g.lon_lotate.
        set_lost_axes($lost_axis).
        set_att("long_name","").set_att("unit","")
    end

    GGraph.next_fig('viewport'=>[0.10,0.80,0.20,0.45], 'new_frame'=>false )
    GGraph.set_axes('ytitle'=>"sigma",'yunits'=>"1")

    if $variable == "DTCUM"
      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
    GGraph.set_axes('ytitle'=>nil,'yunits'=>nil)

    if (120 > num*timestep)
      vect_plot(gphys_u1,gphys_sig1,num,timestep,0) if $vector 
    elsif (120 <= num*timestep) && (240 > num*timestep)
      vect_plot(gphys_u2,gphys_sig2,num,timestep,-120) if $vector 
    else
      vect_plot(gphys_u3,gphys_sig3,num,timestep,-240) if $vector 
    end

    DCL::uzinit 
    DCL::uxsttl("t", $title, -1.0 ) 
    DCL::uzinit 
    DCL::uxsttl("t", "lat=#{$latitude}" , 1.0 )
    $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 --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
  $exp = "T39L48_eml"
  $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")
  $exp = ARGV[num+1]           if num = ARGV.index("-exp")

  print "$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


  # 経度 180 度回転
  def lon_lotate

    # 軸の回転 (180 度引き算)
#    lon = VArray.new(@grid.axis("lon").pos.val-180).rename("lon")
    lon = VArray.new(@grid.axis(0).pos.val-180).rename("lon").
      put_att("units",@grid.axis(0).pos.get_att("units"))
    lon = Axis.new().set_pos(lon)
    @grid = @grid.change_axis(0, lon)

#    lat = VArray.new(@grid.axis(1).pos.val).rename("lat")
#    lat = Axis.new().set_pos(lat)
#    @grid = @grid.change_axis(1, lat)

    # データの回転
    dim = @data.val.shape[0]
    data_val =  @data.copy.val
    data_val[0..dim/2-1, false] = @data.val[dim/2..dim-1, false]
    data_val[dim/2..dim-1, false] = @data.val[0..dim/2-1, false]
    data = @data.copy
    data.val = (data_val)

    # gphys class の値を返す
    GPhys.new( @grid, data )

  end

end

# ----------------------------------------------
def vect_plot(gphys_u,gphys_sig,num,timestep,numdiff)

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

      g_u = g_u.lon_lotate
      g_sig = g_sig.lon_lotate

    # ベクトルの間引
    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_l = 3 * 2*2 if $exp =~ /T159/
    num_z = 2

    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::ugrset('XFACT1', 10E-6*100)
    DCL::ugrset('YFACT1', 10.0/100)
    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..-1] ) 

end