#!/usr/bin/env ruby
# ----------------------------------------------
# local load path

# $local_path = '/work11/ape/yukiko/lib'
# $local_path = '/home/yukiko/tmp/ape-data/lib'
 $local_path = '/home/yukiko/work/ape/yukiko/lib'
$: << $local_path

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

load "#{$local_path}/ape-view.rb"

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

def set_path

#  $cum_parameta = "kuo"
#  $exp = "mradlAd"
#  $fig_path = "/work11/ape/yukiko/figs/tmp/"
  $fig_path = "/home/yukiko/work/yukiko/eva01-home/work/aqua-mradl-ronbun/figs/tmp/"
  $ncfile_path = "/home/yukiko/work/aqua3d/yukiko/mradl/data/wk-plot-data/"
#  $file_label = "AGCM5_#{$cum}_#{$expid}"
  $groupid = "AGCM5-#{$cum}-#{$expid}"
  $file_label = "#{$cum_parameta}_#{$exp_lavel[$exp]}"
end

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

host = "eva01"
 a = Ape_mkfig.new 3
# a = Ape_mkfig.new 2
# a = Ape_mkfig.new 1

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

END{

#  $exp = "mradlAa"
#  $exp = "con"
#  $exp = "con-east"
#  $exp = "con-west"
#  $cum_parameta = "kuo"
#  $cum_parameta = "adj"

#  set_path
#  a.gr_comp_point
#  a.gr_comp_tuom

  $cum_parameta = "kuo"
  $exp = "mradlAa"

 ["mradlAa", "con","mradlAb","mradlAc","mradlAd"].each{ |exp|
   $exp = exp
   set_path
   a.gr_comp_point
   a.gr_comp_tuom
  }

  $cum_parameta = "adj"
  $exp = "con-east"
  set_path
  a.gr_comp_point
  a.gr_comp_tuom

  $exp = "con-west"
  set_path
  a.gr_comp_point
  a.gr_comp_tuom


}

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

  $exp_lavel ={ 
    "con" => "control", 
    "mradlAa" => "a",
    "mradlAb" => "b",
    "mradlAc" => "c",
    "mradlAd" => "d",
#    "mradlAa" => "mradlAa",
#    "mradlAb" => "mradlAb",
#    "mradlAc" => "mradlAc",
#    "mradlAd" => "mradlAd",
    "con-east" => "control_eastward",
    "con-west" => "control_westward"
  }

  $exp_file_lavel ={ 
    "con" => "con", 
    "mradlAa" => "mradlAa",
    "mradlAb" => "mradlAb",
    "mradlAc" => "mradlAc",
    "mradlAd" => "mradlAd",
    "con-east" => "con",
    "con-west" => "con"
  }

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


class Ape_mkfig

  def mkfig_plot_comp(gphys, gphys_comp)

    if @fig_num == 3 then
      @data.mkdump_comp(gphys, gphys_comp)
    elsif @fig_num == 2 then
      @data.mkps_comp(gphys, gphys_comp)
    else
      @data.plot_comp(gphys, gphys_comp)
    end

  end

  def gr_comp_point

    lost_axes = ["latitude=1.395307"]
    # t_conv composite
    file_path = "/home/yukiko/work/yukiko/eva01-home/work/aqua-mradl-ronbun/composite_nc/#{$cum_parameta}700#{$exp}_composite_tuom_y32.nc"
    @data = Ape.new(file_path)
    comp_rain = @data.gphys_open("comp_point")

    # tr_tppn
    file_path = "/home/yukiko/work/aqua3d/yukiko/mradl/data/wk-plot-data/#{$cum_parameta}700#{$exp_file_lavel[$exp]}400/rain_sum.nc"
    @data = Ape.new(file_path)
    gphys = @data.gphys_open("rain").rename("tr_tppn").
      set_att("ape_name","precipitation flux"). 
      set_att("units","kg m-2 s-1")*0.4*10**(-6)
    gphys = gphys.cut(true,0,true)[true,-400..-1]

    grid_1 =
      Axis.new().
      set_pos(VArray.new(NArray.sfloat(400).indgen!/4).rename("time").
            put_att("units","days"))

    grid_0 = comp_rain.grid_copy.axis(0)
    grid = Grid.new(grid_0,grid_1)
    comp_rain = GPhys.new(grid, comp_rain.data ).set_lost_axes(lost_axes)

    grid_0 = gphys.grid_copy.axis(0)
    grid = Grid.new(grid_0,grid_1)
    gphys = GPhys.new(grid, gphys.data ).set_lost_axes(lost_axes)

    # 描画
