#!/usr/bin/env ruby

# ----------------------------------------------
# local load path

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

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

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

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


END{

#  a = Ape_mkfig.new 3
  a = Ape_mkfig.new 2

#  [ "T39L48_eml", "T39L48_non", "T159L48_eml", "T159L48_non"].each { |rezol| 
  ["T39L48_non"].each { |rezol| 
    $rezol = rezol
    set_path
    a.gr_comp_tuv
#    a.gr_comp_tuv("tconv")
#    a.gr_comp_tuv("novect")

  }    

}

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


def set_path
  $expID = "control"
  $fig_path = "/home/yukiko/work/ape/yukiko/figs/tmp/"
  $gr2ncfile_path = "/home/yukiko/work/ape/yukiko/data/spctfilt/"
#  $compfile_path = "./netcdf/"
  $compfile_path = "/home/yukiko/work/ape/yukiko/data/spctfilt/"
  $file_label = "#{$rezol}_control"
  $groupid = "AGUforAPE-03a"
end

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


class Ape_mkfig
  
  def gr_comp_tuv(compmethod=nil)
    
    if compmethod == nil 
      compmethod="threshold"
    end
    
    filt_array = ["ad","k","g"]
    filt_hash =  { 
      "ad" => "advect", 
      "k"  => "kelvin",
      "g"  => "grav"
    }

    filt_title_hash =  { 
      "ad" => "advection", 
      "k"  => "kelvin wave",
      "g"  => "gravity wave"
    }


    filt_array.each{ |filter|

      $filt = filter

      file_name = 
        "#{$compfile_path}/#{$rezol}_control_#{filt_hash[filter]}filt-composite.nc"
    
      @data = Ape.new(file_name)
      t_comp = @data.gphys_open("comp_t_#{filter}filt").
        rename("comp_tuom_#{filter}filt_mono")
      tconv_comp = @data.gphys_open("comp_tconv_#{filter}filt").
        rename("comp_tconvuom_#{filter}filt_mono") 
      u_comp = @data.gphys_open("comp_u_#{filter}filt").
        rename("comp_uuom_#{filter}filt_mono") 
      om_comp = @data.gphys_open("comp_om_#{filter}filt")

      lost_axes_diff = t_comp.get_att("lost_axis").to_s.split("\n")
      lost_axes = tconv_comp.get_att("lost_axis").to_s.split("\n")
      
      t_comp = t_comp.set_lost_axes(lost_axes_diff)
      u_comp = u_comp.set_lost_axes(lost_axes_diff)
      tconv_comp = tconv_comp.set_lost_axes(lost_axes)
      
      # ベクトルの間引
      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)
      
      if compmethod =~ /_large/ 
        num_l = u_comp.axis(0).pos.val.size/60
      else
        num_l = u_comp.axis(0).pos.val.size/120
      end
      num_z = u_comp.axis(1).pos.val.size/48
    
      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
      }

      # 描画
      dim = t_comp.axis(0).pos.val.size

      if compmethod =~ /_large/
        mkfig_plot(t_comp,
	           u_data,
	           -om_data )

        mkfig_plot(u_comp,
	           u_data,
	           -om_data )

        mkfig_plot(tconv_comp, 
	           u_data, 
	           -om_data )
      elsif compmethod == "novect"

        t_comp = t_comp.
       set_att("ape_name","#{filt_title_hash[filter]} filter composite temperature").
       set_att("units","K").
        rename("comp_t_#{filter}filt_mono")
        mkfig_plot(t_comp[(dim/3)..(dim*2/3),true] )

        u_comp = u_comp.
       set_att("ape_name","#{filt_title_hash[filter]} filter composite eastward_wind").
       set_att("units","m s-1").
        rename("comp_u_#{filter}filt_mono")
        mkfig_plot(u_comp[(dim/3)..(dim*2/3),true] )

      elsif compmethod == "tconv"

        tconv_comp = tconv_comp.rename("comp_tconv_cont")

        t_comp = t_comp.
       set_att("ape_name","composite (temperature, condensation_heating)").
       set_att("units","K, K s-1").
        rename("comp_t-tconv_#{filter}filt_mono")
        mkfig_plot_tconv(t_comp[(dim/3)..(dim*2/3),true],
                   tconv_comp[(dim/3)..(dim*2/3),true] )

        u_comp = u_comp.
       set_att("ape_name","composite (eastward_wind, condensation_heating)").
       set_att("units","m s-1, K s-1").
        rename("comp_u-tconv_#{filter}filt_mono")
        mkfig_plot_tconv(u_comp[(dim/3)..(dim*2/3),true],
                   tconv_comp[(dim/3)..(dim*2/3),true] )

      else
        mkfig_plot(t_comp[(dim/3)..(dim*2/3),true],
	           u_data[(dim/3)..(dim*2/3),true], 
	           -om_data[(dim/3)..(dim*2/3),true])

        mkfig_plot(u_comp[(dim/3)..(dim*2/3),true],
	           u_data[(dim/3)..(dim*2/3),true], 
	           -om_data[(dim/3)..(dim*2/3),true])

        mkfig_plot(tconv_comp[(dim/3)..(dim*2/3),true],
	           u_data[(dim/3)..(dim*2/3),true], 
	           -om_data[(dim/3)..(dim*2/3),true])
      end

     }

  end

