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

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

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

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

END{

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

#  ["con","mradlAa","mradlAb","mradlAc","mradlAd"].each{ |item|
  ["con","a","b","c","d"].each{ |item|
    $expid = item
    $cum = "adj"
    set_path
    a.nc_tr_agcm5
    a.nc_tr_wkplot_agcm5

    $expid = item
    $cum = "kuo"
    set_path
    a.nc_tr_agcm5
    a.nc_tr_wkplot_agcm5
  }


}

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


def set_path

  $fig_path = "/home/yukiko/work/aqua3d/yukiko/sh/ruby/mradlstdexp/mlshtr/figs/"
  $ncfile_path = "/home/yukiko/work/aqua3d/yukiko/mradl/data/wk-plot-data/"
#  $file_label = "AGCM5_#{$cum}_#{$expid}"

  if $expid == "con"
    $file_label = "#{$cum}_control"
#    $file_label = "#{$cum}_con"
  else
    $file_label = "#{$cum}_#{$expid}"
    $expid = "mradlA" + $expid
  end

  $groupid = "AGCM5-#{$cum}-#{$expid}"
end

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

class Ape_mkfig
  # #{$groupid}_TR_control.nc
  def nc_tr_agcm5

    lat_lost_axis = ["latitude=1.395307"]
      
    (1..16).each { |varname|
        
      @data = 
        Ape.new("#{$ncfile_path}/#{$cum}700#{$expid}400/u_out_z#{varname}.nc")
        
      name = "tr_u_z#{varname}"
        
      if name =~ /tr_u/
      
        if name == "tr_u_z1"
          lost_axis = ["sigma=0.994997"]
        elsif name == "tr_u_z2"
          lost_axis = ["sigma=0.9799879"]
        elsif name == "tr_u_z3"
          lost_axis = ["sigma=0.9499499"]
        elsif name == "tr_u_z4"
          lost_axis = ["sigma=0.8998809"]
        elsif name == "tr_u_z5"
          lost_axis = ["sigma=0.8297704"]
        elsif name == "tr_u_z6"
          lost_axis = ["sigma=0.7446762"]
        elsif name == "tr_u_z7"
          lost_axis = ["sigma=0.6495416"]
        elsif name == "tr_u_z8"
          lost_axis = ["sigma=0.549458"]
        elsif name == "tr_u_z9"
          lost_axis = ["sigma=0.4544692"]
        elsif name == "tr_u_z10"
          lost_axis = ["sigma=0.3694841"]
        elsif name == "tr_u_z11"
          lost_axis = ["sigma=0.2945043"]
        elsif name == "tr_u_z12"
          lost_axis = ["sigma=0.2295326"]
        elsif name == "tr_u_z13"
          lost_axis = ["sigma=0.1745732"]
        elsif name == "tr_u_z14"
          lost_axis = ["sigma=0.1244002"]
        elsif name == "tr_u_z15"
          lost_axis = ["sigma=0.07398598"]
        elsif name == "tr_u_z16"
          lost_axis = ["sigma=0.02074752"]
        end
        
      # x-t ダイヤグラム
      gphys = @data.gphys_open("u").rename(name)
      gphys =  gphys.cut(true,0,true)[true,-400..-1]
      
      gphys = gphys.
        set_att("ape_name","eastward_wind").
        set_att("units","m s-1")
        
	grid_1 =
	  Axis.new().
	  set_pos(VArray.new(NArray.sfloat(400).indgen!/4).rename("time").
            put_att("units","days"))
	grid_0 = gphys.grid_copy.axis(0)
        grid = Grid.new(grid_0,grid_1)
        gphys = GPhys.new(grid, gphys.data )

	mkfig_plot((gphys - gphys.mean(0)).
		   add_lost_axes(lat_lost_axis).
		   add_lost_axes(lost_axis).
		   add_lost_axes("(diff) from (mean) zonal"))

	# 時空間スペクトル
	gphys = gphys.stspct_fft("#{name}_spct").
	  set_att("ape_name",
		  "#{gphys.data.get_att("ape_name")}_S-T_power_spectrum").
	  set_att("units","(#{gphys.data.get_att("units")} day-1)2").
          add_lost_axes(lat_lost_axis).
	  add_lost_axes(lost_axis)
	
	dim1 = gphys.coord(0).shape.to_s.to_i
	dim2 = gphys.coord(1).shape.to_s.to_i
	mkfig_plot(gphys[((dim1+1)/2-31)..((dim1+1)/2+29),0..80] )
	
      end
    }
  end

  def nc_tr_wkplot_agcm5

    @data = Ape.new("#{$ncfile_path}/data/#{$groupid}-u_wk_40smooth.nc")
    @data.netcdf_open.var_names.each { |name|
      unless name == "wvn" || name == "freq" || name == "freq0" || name == "noname"
      
        if name == "tr_u_z1"
          lost_axis = ["sigma=0.994997"]
        elsif name == "tr_u_z2"
          lost_axis = ["sigma=0.9799879"]
        elsif name == "tr_u_z3"
          lost_axis = ["sigma=0.9499499"]
        elsif name == "tr_u_z4"
          lost_axis = ["sigma=0.8998809"]
        elsif name == "tr_u_z5"
          lost_axis = ["sigma=0.8297704"]
        elsif name == "tr_u_z6"
          lost_axis = ["sigma=0.7446762"]
        elsif name == "tr_u_z7"
          lost_axis = ["sigma=0.6495416"]
        elsif name == "tr_u_z8"
          lost_axis = ["sigma=0.549458"]
        elsif name == "tr_u_z9"
          lost_axis = ["sigma=0.4544692"]
        elsif name == "tr_u_z10"
          lost_axis = ["sigma=0.3694841"]
        elsif name == "tr_u_z11"
          lost_axis = ["sigma=0.2945043"]
        elsif name == "tr_u_z12"
          lost_axis = ["sigma=0.2295326"]
        elsif name == "tr_u_z13"
          lost_axis = ["sigma=0.1745732"]
        elsif name == "tr_u_z14"
          lost_axis = ["sigma=0.1244002"]
        elsif name == "tr_u_z15"
          lost_axis = ["sigma=0.07398598"]
        elsif name == "tr_u_z16"
          lost_axis = ["sigma=0.02074752"]
        end
	
	# 時空間スペクトル
	gphys = @data.gphys_open(name). 
	  add_lost_axes(lost_axis)
	dim1 = gphys.coord(0).shape_current.to_s.to_i
	dim2 = gphys.coord(1).shape_current.to_s.to_i

	if name =~ /bg_spct/
	  mkfig_plot(gphys[((dim1+1)/2-31)..((dim1+1)/2+29),0..80].log.
                       add_lost_axes("logarithm of data") )
	else
	  mkfig_plot(gphys[((dim1+1)/2-31)..((dim1+1)/2+29),0..80])
	end
        
      end
    }
    
  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","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')]