#    mkfig_plot(comp)
    mkfig_plot_comp(gphys, comp_rain)

  end


  def gr_comp_tuom

    file_path = "/home/yukiko/work/yukiko/eva01-home/work/aqua-mradl-ronbun/composite_nc/#{$cum_parameta}700#{$exp}_composite_tuom_y32.nc"

    @data = Ape.new(file_path)
    lost_axes = ["latitude=1.395307"]

    # T composite
    t_comp = @data.gphys_open("t_anm")
    t_comp = t_comp.
      set_att("ape_name", "t (u, -sigdot) composite").
      set_att("units", "K, (m s-1, s-1)" ).
      set_lost_axes(lost_axes).
      add_lost_axes("(diff) from (mean) zonal").
      rename("agcm_comp_tuom")

    # U composite
    u_comp = @data.gphys_open("u_anm")

    # sigdot composite
    om_comp = @data.gphys_open("sigdot_anm")

    # ベクトルの間引
    u_data  = NArray.sfloat( u_comp.axis(0).pos.val.size , 
			   u_comp.axis(1).pos.val.size ).fill!(0.0)
 
    om_data = NArray.sfloat( om_comp.axis(0).pos.val.size , 
			    om_comp.axis(1).pos.val.size ).fill!(0.0)

    num_l = u_comp.axis(0).pos.val.size/128
#    num_l = u_comp.axis(0).pos.val.size/128*3
    num_z = u_comp.axis(1).pos.val.size/16
    
    u_comp.axis(0).pos.val.size.times{ |num|
      if (num % num_l) == 0
	u_data[num,true]  =  u_comp.data.val[num,true]
	om_data[num,true] =  om_comp.data.val[num,true]
      elsif
	u_data[num,true]  = 0.0
	om_data[num,true] = 0.0
      end
    }
    
    num_z =1 if num_z < 1

    u_comp.axis(1).pos.val.size.times{ |num|
      if (num % num_z) == 0
	u_data[true,num]  =  u_data[true,num]
	om_data[true,num] =  om_data[true,num]
      elsif
	u_data[true,num]  = 0.0
	om_data[true,num] = 0.0
      end
    }

    # 軸の long_name 直し (pressure level => sigma)
    z = t_comp.grid_copy.axis(1).pos.
      set_att("long_name","sigma").
      set_att("units","1")
    z = Axis.new().set_pos(z)
    grid = t_comp.grid_copy.change_axis(1, z)

    # 配列数取得
    dim = t_comp.axis(0).pos.val.size


=begin
    # 描画 (T)
    t_comp = GPhys.new( grid, t_comp.data )
    mkfig_plot(t_comp, u_data, -om_data)
=end
    # 描画 (T)
    t_comp = GPhys.new( grid, t_comp.data )
    mkfig_plot(t_comp[(dim/3 +1)..(dim*2/3 -1),true],
	       u_data[(dim/3 +1)..(dim*2/3 -1),true], 
	       -om_data[(dim/3 +1)..(dim*2/3 -1),true])

  end


end

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

class Ape

  def plot_comp(gphys, gphys_u=nil, gphys_v=nil)
    gropn(1) if $gropn == nil
    
    if gphys.class == Array
      plot_set(gphys[0]) 
    else 
      plot_set(gphys) 
    end

    plot_main_comp(gphys, gphys_u)

    # タイトルを標準出力へ
    if gphys.class == Array then
      if  gphys[0].data.get_att("ape_name") then
	puts "#{gphys[0].name} (#{gphys[0].data.get_att("ape_name")})"
      end
    elsif gphys.data.get_att("ape_name") then
      puts "#{gphys.name} (#{gphys.data.get_att("ape_name")})"
    else
      puts "#{gphys.name}"
    end
  end

  def mkdump_comp(gphys, gphys_u=nil, gphys_v=nil)
    gropn(3) if $gropn == nil

    plot_comp(gphys, gphys_u)
    
    file_name = "#{$file_label.gsub("control_eastward","con_east").gsub("control_westward","con_west").gsub("control","con")}_#{gphys_u.name}.gif"
    `rm dcl_* tmp.pnm`  
    plot(gphys_u)
    `for file in *.xwd;do xwdtopnm ${file} | pnmcut 2 2 904 654 > tmp.pnm;done`