end


class Ape
  
  def plot_tconv(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_tconv(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


  # ps ファイル出力
  def mkps_tconv(gphys, gphys_u=nil, gphys_v=nil)
    'rm dcl*.ps'
    gropn(2) if $gropn == nil
    plot_tconv(gphys, gphys_u)

    gphys = gphys[0]  if gphys.class == Array
      
    file_name = "#{$file_label}_#{gphys.name}.ps"
    `mv dcl*.ps  #{$fig_path}#{file_name}` 
#    `dclpsrmcm dcl*.ps | dclpsrot > tmp.ps`
#    `eps2eps tmp.ps  #{$fig_path}#{file_name}` 
    `rm tmp.ps dcl*.ps`

  end

  # ps ファイル出力
  def mkps(gphys, gphys_u=nil, gphys_v=nil)
    'rm dcl*.ps'
    gropn(2) 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}_#{gphys.name}.ps"
    `mv dcl*.ps  #{$fig_path}#{file_name}` 

#    `dclpsrmcm dcl*.ps | dclpsrot > tmp.ps`
#    `eps2eps tmp.ps  #{$fig_path}#{file_name}` 
    `rm tmp.ps dcl*.ps`
  end

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

    plot_tconv(gphys, gphys_u)

    gphys = gphys[0]     if gphys.class == Array
    
    file_name = "#{$file_label}_#{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 830 555 > tmp.pnm;done`
    `ppmtogif tmp.pnm > #{$fig_path}#{file_name}`
    `rm dcl_* tmp.pnm`  

    $file_name = $pre_file

  end

  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}_#{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 830 555 > tmp.pnm;done`
    `ppmtogif tmp.pnm > #{$fig_path}#{file_name}`
    `rm dcl_* tmp.pnm`  

    $file_name = $pre_file

  end


  def plot_main(gphys, gphys_u=nil, gphys_v=nil)

    # 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


    # 風 (u,-sigdot) ベクトル
    DCL::ugvect( gphys_u, gphys_v ) unless gphys_u == nil || gphys_v == nil 

    # タイトル
    DCL::uzinit 
    if $filt == "ad"
      DCL::uxmttl("t", "Adv. filter", -1.0 ) 
    elsif $filt == "g"
      DCL::uxmttl("t", "GW filter", -1.0 ) 
    elsif $filt == "k"
      DCL::uxmttl("t", "KW filter", -1.0 ) 
    else
      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
      DCL::uzinit 
      $file_label = "#{@file}@#{gphys.name}"  if $file_label == "filename"
      DCL.uxpttl("t", 0, $file_label, 1.0)
    end

    # トーンバー
    unless gphys.axnames.size == 1
      unless $tone_hash['levels'][0] == $tone_hash['levels'][1] 
	if gphys.name == "comp_tuom_mono"
	else
