#!/usr/bin/env ruby

# ----------------------------------------------
# 必要なライブラリ, モジュールの読み込み
require "numru/ggraph"
require "numru/dcl"
include NumRu
include NMath

$fig_path = "./"
$local_path = "/work11/ape/yukiko/lib"
load "#{$local_path}/lib-ape-base.rb"
load "#{$local_path}/lib-ape-view.rb"

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

END{

  $ape = Ape.new("")

 $eqwave = "kelvin"
# $eqwave = "mixed-rossby" 
# $eqwave = "grav-east0"
# $eqwave = "rossby1"
# $eqwave = "grav-east1" 
# $eqwave = "grav-west1" 
# $eqwave = "kelvin_wall1"
# $eqwave = "kelvin_wall2"

  eqwave_calc
  $ape.gropn(1)
  eqwave_plot(1)
  eqwave_var_plot(1)
  $ape.grcls

}


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

# $x_size = 21
$x_size = 161

$eqwave_array = [
  "kelvin",  "mixed-rossby", "grav-east0", 
  "rossby1", "grav-east1", "grav-west1" ,
  "kelvin_wall1", "kelvin_wall2"  
  ]

$title_hash = { 
  "kelvin" => "kelvin wave",
  "mixed-rossby" => "mixed rossby-gravity wave",
  "grav-east0" => "eastward inertio-gravity wave (n=0)",
  "rossby1" => "rossby wave (n=1)",
  "grav-east1" => "eastward inertio-gravity wave (n=1)",
  "grav-west1" => "westward inertio-gravity wave (n=1)",
  "kelvin_wall1" => "kelvin wave",
  "kelvin_wall2" => "kelvin wave"
}

=begin

$eig_hash = {
  "kelvin" => 158,
  "mixed-rossby" => 160,
  "grav-east0" => 156,
  "rossby1" => 161,
  "grav-east1" => 154,
  "grav-west1" => 155,
  "kelvin_wall1" => 157, 
  "kelvin_wall2" => 159
}

=end

$eig_hash = {
  "kelvin" => 164,  
#  "kelvin" => 158,  
#  "kelvin" => 164,  
  "mixed-rossby" => 79, 
  "grav-east0" => 78,
  "rossby1" => 83,
  "grav-east1" => 76,
  "grav-west1" => 77, 
  "kelvin_wall1" => 81, 
  "kelvin_wall2" => 82  
}


# lm=81 2 が外側
# DCL.uumrk(k[4..5],eig[4..5,323]) # 赤道ケルビン1
# DCL.uumrk(k[4..5],eig[4..5,319])  # 赤道ケルビン2
# DCL.uumrk(k[4..5],eig[4..5,325]) # 南北壁のケルビン1
# DCL.uumrk(k[4..5],eig[4..5,317])  # 南北壁のケルビン1
# DCL.uumrk(k[4..5],eig[4..5,321]) # 南北壁のケルビン2
# DCL.uumrk(k[4..5],eig[4..5,331]) # 南北壁のケルビン2
# DCL.uumrk(k[4..5],eig[4..5,322]) # 混合1
# DCL.uumrk(k[4..5],eig[4..5,316])  # 混合2
# DCL.uumrk(k[4..5],eig[4..5,320]) # n=0 東進重力1
# DCL.uumrk(k[4..5],eig[4..5,311])  # n=0 東進重力2
# DCL.uumrk(k[4..5],eig[4..5,314])  # n=1 西進重力1 ??
# DCL.uumrk(k[4..5],eig[4..5,318])  # n=1 西進重力2
# DCL.uumrk(k[4..5],eig[4..5,315])  # n=1 東進重力1
# DCL.uumrk(k[4..5],eig[4..5,313])  # n=2 東進重力1
# DCL.uumrk(k[4..5],eig[4..5,312])  # n=3 西進重力1
# DCL.uumrk(k[4..5],eig[4..5,310])  # n=3 東進重力1
#DCL.uumrk(k[4..5],eig[4..5,327]) # ロスビー
# DCL.uumrk(k[4..5],eig[4..5,326]) # ロスビー
# DCL.uumrk(k[4..5],eig[4..5,324]) # ロスビー1
# DCL.uumrk(k[4..5],eig[4..5,329]) # ロスビー2

$knum = 0.5

# ----------------------------------------------
  
