#!/usr/bin/env ruby

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

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

require "getopts"
require "numru/netcdf"
module NumRu
  class NetCDFVar
    alias put scaled_put
    alias get scaled_get
  end
end

require "numru/ggraph"
require "colorbar"
require 'rational'
require 'complex'
require 'mathn'
include NumRu
include NMath


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

END{

bunsan_calc
bunsan_plot2


}

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

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.12 - 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

end


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

def bunsan_plot

  DCL::gropn(-1)
  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 軸単位を書く位置

  $fig_set_hash = { "window" => [ $x_size[0].real,$x_size[16].real,0.0,4.0 ] }
  GGraph.set_fig($fig_set_hash)
  
  line_index_ary = Array.new; name_ary = Array.new
  
  name_ary = ["kelvin", "mixed", "n=1", "n=2","n=3" ]
  line_index_ary = [12,22,92,32,42]

  # 描画
  GGraph.line( $kelvin.rename("freq"), true, "index" => line_index_ary[0] )
  GGraph.line( $mixed, false, "index" => line_index_ary[1] )
  3.times{ |num|
    GGraph.line( $grav_array[num], false, "index" => line_index_ary[num+2] )
    GGraph.line( $rossby_array[num], false, "index" => line_index_ary[num+2] )
  }
  GGraph.line( $axis0, false, "index" => line_index_ary[0] )

  # タイトル
  tani = "(#{$kelvin.data.get_att("units")})"
  DCL::uxsttl("t", tani , -1.0 )
  DCL::uxmttl("t", "equatorial wave", -1.0 ) 

  $file_label = ""
  
  charsize = 0.81 - $file_label.size * 0.0095 
  DCL::sgtxzv(charsize,0.6,$file_label,0.8 * DCL.uzpget('rsizec1'),0,-1,1) 

  line_index(name_ary,line_index_ary)
  
  DCL::grcls
  
end


def bunsan_plot2


  DCL.swpset('IWIDTH' ,640)
  DCL.swpset('IHEIGHT' ,480)
  DCL::gropn(-1)
  GGraph.set_fig('viewport'=>[0.15,0.55,0.23,0.88]) # 縦横比は 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 軸単位を書く位置

  $fig_set_hash = { "window" => [ $x_size[8].real,$x_size[16].real,-4.0,4.0 ] }
  GGraph.set_fig($fig_set_hash)
  
  line_index_ary = Array.new; name_ary = Array.new
  
  name_ary = ["kelvin", "mixed", "n=1", "n=2","n=3" ]
  line_index_ary = [12,22,92,32,42]

  # 描画
  GGraph.line( $kelvin.rename("freq")[8..16], true, 
	      "index" => line_index_ary[0] )
  GGraph.line( $mixed[8..16], false, "index" => line_index_ary[1] )
  GGraph.line( $mixed2[8..16], false, "index" => line_index_ary[1] )
  3.times{ |num|
    GGraph.line( ($grav_array[num])[8..16], false, 
		"index" => line_index_ary[num+2] )
    GGraph.line( ($grav2_array[num])[8..16], false, 
		"index" => line_index_ary[num+2] )
    GGraph.line( ($rossby2_array[num])[8..16], false, 
		"index" => line_index_ary[num+2] )
  }
  GGraph.line( $axis0[8..16], false, "index" => line_index_ary[0] )

  # タイトル
  tani = "(#{$kelvin.data.get_att("units")})"
  DCL::uxsttl("t", tani , -1.0 )
  DCL::uxmttl("t", "equatorial wave", -1.0 ) 

  $file_label = ""
  
  charsize = 0.81 - $file_label.size * 0.0095 
  DCL::sgtxzv(charsize,0.6,$file_label,0.8 * DCL.uzpget('rsizec1'),0,-1,1) 

  line_index(name_ary,line_index_ary)
  
  DCL::grcls
  
end

def line_index(str_ary,line_index_ary)
  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.15
  vy = 0.12 - charsize/2
  str_ary.size.times{|num|
    DCL::sgtxzv(vx,vy,"--- #{str_ary[num]}", 
		charsize, 0, -1,line_index_ary[num])
    vy -= dvy
    if num == 3
      vx = 0.30
      vy = 0.12 - charsize/2
    end
  }
  nil
end