#	  DCL::Util::color_bar($cbar_conf)  unless $tone_flag == false
###          GGraph.color_bar($colorbar_hash)
	end
      end
    end

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

  end


  def plot_main_tconv(gphys, gphys_u=nil)

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

    GGraph.set_fig($fig_set_hash)

    plot_set(gphys_u)  
    GGraph.tone( gphys_u, true , $tone_hash ) 
    GGraph.color_bar($colorbar_hash)

    plot_set(gphys)  
    GGraph.contour( gphys, false, $cont_hash )



    # タイトル
    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
    DCL::uzinit 
    $file_label = "#{@file}@#{gphys.name}"  if $file_label == "filename"
    DCL.uxpttl("t", 0, $file_label, 1.0)
    
  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 
      NumRu::DCL.udpset("LMSG",false)
    else
      $tone_flag = true
      NumRu::DCL.udpset("LMSG",true)
    end
    NumRu::DCL.udpset("LMSG",false)

    # 初期化
    $fig_set_hash = {"window" => nil}
    $line_hash = { 'exchange'=>false, 'index'=>3 }
    DCL::ugrset('UXUNIT', -999)
    DCL::ugrset('UYUNIT', -999)

    rmiss = DCL.glpget('rmiss')

    if gphys.name =~ /comp_tuom/ || gphys.name =~ /comp_t-tconv/ || gphys.name =~ /comp_t_/
      levels = NArray[rmiss, -1.5, -1, -0.5, 0, 0.5, 1,rmiss]
      patterns = NArray[603, 603, 603, 603, 99001, 99001, 99001]
#      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)
      DCL::uglset('LUNIT', true)
      DCL::uglset('LUNIT', false)
      DCL::ugrset('VXUNIT', 0.06)
      DCL::ugrset('VYUNIT', 0.06)
      DCL::ugrset('VXUOFF', 0.01)
      DCL::uglset('LMSG', false)

    elsif gphys.name =~ /comp_tconvuom/
      levels = NArray[rmiss, -3e-4, 0, 3e-4, 6e-4, 9e-4, 12e-4,rmiss]
      patterns = NArray[603, 603, 99001, 99001, 99001, 99001, 99001]
      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)
      DCL::uglset('LUNIT', true)
      DCL::uglset('LUNIT', false)
      DCL::ugrset('VXUNIT', 0.06)
      DCL::ugrset('VYUNIT', 0.06)
      DCL::ugrset('VXUOFF', 0.01)
      DCL::uglset('LMSG', false)

    elsif gphys.name =~ /comp_tconv_cont/
#      levels = NArray[1e-5, 2e-5, rmiss]
#      patterns = NArray[603,603,603]
      levels = NArray[rmiss, 1e-5, rmiss]
      patterns = NArray[99000, 603]

#      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)
      DCL::uglset('LUNIT', true)
      DCL::uglset('LUNIT', false)
      DCL::ugrset('VXUNIT', 0.06)
      DCL::ugrset('VYUNIT', 0.06)
      DCL::ugrset('VXUOFF', 0.01)
      DCL::uglset('LMSG', false)

    elsif gphys.name =~ /comp_uuom/ || gphys.name =~ /comp_u-tconv/ || gphys.name =~ /comp_u_/
      levels = NArray[rmiss, -9, -6, -3, 0, 3, 6, rmiss]
      patterns = NArray[603, 603, 603, 603, 99001, 99001, 99001]
#      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
      }

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

    else
      plot_def(gphys)
    end
    
  end

end

class GPhys
  def cshift(n) 
    y = self.dup
    y[n..-1,true] = self[0..-1-n,true]
    y[0..n-1,true] = self[-n..-1,true]
    return y
  end
end


class Ape_mkfig

  def mkfig_plot_tconv(gphys, gphys_u=nil, gphys_v=nil)

    if @fig_num == 3 then
        @data.mkdump_tconv(gphys, gphys_u)
    elsif @fig_num == 2 then
        @data.mkps_tconv(gphys, gphys_u)
    end

  end

end