def eqwave_calc

  # ファイル名, 変数の読み込み
  
  filename = "wavecisk.nc"
  p = GPhys::NetCDF_IO.open(filename,"H")
  u = GPhys::NetCDF_IO.open(filename,"U")
  v = GPhys::NetCDF_IO.open(filename,"V")
  eig = GPhys::NetCDF_IO.open(filename,"eigen_val")

  $p = NArray.float($x_size,p.shape[0],p.shape[1]).fill!(1.0)
  $u = NArray.float($x_size,p.shape[0],p.shape[1]).fill!(1.0)
  $v = NArray.float($x_size,p.shape[0],p.shape[1]).fill!(1.0)
  x = NArray.float($x_size).indgen!(-80.0) /5.0/8.0 * PI

  $x_size.times{ |num|
    $p[num,true,true] =  p[true,true,$eig_hash[$eqwave]].data.val * sin(x[num]/2.0)
    $u[num,true,true] =  u[true,true,$eig_hash[$eqwave]].data.val * sin(x[num]/2.0)
    $v[num,true,true] =  v[true,true,$eig_hash[$eqwave]].data.val * cos(x[num]/2.0)
  }

  p.shape[1].times{ |znum|
    $p[true,true,znum] = $p[true,true,znum] / (p[true,znum,$eig_hash[$eqwave]].max - p[true,znum,$eig_hash[$eqwave]].min )
    $u[true,true,znum] = $u[true,true,znum] / (u[true,znum,$eig_hash[$eqwave]].max - u[true,znum,$eig_hash[$eqwave]].min )
    $v[true,true,znum] = $v[true,true,znum] / (v[true,znum,$eig_hash[$eqwave]].max - v[true,znum,$eig_hash[$eqwave]].min )
  }



  x = VArray.new(x).rename("x").put_att("units","1")
  x = Axis.new().set_pos(x)
  y = p.grid_copy.axis(0)
  z = p.grid_copy.axis(1)
  
  grid = Grid.new(x,y,z)
  
  $p = VArray.new($p).rename($eqwave). 
    put_att("units","1").put_att("ape_name",$title_hash[$eqwave])
  $u = VArray.new($u)
  $v = VArray.new($v)
  $p = GPhys.new(grid,$p)
  $u = GPhys.new(grid,$u)
  $v = GPhys.new(grid,$v)


  $u.shape[1].times{ |num|
    if (num % 2) == 0
      $u[true,num,true]  = $u[true,num,true]
      $v[true,num,true]  = $v[true,num,true]
    elsif   
      $u[true,num,true]  = 0.0
      $v[true,num,true]  = 0.0
    end
  }  

  $u.shape[0].times{ |num|
    if (num % 8) == 0
      $u[num,true,true]  = $u[num,true,true]
      $v[num,true,true]  = $v[num,true,true]
    elsif   
      $u[num,true,true]  = 0.0
      $v[num,true,true]  = 0.0
    end
  }  


  $omg = eig[4,$eig_hash[$eqwave]].data.val[0]

end

def eqwave_plot(mkfignum)

  # Hosaka ps の風ベクトル
  DCL::ugpset('LNRMAL', false)
  DCL::ugrset('VXULOC', 0.85)
  DCL::ugrset('VYULOC', 0.23)
  DCL::ugpstx('VXUNIT', 0.05)
  DCL::ugpstx('VYUNIT', 0.05)
  DCL::uglset('LUNIT', true)
  DCL::ugpstx('XFACT1', 0.02/$u.max)
  DCL::ugpstx('YFACT1', 0.02/$u.max)

  $file_label = "k = #{$knum}, omg = #{format("%#.3g", $omg)}"
  $fig_set_hash = { "window" => nil }

  plot_def($p)
  $ape.plot_main($p,$u.data.val,$v.data.val)
#  $ape.plot_main($p.cut(true,true,2),$u.cut(true,true,2).data.val,$v.cut(true,true,2).data.val)
#  $ape.grcls
  
  $file_name = "eqwave_#{$eqwave}"
  if mkfignum == 3 then 
    #  `xwdtopnm dcl*xwd | pnmcut 2 2 904 654  > tmp.pnm`
    `xwdtopnm dcl*xwd  > tmp.pnm`
    `ppmtogif tmp.pnm > #{$fig_path}#{$file_name}.gif`
    `rm tmp.pnm  *xwd `
  elsif mkfignum == 2 then 
    `dclpsrmcm dcl*.ps | dclpsrot > tmp.ps`
    `eps2eps tmp.ps #{$fig_path}#{$file_name}.eps`
    `mv dcl*.ps tmp.ps`
  end
  
end

def eqwave_var_plot(mkfignum)

  $file_label = "k = #{$knum}, omg = #{format("%#.3g", $omg)}"
  $fig_set_hash = { "window" => nil }

  plot_def($p)
  $ape.plot_main($p.cut(true,0,true))