#    `for file in *.xwd;do xwdtopnm ${file} | pnmcut 2 2 840 630 > tmp.pnm;done`
    `ppmtogif tmp.pnm > #{$fig_path}#{file_name}`
    `rm dcl_* tmp.pnm`  

    $file_name = $pre_file

  end

  def mkps_comp(gphys, gphys_u=nil, gphys_v=nil)
    'rm dcl*.ps'
    gropn(2) if $gropn == nil

    plot_comp(gphys, gphys_u)

    file_name = "#{$file_label}_#{gphys_u.name}.ps"
    gphys = gphys[0]  if gphys.class == Array
#    file_name = "#{$file_label}_#{file_name}.ps"
    `mv dcl*.ps  #{$fig_path}#{file_name}` 

  end

  def plot_main_comp(gphys, gphys_u=nil)

    $colorbar_hash = {
      'vlength'=>0.2,
      "vwidth"=>0.015,
      'labelintv'=>1,
      'voff'=>0.02,
      'vcent'=> 0.50, 
      'charfact' => 1.0,
      "constwidth" => true
    }

    # attribute の long_name を消す (タイトル描画位置の都合)
    if gphys.data.get_att("long_name")
      gphys = gphys.set_att("long_name","")
    end

    GGraph.set_fig($fig_set_hash)

    # 描画
    if gphys.axnames.size == 1 || gphys.name =~ /gt_/
      
      GGraph.line( gphys, true, $line_hash )
    elsif  $tone_flag == false
      GGraph.contour( gphys, true, $cont_hash )
    elsif $tone_hash['levels'][0] == $tone_hash['levels'][1] 
      GGraph.tone( gphys, true ) 
      GGraph.contour( gphys, false ) unless $cont_flag == false
    else
      GGraph.tone( gphys, true , $tone_hash ) 
      GGraph.contour( gphys, false, $cont_hash ) unless $cont_flag == false
    end

#    plot_set(gphys_u) 
#    GGraph.contour( gphys_u, false, $cont_hash )

    # comp_point のマーカー
    gphys_u.grid_copy.axis(1).pos.val.size.times{ |num|

      marker = NArray.sfloat(gphys_u.grid_copy.axis(0).pos.val.size).fill!(-100.0)
      gphys_u.grid_copy.axis(0).pos.val.size.times{ |numnum|

	if gphys_u.data.val[numnum,num] == 1.0
	  marker[numnum] = gphys_u.grid_copy.axis(1).pos.val[num]
	  puts num, gphys_u.grid_copy.axis(1).pos.val[numnum]
	end
      }

      DCL::uusmki(2)
      DCL::uusmkt(5)
#      DCL::uusmks(2)
      DCL::uumrk(gphys_u.grid_copy.axis(0).pos.val, marker)
      print "num = #{num}\n"

    }


    # タイトル
    DCL::uzinit  
    tani = "(#{gphys.data.get_att("units")})"
    DCL::uxsttl("t", tani , -1.0 )
    if  gphys.data.get_att("ape_name") then
      DCL::uxmttl("t", gphys.data.get_att("ape_name").gsub("_", " "), -1.0 ) 
    end

    $file_label = "#{@file}@#{gphys.name}"  if $file_label == "filename"

    # nc ファイル名
    if $gropn == 2
      charsize = 0.79 - $file_label.size * 0.0115 
      DCL::sgtxzv(charsize,0.62,$file_label,0.8 * DCL.uzpget('rsizec1'),0,-1,1)
    else
