#!/usr/bin/env ruby

$local_path = "/home/yukiko/tmp/ape-data/lib/"
$fig_path = "/home/yukiko/tmp/ape-data/figs/tmp/"

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

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

END{

  $ape = Ape.new("")

 $eqwave = "kelvin"
# $eqwave = "mixed-grav" 
# $eqwave = "mixed-rossby" 
# $eqwave = "grav-east0"
# $eqwave = "rossby1"
# $eqwave = "grav-east1" 
# $eqwave = "grav-west1" 
# $eqwave = "rossby2" 
# $eqwave = "grav-east2" 
# $eqwave = "grav-west2" 

  eqwave_calc
  eqwave_plot(1)

}

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

  $eqwave_array = [
    "rossby2", "grav-east2", "grav-west2", 
    "kelvin", "mixed-grav", "mixed-rossby", "grav-east0", 
    "rossby1", "grav-east1", "grav-west1" 
  ]

$title_hash = { 
  "kelvin" => "kelvin wave",
  "mixed-grav" => "mixed rossby-gravity 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)",
  "rossby2" => "rossby wave (n=2)",
  "grav-east2" => "eastward inertio-gravity wave (n=2)",
  "grav-west2" => "westward inertio-gravity wave (n=2)"
}

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

def eqwave_calc