#  $ape.grcls
  
  $file_name = "eqwave_#{$eqwave}"
  if mkfignum == 3 then 
    #  `xwdtopnm dcl*xwd | pnmcut 2 2 904 654  > tmp.pnm`
    `xwdtopnm dcl*xwd  > tmp.pnm`
    `ppmtogif tmp.pnm > #{$fig_path}#{$file_name}.gif`
    `rm tmp.pnm  *xwd `
  elsif mkfignum == 2 then 
    `dclpsrmcm dcl*.ps | dclpsrot > tmp.ps`
    `eps2eps tmp.ps #{$fig_path}#{$file_name}.eps`
    `mv dcl*.ps tmp.ps`
  end
  
end

# トーン, コンターパターンのデフォルト設定
def plot_def(gphys)
  
  # トーンパターンのデフォルト
  patterns = NArray[30999, 40999, 55999, 70999, 75999, 85999]
  
  # トーンレベルのデフォルト: ((最大値 - 最小値) / 7 )
  levels = Array.new
  dx = (gphys.max.val - gphys.min.val)/6 ;  x  = gphys.min.val 
  7.times{ |num| ; levels.push(x) ;  x = dx + x }
  levels[6] = levels[6] + dx ; levels[0] = levels[0] - dx
  levels = NArray.to_na(levels)
  6.times{ |num| 
    if levels[num] < 0 && levels[num+1] > 0 
      if levels[num].abs < levels[num+1].abs
	levels = levels - levels[num]
      else
	levels = levels - levels[num+1]
      end
    end
  }
  
  # カラーバーセットのデフォルト: colorbar.rb 参照
  $cbar_conf = {
    "levels"=>levels, 
    "colors"=>patterns,
    "eqlev"=>true, 
    "nobound"=>0, 
    "tick1"=>20,"tick2"=>1
  }
  
  # コンター, トーンセット
  $tone_hash = { 'patterns'=> patterns, 'levels'=>levels }
  $cont_hash = {'int' => dx/2 }
  
  # ラインセット
  $line_hash = { 'exchange'=>false }
  
  # window セット
  $fig_set_hash = { "window" => nil }
  
end

def gomi
# ----------------------------------------------
# 窓開く

iws = (ARGV[0] || (puts ' WORKSTATION ID (I)  ? ;'; DCL::sgpwsn; gets)).to_i
DCL::gropn(iws)

GGraph.set_fig('viewport'=>[0.13,0.83,0.23,0.58]) # 縦横比は 2 : 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 軸単位を書く位置

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

hoge = NArray.float(21,80).fill!(1.0)
lon = NArray.float(21).indgen!(0.0)/20.0
21.times{ |num|
  hoge[num,true] =  gphys[true,161].data.val * sin(2*PI*lon[num])
}
lon = VArray.new(lon).rename("lon")
hoge = VArray.new(hoge).rename("H")
lon = Axis.new().set_pos(lon)
lat = gphys.grid_copy.axis(0)
grid=Grid.new(lon,lat)
gphys = GPhys.new(grid,hoge)

GGraph.tone( gphys)

#GGraph.line( gphys2[true,-154],false)
#GGraph.line( gphys2[true,154],false,"index"=>21)

# lm=80
# DCL.uumrk(k[9..10],eig[9..10,154])  # n=1 東進重力
# DCL.uumrk(k[9..10],eig[9..10,155])  # n=1 西進重力
# DCL.uumrk(k[8..10],eig[8..10,156])  # n=0 東進重力
# DCL.uumrk(k[8..9],eig[8..9,157])    # 南北壁のケルビン
# DCL.uumrk(k[9..10],eig[9..10,158])    # 赤道ケルビン
# DCL.uumrk(k[8..9],eig[8..9,159])    # 南北壁のケルビン
# DCL.uumrk(k[9..10],eig[9..10,160])  # 混合
# DCL.uumrk(k[9..10],eig[9..10,161])  # n=1 ロスビー

# lm=40
# DCL.uumrk(k[9..10],eig[9..10,34])  # n=1 東進重力
# DCL.uumrk(k[8..9],eig[8..9,36])  # n=1 西進重力
# DCL.uumrk(k[9..10],eig[9..10,37])  # n=0 東進重力
# DCL.uumrk(k[9..10],eig[9..10,38])  # 赤道ケルビン
# DCL.uumrk(k[9..10],eig[9..10,39]) # 混合
# DCL.uumrk(k[9..10],eig[9..10,40]) # 南北壁のケルビン
# DCL.uumrk(k[9..10],eig[9..10,41])  # n=1 ロスビー

end