#
# 表題: 時空間スペクトル(100日)生成
#
# 履歴: 2003/07/05 やまだ由
#
# メモ:
#
# u-z13-eq-spct-con.ps
# u-z12-eq-spct-con.ps
# v-z12-eq-spct-con.ps
# rain-sum-eq-spct-con-sym.ps
# rain-sum-eq-spct-con-asym.ps
#
# の生成に利用. 
#
# gtool3 形式 -(gt3conv)-> NetCDF形式 --> 
# gdcl である程度成形, コードをセーブ
#
# カラーバーがうまくいってない. 
# プログラム内の説明書きがなにもない. 
# 配列の書き方がよく分からん. 
#
#
# gtool3 コマンドメモ
#
# gtsel u.out z=13  y=32 out:stp.out
# gtstspct stp.out out:stp.out 
# gtext ystr=1 yend=50 xstr=-30 xend=30 stp.out out:u-z13-eq.spct
# gtcont tone=,0.001,0.01,0.1,1,10,100 pat=25999,35999,45999,55999,70999,75999,85999
#
# gif ファイル生成
#
# pstopnm dcl.ps 
# ppmtogif dcl001.ppm >  dcl.gif
# gifsicle --rotate-90 dcl.gif > rain-sum-eq-spct-con-asym.gif
# mv dcl.ps rain-sum-eq-spct-con-asym.ps 
#


require "gtk"
require "narray"
require "numru/dcl"
require "numru/netcdf"

include NumRu


filename = "/home/yukiko/tmp/inshi/ruby/kuo/u_z12_y30_kuo.nc"
#filename = "/home/yukiko/tmp/inshi/ruby/kuo/v_z12_eq_kuo.nc"
#filename = "/home/yukiko/tmp/inshi/ruby/kuo/rain_sum_eq_kuo.nc"


file = NetCDF.open(filename,mode="r")

varname = "u"
#varname = "v"
#varname = "rain"
varobj = file.var(varname)
st = [0,0,0]; en = [-1,-1,0]
var = varobj.get("start"=>st,"end"=>en)

axisxname = "wvnwe"
axisxobj = file.var(axisxname)
axisx = axisxobj.get

axisyname = "freq200_100"
axisyobj = file.var(axisyname)
axisy = axisyobj.get

file.close

xx = NArray.sfloat(16).indgen!
 
h = [ 12, 25, 50, 100, 200 ]
hh = NArray.to_na(h)

c = ((hh*9.8)**0.5)*60*60*24/(4e+7)

kelvin = [xx,xx,xx,xx,xx]
rossby = [xx,xx,xx,xx,xx]
grab   = [xx,xx,xx,xx,xx]
mixa  = [xx,xx,xx,xx,xx]
mixb  = [xx,xx,xx,xx,xx]

for i in (0..4)
  kelvin[i] = xx*c[i]
  rossby[i] = -xx/(xx**2+(2 +1)/c[i] )
  grab[i] = -(c[i]*( c[i]*xx**2 + (2+1) ) )**0.5
  mixa[i] = c[i]*xx/2 + ((c[i]*xx/2)**2 + c[i] )**0.5 
  mixb[i] = c[i]*xx/2 - ((c[i]*xx/2)**2 + c[i] )**0.5 
end


DCL::gropn(1)
DCL::sglset("LCNTL", false )
DCL::udlset("LMSG", false )
DCL::gllset("LMISS", true )
DCL::glrset("RMISS", -999.0)

xmin = -15.0
xmax = 15.0
#ymin = 0.009999999776
ymin = 0.0
ymax = 0.5

vxmin = 0.15
vxmax = 0.85
vymin = 0.15
vymax = 0.65


DCL::sglset('LFULL', true)
DCL::sglset('LCORNER', false)
DCL::grfrm
DCL::grswnd(xmin,xmax,ymin,ymax)
DCL::grsvpt(vxmin,vxmax,vymin,vymax)
DCL::grstrn(1)
DCL::grstrf

x_title = "wavenumber"
x_unit = "1"
y_title = "frequency"
y_unit = "1/DAY"

DCL::sglset("LCLIP", true )

DCL::uwsgxa(axisx)
DCL::uwsgya(axisy)
DCL::ueitlv
DCL::uegtla(0.0, 0.049, 0.001 )

# u,  v
DCL::uestlv(0,0.001,01000)  
DCL::uestlv(0.001,0.01,40000)  
DCL::uestlv(0.01,0.1,60999)  
DCL::uestlv(0.1,1,70999)
DCL::uestlv(1,10,90999)  

# rain
#DCL::uestlv(0,1,01000)  
#DCL::uestlv(1,10,40000)  
#DCL::uestlv(10,100,60999)  
#DCL::uestlv(100,1000,70999)
#DCL::uestlv(1000,10000,90999)  


DCL::uetonf(var)

#DCL::udgcla(1.036227246e-07, 6.172527313, 0 )
#DCL::udcntz(var)

for i in (0..4)

  DCL::uuslnt(1)
  ii=(i)*20+5
  DCL::uuslni(ii)
  DCL::uulin(xx, kelvin[i])
  DCL::uulin(-xx, -rossby[i])
  DCL::uulin(-xx, -grab[i])
  DCL::uulin(xx, -grab[i])
#  DCL::uulin(xx,mixa[i])
#  DCL::uulin(-xx,-mixb[i])

end

DCL::ussttl(x_title, x_unit, y_title, y_unit)
DCL::usdaxs
DCL::sglset("LCLIP", false )

i=0
#DCL::sgtxzv(vxmax+0.01,vymax-0.03*i,"time=1584000.0",0.02,0,-1,1)
DCL::sgtxzv(vxmax+0.01,vymax-0.03*i,"",0.02,0,-1,1)
DCL::uzrset("ROFFXT", 0.06)
title = "u-velocity s-t"
DCL::uxsttl("t", title, 0 )


DCL::grswnd(0.0,1.0,0,10.0) 
DCL::grsvpt(0.88,0.90,0.40,0.65)  
DCL::grstrn(1)  
DCL::grstrf
#pi=[[0,1],[0,1,2,3,4,5]]
pi = NArray.sfloat(2, 6)
pi[0,5] = 0
pi[1,5] = 0
pi[0,4] = 0.001
pi[1,4] = 0.001
pi[0,3] = 0.01
pi[1,3] = 0.01
pi[0,2] = 0.1
pi[1,2] = 0.1
pi[0,1] = 1
pi[1,1] = 1
pi[0,0] = 10
pi[1,0] = 10

DCL::uestlv(0,0.001,01000)  
DCL::uestlv(0.001,0.01,40000)  
DCL::uestlv(0.01,0.1,60999)  
DCL::uestlv(0.1,1,70999)
DCL::uestlv(1,10,90999)  

#DCL::udcntr(pi)  
#DCL::uetonf(pi)  
#DCL::slpvpr(3)
##DCL::uzfact(0.8)
##DCL::uysfmt('(F4.1)')
#DCL::uyaxdv('R', 0.1, 0.2)



DCL::grcls