# Kelvin
## n = -1,k = 0.5,omega = 0.5

  $x_size = (NArray.sfloat(33).indgen - 16.0) * 3.14159/8.0
  $y_size = (NArray.sfloat(19).indgen - 9.0) / 3.0
  
  # 軸
  x = VArray.new($x_size).rename("x").put_att("units","1")
  y = VArray.new($y_size).rename("y").put_att("units","1")
  x = Axis.new().set_pos(x)
  y = Axis.new().set_pos(y)
  grid = Grid.new(x,y)

  $p = NArray.sfloat(33,19).fill!(0.0)
  $u = NArray.sfloat(33,19).fill!(0.0)
  $v = NArray.sfloat(33,19).fill!(0.0)


  bunsan_calc

  33.times{ |num|

    if $eqwave == "kelvin" then
      
      $p[num,true] = - exp(-0.5*($y_size**2))*sin($x_size[num]*$knum) 
      $u[num,true] = - exp(-0.5*($y_size**2))*sin($x_size[num]*$knum) 

    elsif $eqwave == "mixed-grav" then
    
      $p[num,true] = - 2*$y_size*exp(-0.5*($y_size**2))*
	sin($x_size[num]*$knum) 
      $u[num,true] = - 2*$y_size*exp(-0.5*($y_size**2))*
	sin($x_size[num]*$knum) 
      $v[num,true] = - 2*(-$omg + $knum)*exp(-0.5*($y_size**2))*
	cos($x_size[num]*$knum)
      
    elsif $eqwave == "mixed-rossby" then

      $p[num,true] = -2*$y_size*exp(-0.5*($y_size**2))*
	sin($x_size[num]*$knum) 
      $u[num,true] = -2*$y_size*exp(-0.5*($y_size**2))*
	sin($x_size[num]*$knum) 
      $v[num,true] = - 2*(-$omg + $knum)*exp(-0.5*($y_size**2))*
	cos($x_size[num]*$knum)

    elsif $eqwave == "grav-east0" then

      $p[num,true] = -2*$y_size*exp(-0.5*($y_size**2))*
	sin($x_size[num]*$knum) 
      $u[num,true] = -2*$y_size*exp(-0.5*($y_size**2))*
	sin($x_size[num]*$knum) 
      $v[num,true] = -2*(-$omg + $knum)*exp(-0.5*($y_size**2))*
	cos($x_size[num]*$knum)
     
    elsif $eqwave == "rossby1" then

      $p[num,true] = -(0.5*(-$omg - $knum)*(4*($y_size**2)-2) -
		       (- $omg + $knum))*
	exp(-0.5*($y_size**2))*sin($x_size[num]*$knum) 
      $u[num,true] = -(0.5*(-$omg - $knum)*(4*($y_size**2)-2) +
		       (- $omg + $knum))*
	exp(-0.5*($y_size**2))*sin($x_size[num]*$knum) 
      $v[num,true] = -($omg**2 - $knum**2)*2*$y_size*exp(-0.5*($y_size**2))*
	cos($x_size[num]*$knum)

    elsif $eqwave == "grav-east1" then

      $p[num,true] = -0.25*(0.5*(- $omg - $knum)*(4*($y_size**2)-2)-
			    (- $omg + $knum))*
	exp(-0.5*($y_size**2))*sin($x_size[num]*$knum) 
      $u[num,true] = -0.25*(0.5*(- $omg - $knum)*(4*($y_size**2)-2)+
			    (- $omg + $knum))*
	exp(-0.5*($y_size**2))*sin($x_size[num]*$knum) 
      $v[num,true] = -0.25*($omg**2 - $knum**2)*2*$y_size*
	exp(-0.5*($y_size**2))*cos($x_size[num]*$knum)

    elsif $eqwave == "grav-west1" then

      $p[num,true] = -(0.5*(-$omg - $knum)*(4*($y_size**2)-2) -
		       (-$omg + $knum))*
	exp(-0.5*($y_size**2))*sin($x_size[num]*$knum) 
      $u[num,true] = -(0.5*(-$omg - $knum)*(4*($y_size**2)-2) +
		       (-$omg + $knum))*
	exp(-0.5*($y_size**2))*sin($x_size[num]*$knum) 
      $v[num,true] = -($omg**2 - $knum**2)*2*$y_size*exp(-0.5*($y_size**2))*
	cos($x_size[num]*$knum)
      
    elsif $eqwave == "rossby2" then

      $p[num,true] = -(0.5*(-$omg - $knum)*(8*($y_size**3)-12*$y_size) - 
		       2*(-$omg + $knum)*2*$y_size)*
	exp(-0.5*($y_size**2))*sin($x_size[num]*$knum) 
      $u[num,true] = -(0.5*(-$omg - $knum)*(8*($y_size**3)-12*$y_size) +
		       2*(-$omg + $knum)*2*$y_size)*
	exp(-0.5*($y_size**2))*sin($x_size[num]*$knum) 
      $v[num,true] = -($omg**2 - $knum**2)*(4*($y_size**2)-2)*
	exp(-0.5*($y_size**2))*cos($x_size[num]*$knum)
      
    elsif $eqwave == "grav-east2" then
      
      $p[num,true] = -(0.5*(-$omg - $knum)*(8*($y_size**3)-12*$y_size) -
		       2*(-$omg + $knum)*2*$y_size)*
	exp(-0.5*($y_size**2))*sin($x_size[num]*$knum) 
      $u[num,true] = -(0.5*(-$omg - $knum)*(8*($y_size**3)-12*$y_size) +
		       2*(-$omg + $knum)*2*$y_size)*
	exp(-0.5*($y_size**2))*sin($x_size[num]*$knum) 
      $v[num,true] = - ($omg**2 - $knum**2)*(4*($y_size**2)-2)*
	exp(-0.5*($y_size**2))*cos($x_size[num]*$knum)
      
    elsif $eqwave == "grav-west2" then
      
      $p[num,true] = - (0.5*(-$omg - $knum)*(8*($y_size**3)-12*$y_size) -
			2*(-$omg + $knum)*2*$y_size)*
	exp(-0.5*($y_size**2))*sin($x_size[num]*$knum) 
      $u[num,true] = - (0.5*(-$omg - $knum)*(8*($y_size**3)-12*$y_size) +
			2*(-$omg + $knum)*2*$y_size)*
	exp(-0.5*($y_size**2))*sin($x_size[num]*$knum) 
      $v[num,true] = - ($omg**2 - $knum**2)*(4*($y_size**2)-2)*
	exp(-0.5*($y_size**2))*cos($x_size[num]*$knum)
      
    end

  }

  $p = VArray.new($p).rename($eqwave).
    put_att("units","1").put_att("ape_name",$title_hash[$eqwave])
  $p = GPhys.new( grid, $p)
  
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 }