def bunsan_calc

  $x_size = (NArray.sfloat(17).indgen - 8.0) *0.5
  # omega = 2*3.1415/24/60/60
  # beta =  2*(omega)/(4e+7) * ( 60*60*24*(6370e+3) )  # (2*Omega/a)*cos(0)
  # c_g = ((200*9.8)**0.5)*60*60*24/(4e+7)

  beta = 1.0
  c_g = 1.0
  
  # 軸
  knum = VArray.new($x_size).rename("knum").put_att("units","1")
  knum = Axis.new().set_pos(knum)
  grid = Grid.new(knum)

  # 軸
  $axis0 = $x_size * 0.0
  $axis0 = VArray.new($axis0).rename("axis0").
    put_att("units","1").put_att("ape_name","axis0")
  $axis0 = GPhys.new( grid, $axis0)

  ## 虚数を含む軸に変換
  x_size = (NArray.complex(17).indgen - 8.0) *0.5

  # 初期化
  $rossby_array = Array.new
  $rossby2_array = Array.new
  $grav_array = Array.new
  $grav2_array = Array.new

  3.times{ |num|

    num = num +1
    
    ## 3 次方程式解法: カルダノの公式
    ## x^3 + 3px + q
    ## t^2 + qt - p =0, t = -q/2 +- (q**2/4 + p)**0.5
    p = - ( (c_g*x_size)**2 + c_g*beta*(2*num+1)) / 3
    q = - c_g**2*beta*x_size
    t1  =  -q/2 + (q**2/4 + p**3)**0.5
    t2  =  -q/2 - (q**2/4 + p**3)**0.5
    t1 = t1**(1.0/3)
    t2 = t2**(1.0/3)

    # 重力波
    grav1 = (t1 + t2).real
    grav1 = VArray.new(grav1).rename("grav1").
      put_att("units","1").put_att("ape_name","gravity")
    grav1 = GPhys.new( grid, grav1)
    
    # 反対重力波 (図示しない)
    grav2 = (t1 *  Complex.polar(1, Math::PI*2/3) + 
	     t2 * Complex.polar(1, Math::PI*4/3)).real
    grav2 = VArray.new(grav2).rename("grav2").
      put_att("units","1").put_att("ape_name","gravity")
    grav2 = GPhys.new( grid, grav2)
    
    # ロスビー波
    rossby = (t2 * Complex.polar(1, Math::PI*2/3) + 
	      t1 * Complex.polar(1, Math::PI*4/3) ).real
    rossby2 = rossby
    rossby = (rossby > 0) * rossby
    rossby = VArray.new(rossby).rename("rossby").
      put_att("units","1").put_att("ape_name","rossby")
    rossby = GPhys.new( grid, rossby)

    rossby2 = VArray.new(rossby2).rename("rossby").
      put_att("units","1").put_att("ape_name","rossby")
    rossby2 = GPhys.new( grid, rossby2)

    $rossby_array.push(rossby)
    $rossby2_array.push(rossby2)
    $grav_array.push(grav1)
    $grav2_array.push(grav2)
  
  }

  # 混合ロスビー
  $mixed = c_g*x_size/2 + ((c_g*x_size/2)**2 + c_g*beta )**0.5
  $mixed = VArray.new($mixed).rename("mixed").
    put_att("units","1").put_att("ape_name","mixed")
  $mixed = GPhys.new( grid, $mixed)

  # 反対混合ロスビー (図示しない)
  $mixed2 = c_g*x_size/2 - ((c_g*x_size/2)**2 + c_g*beta )**0.5
  $mixed2 = VArray.new($mixed2).rename("mixed").
    put_att("units","1").put_att("ape_name","mixed")
  $mixed2 = GPhys.new( grid, $mixed2)

  # ケルビン波
  $kelvin = (c_g*x_size).real
  $kelvin = ($kelvin > 0) * $kelvin
  $kelvin = VArray.new($kelvin).rename("kelvin").
    put_att("units","1").put_att("ape_name","kelvin")
  $kelvin = GPhys.new( grid, $kelvin)
  
end


=begin
class Ape

  def plot_spct_bunsan(gphys)

    gropn(1) if $gropn == nil 

    if gphys.data.get_att("long_name")
      gphys = gphys.set_att("long_name","")
    end


    GGraph.line( gphys, true , $tone_hash ) 


    # 分散曲線
    x_size = gphys.grid_copy.axis(0).pos.val.size
    x = gphys.grid_copy.axis(0).pos.val[x_size/2..-1]
    h_array = NArray.to_na([ 12, 25, 50, 100, 200 ])

    omega = 2*3.1415/24/60/60
    beta =  2*(omega)/(4e+7) * ( 60*60*24*(6370e+3) )  # (2*Omega/a)*cos(0)


    for i in (0..4)
      
      h = h_array[i]
      c_g = ((h*9.8)**0.5)*60*60*24/(4e+7)    # 重力波速度
      
      kelvin = x*c_g
      rossby = -x/(x**2/beta +(2+1)/c_g )
      grav = -(c_g*( c_g*x**2 + beta*(2+1) ) )**0.5

      DCL::uuslnt(1)
      ii=(i)*20+4
      ii = 5*20+4 if ii == 2*20+4
      ii = 720+4 if ii == 4*20+4
      DCL::uuslni(ii)
      DCL::uulin(x, kelvin)
      DCL::uulin(-x, -rossby)
      DCL::uulin(-x, -grav)
      DCL::uulin(x, -grav)

    end

    name_ary = ["h = 12 ", "h = 25 ","h = 50 ","h =100", "h =200 "]
    line_index_ary = [4,24,104,64,724]

    __line_index(name_ary,line_index_ary) 

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

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

    # トーンバー
    unless gphys.axnames.size == 1
      DCL::Util::color_bar($cbar_conf)  unless $tone_flag == false
    end

  end  
  

end
=end