# 青色
#      levels = NArray[0.0000, 0.0001, 0.0002, 0.0006, 1000]
#      patterns = NArray[99001, 40999, 30999, 20999]
      $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 == "tr_u_z1" || gphys.name == "tr_u_z2" || gphys.name == "tr_u_z3" || gphys.name == "tr_u_z4" || gphys.name == "tr_u_z5" || gphys.name == "tr_u_z6" || gphys.name == "tr_u_z7" || gphys.name == "tr_u_z8" || gphys.name == "tr_u_z9" || gphys.name == "tr_u_z10" || gphys.name == "tr_u_z11" || gphys.name == "tr_u_z12" || gphys.name == "tr_u_z13" || gphys.name == "tr_u_z14" || gphys.name == "tr_u_z15" || gphys.name == "tr_u_z16" 
      levels = NArray[DCL.glpget('rmiss'), -10, -5, 0, 5, 10, 20, DCL.glpget('rmiss')]
      patterns = NArray[30999, 35999, 40999, 55999, 70999, 75999, 85999]
      $tone_hash = { 'patterns'=> patterns, 'levels'=>levels }
      $cbar_conf = {
	"levels"=>levels, 
	"colors"=>patterns,
	"eqlev"=>true,
	"nobound"=>0,
	"tick1"=>20,"tick2"=>1
      }

    elsif gphys.name == "tr_u_z1_spct" || gphys.name == "tr_u_z2_spct" || gphys.name == "tr_u_z3_spct" || gphys.name == "tr_u_z4_spct" || gphys.name == "tr_u_z5_spct" || gphys.name == "tr_u_z6_spct" || gphys.name == "tr_u_z7_spct" || gphys.name == "tr_u_z8_spct" || gphys.name == "tr_u_z9_spct" || gphys.name == "tr_u_z10_spct" || gphys.name == "tr_u_z11_spct" || gphys.name == "tr_u_z12_spct" || gphys.name == "tr_u_z13_spct" || gphys.name == "tr_u_z14_spct" || gphys.name == "tr_u_z15_spct" || gphys.name == "tr_u_z16_spct" 
      levels = NArray[1e-3, 1e-2, 1e-1, 1, 10 ]
      patterns = NArray[40999, 30999, 20999,17999]
      $tone_hash = { 'patterns'=> patterns, 'levels'=>levels }
      $cont_hash = {'int'=>1e-1}
      $cbar_conf = {
	"levels"=>levels, 
	"colors"=>patterns,
	"eqlev"=>true,
	"nobound"=>2,
	"tick1"=>7,"tick2"=>1,
      "vx0"=>0.05,
      "vy0"=>0.1,
      "label_size"=>0.015
      }

      # W-K plot
    elsif gphys.name =~ /_sym_/
      levels = NArray[1.2, 1.4, 1.6, 2.0, DCL.glpget('rmiss')]
      patterns = NArray[40999, 30999, 20999,17999]
      $tone_hash = { 'patterns'=> patterns, 'levels'=>levels }
      $cbar_conf = {
	"levels"=>levels, 
	"colors"=>patterns,
	"eqlev"=>true,
	"nobound"=>2,
	"tick1"=>7,"tick2"=>1,
      "vx0"=>0.05,
      "vy0"=>0.1,
      "label_size"=>0.015
      }

    elsif gphys.name =~ /_asym_/ 
      levels = NArray[0.1, 0.3, 0.5, 0.7, DCL.glpget('rmiss')]
      patterns = NArray[40999, 30999, 20999,17999]
      $tone_hash = { 'patterns'=> patterns, 'levels'=>levels }
      $cbar_conf = {
	"levels"=>levels, 
	"colors"=>patterns,
	"eqlev"=>true,
	"nobound"=>2,
	"tick1"=>7,"tick2"=>1,
      "vx0"=>0.05,
      "vy0"=>0.1,
      "label_size"=>0.015
      }

    else
      plot_def(gphys)
    end

  end

end