#!/usr/bin/env ruby

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

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

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

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

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

class Ape

  def plot_spct_bunsan(gphys)

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

    GGraph.tone( gphys, true , $tone_hash ) 
#    GGraph.contour( gphys, false, $cont_hash ) unless $cont_flag == false

    # 分散曲線
    $x_size = gphys.grid_copy.axis(0).pos.val
    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)

    name_ary = ["h = 12 ", "h = 25 ","h = 50 ","h =100", "h =200 "]
    if $gropn == 2
      line_index_ary = [4,24,104,64,724]
#      line_index_ary = [1,1,1,1,1]
      line_type_ary =  [1,1,1,1,1]
    else
      line_index_ary = [1,21,101,61,721]
#      line_index_ary = [1,1,1,1,1]
      line_type_ary =  [1,1,1,1,1]
    end

    for i in (0..4)
      
      h = h_array[i]

      if gphys.name =~ /asym/ 
	$sym_num = 2
      else
	$sym_num = 1
      end

      DCL.sgpset('LCLIP', true) 
      bunsan_calc(h)

      GGraph.line( $grav1, false, "index" => line_index_ary[i] )

      if gphys.name =~ /asym/ 
	GGraph.line( $mixed, false, "index" => line_index_ary[i] )
      else
	GGraph.line( $rossby, false, "index" => line_index_ary[i] )
	GGraph.line( $kelvin, false, "index" => line_index_ary[i] )
      end

    end
    GGraph.line( $grav1*0.0, false, "index" => line_index_ary[0] )
#    __line_index(name_ary,line_index_ary) 
    DCL.sgpset('LCLIP', false) 
    __show_line_index(name_ary,line_index_ary, line_type_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 ) 

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

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

  end  


  #  ラインインデックスの種類を表示
  def __show_line_index(str_ary,line_index_ary, line_type_ary, x=0.15, y=0.1)
    charsize = 0.8 * DCL.uzpget('rsizec1')
    len = 0.04     # 線の長さ
    dvx = 0.03
    dvy = charsize*1.6
    
    raise TypeError,"Array expected" if ! str_ary.is_a?(Array)
    vxmin,vxmax,vymin,vymax = DCL.sgqvpt
    vx = 0.15
    vy = 0.13 - charsize/2
    str_ary.size.times{|num|
      DCL::sgplzv([vx, vx + len],[vy, vy], line_type_ary[num],
		  line_index_ary[num])
      DCL::sgtxzv(vx + len + 0.01 ,vy,str_ary[num],
                  charsize, 0, -1, 13 )
      vy -= dvy
      if num == 2
#      if num == 3
        vx = 0.34
        vy = 0.13 - charsize/2
      elsif num == 5
#      elsif num == 7
        vx = 0.53
        vy = 0.13 - charsize/2
      end
    }
    nil
  end


  def bunsan_calc(h=200)

    omega = 2*3.1415/24/60/60
    beta =  2*(omega)/(4e+7) * ( 60*60*24*(6370e+3) )  # (2*Omega/a)*cos(0)
    c_g = ((h*9.8)**0.5)*60*60*24/(4e+7)
  
    # 軸
    knum = VArray.new($x_size).rename("knum").put_att("units","1")
    knum = Axis.new().set_pos(knum)
    grid = Grid.new(knum)
        
    ## 虚数を含む軸に変換
    x_size = NArray.complex($x_size.size).fill!(1.0) * $x_size

    # モード
    num = $sym_num
    
    ## 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)
        
    # ケルビン波
    $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)

    # 混合ロスビー
    $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)
  
  end

end