=begin
  levels = NArray[-1.0,-0.6,-0.3, 0, 0.3, 0.6, 1.0]
  patterns = NArray[30999, 40999, 55999, 70999, 75999,85999]
  cont_lev, line_type, label  = t.cont_lev_set(levels)
  $tone_hash = { 'patterns'=> patterns, 'levels'=>levels }
    $cont_hash = {'levels' => cont_lev, 'index'=> [2,1], 'label'=> label, 
      'line_type'=> line_type }
  $cont_hash = {'int' => 0.15 }
  
  $cbar_conf = {
    "levels"=>levels, 
    "colors"=>patterns,
    "eqlev"=>true,
    "nobound"=>3,
    "tick1"=>20,"tick2"=>1
  }
=end


  plot_def($p)
  $ape.gropn(mkfignum)
  $ape.plot_main($p,$u,$v)
  $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 bunsan_calc

  ## 虚数を含む軸
  x_size = (NArray.complex(17).indgen - 8.0) *0.5
  
  if $eqwave == "kelvin" then
    
    $knum = 0.5
    $omg = $knum
    
  elsif $eqwave == "mixed-grav" then
    
    $knum = 0.5
    $omg =  ( $knum/2 - ( ($knum/2)**2 + 1 )**0.5 )
    
  elsif $eqwave == "mixed-rossby" then
    
    $knum = 1.0
    $omg =  ( $knum/2 - ( ($knum/2)**2 + 1 )**0.5 )

  elsif $eqwave == "grav-east0" then

    $knum = 0.5
    $omg = ( $knum/2 + ( ($knum/2)**2 + 1 )**0.5 )  
    
  elsif $eqwave == "rossby1" then
    
    $knum = 0.5
    t1, t2 = caldun(1)
    $omg = (t2 * Complex.polar(1, Math::PI*2/3) + 
	     t1 * Complex.polar(1, Math::PI*4/3) ).real
    
  elsif $eqwave == "grav-east1" then
    
    $knum = 0.5
    t1, t2 = caldun(1)
    $omg = (t1 + t2).real

  elsif $eqwave == "grav-west1" then

    $knum = 0.5
    t1, t2 = caldun(1)
    $omg = (t1 *  Complex.polar(1, Math::PI*2/3) + 
             t2 * Complex.polar(1, Math::PI*4/3)).real
    
  elsif $eqwave == "rossby2" then

    $knum = 0.5
    t1, t2 = caldun(2)
    $omg = (t2 * Complex.polar(1, Math::PI*2/3) + 
	   t1 * Complex.polar(1, Math::PI*4/3) ).real
    
  elsif $eqwave == "grav-east2" then

    $knum = 0.5
    t1, t2 = caldun(2)
    $omg = (t1 + t2).real
    
  elsif $eqwave == "grav-west2" then
    
    $knum =  0.5    
    t1, t2 = caldun(2)
    $omg = (t1 *  Complex.polar(1, Math::PI*2/3) + 
	   t2 * Complex.polar(1, Math::PI*4/3)).real
    
  end

end


def caldun(num)
  
  ## 3 次方程式解法: カルダノの公式
  ## x^3 + 3px + q
  ## t^2 + qt - p =0, t = -q/2 +- (q**2/4 + p)**0.5
  p = - ( ($knum)**2 + (2*num+1)) / 3
  p = Complex(p,0)
  q = - $knum
  q = Complex(q,0)
  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)
  
  return  t1, t2

end


# トーン, コンターパターンのデフォルト設定
def plot_def(gphys)
  
  # トーンパターンのデフォルト
  patterns = NArray[30999, 40999, 55999, 70999, 75999, 85999]
  
  # トーンレベルのデフォルト: ((最大値 - 最小値) / 7 )
  levels = Array.new
  dx = (gphys.max - gphys.min)/6 ;  x  = gphys.min 
  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