#      charsize = 0.81 - $file_label.size * 0.0095 
#      DCL::sgtxzv(charsize,0.6,$file_label,0.8 * DCL.uzpget('rsizec1'),0,-1,1) 
      charsize = 0.78 - $file_label.size * 0.0095 
      DCL::sgtxzv(charsize,0.65,$file_label,0.8 * DCL.uzpget('rsizec1'),0,-1,1) 
    end

#    $file_label = "filename" if $file_label == "#{@file}@#{gphys.name}"  

    # トーンバー
#    plot_set(gphys) 
#    DCL::Util::color_bar($cbar_conf)  
    GGraph.color_bar($colorbar_hash)

  end

end

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


class Ape

  def mkdump(gphys, gphys_u=nil, gphys_v=nil)
    gropn(3) if $gropn == nil

    if gphys_u == nil || gphys_v == nil then
      plot(gphys)
    else
      plot(gphys, gphys_u, gphys_v)
    end
    
    gphys = gphys[0]     if gphys.class == Array
    
    file_name = "#{$file_label.gsub("control_eastward","con_east").gsub("control_westward","con_west").gsub("control","con")}_#{gphys.name}.gif"
    `rm dcl_* tmp.pnm`  
    plot(gphys)
    `for file in *.xwd;do xwdtopnm ${file} | pnmcut 2 2 904 654 > tmp.pnm;done`
#    `for file in *.xwd;do xwdtopnm ${file} | pnmcut 2 2 840 630 > tmp.pnm;done`
    `ppmtogif tmp.pnm > #{$fig_path}#{file_name}`
    `rm dcl_* tmp.pnm`  

    $file_name = $pre_file

  end

  # トーン, コンターパターンの設定 (変数固有)
  def plot_set(gphys)

    # コンターを描かない条件
    if gphys.name =~ /tr_/ || gphys.name =~ /horinout/  
      $cont_flag = false 
    else
      $cont_flag = true
    end

    # トーンを描かない条件
    if gphys.name == "comp_point"
      $tone_flag = false 
    else
      $tone_flag = true
    end

    # 初期化
    $fig_set_hash = {"window" => nil}
    $line_hash = { 'exchange'=>false, 'index'=>3 }


    if gphys.name == "tr_tppn" 
      # 塚原カラー
      patterns = NArray[35999, 45999, 55999, 70999, 75999, 85999]
# APE 
#      levels = NArray[0.0001, 0.0002, 0.0003, 0.0004, 0.0005, 0.0006, 1000]
# log_plot
#      levels = NArray[0.000001,0.000005,0.00001,0.00005, 0.0001, 0.0005, 1000]
      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 }
      $cont_hash = {'int'=>0.5}
      $cbar_conf = {
	"levels"=>levels, 
	"colors"=>patterns,
	"eqlev"=>true,
	"nobound"=>2,
	"tick1"=>7,"tick2"=>1
      }

    elsif gphys.name == "agcm_comp_tuom"
      levels = NArray[DCL.glpget('rmiss'), -1.5, -1, -0.5, 0, 0.5, 1,DCL.glpget('rmiss')]
      patterns = NArray[30999, 35999, 40999, 55999, 70999, 75999, 85999]
      cont_lev, line_type, label  = cont_lev_set(levels)
      $tone_hash = { 'patterns'=> patterns, 'levels'=>levels }
      $cont_hash = {'levels' => cont_lev, 'index'=> [2,1], 'label'=> label, 
	'line_type'=> line_type }
      $cbar_conf = {
	"levels"=>levels, 
	"colors"=>patterns,
	"eqlev"=>true,
	"nobound"=>0,
	"tick1"=>20,"tick2"=>1
      }

      # 風 (u,-sigdot) ベクトル
      DCL::uglset('LNRMAL', false)
      DCL::ugrset('XFACT1', 2*10E-6*100)
      DCL::ugrset('YFACT1', 15.0/100*100000)
      DCL::uglset('LUNIT', true)
      DCL::ugrset('VXUNIT', 0.06)
      DCL::ugrset('VYUNIT', 0.06)
      DCL::ugrset('VXUOFF', 0.01)

    else
      plot_def(gphys)
    end

  end

end