#!/usr/bin/env ruby

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

#$local_path = '/work11/ape/yukiko/lib' 
$local_path = '/GFD_Dennou_Work2/yukiko/eva01/work11/ape/yukiko/lib' 
$: << $local_path

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

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

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


END{

$rezol = ["T39L48_eml",
  "T39L24_eml",
  "T39L96_eml",
  "T79L48_eml" ]

rezolid = [
  "T39L24_eml",
  "T39L96_eml",
  "T79L48_eml",
  "T39L24_non",
  "T39L96_non",
  "T79L48_non"
 ]


# $rezol = ["T79L48_non"]
#$rezol = ["T39L48_eml"]
#$resol = "T39L48_eml"
#$resol = "T159L48_non"
#$resol = "T159L48_eml"
#$resol = "T319L48_non"
#$sstid = "control"

 sstid_all = ["control","flat","peaked","1keq","3keq","3kw1","Qobs","control-5N","flat3keq","Qobs3keq","H1998con","H1998pa"]
 sstid = ["control","flat","peaked","1keq","3keq","3kw1","Qobs","control-5N"]
sstid_tmp = ["peaked", "control-5N", "1keq", "3kw1"]
 sstid_tmp = ["control","flat","3keq","Qobs","flat3keq","Qobs3keq","H1998con","H1998pa"]


cumulusid = ["ias","ksc","kuo","mca","eml","non"]
sstid_all = ["control","flat","peaked","1keq","3keq","3kw1","Qobs","control-5N","flat3keq","Qobs3keq","H1998con","H1998pa"]

#  cumulusid.each{ |item|
#  rezolid.each{ |item|
  $rezol.each{ |item|
#    sstid.each{ |sst|
#      $resol = "T39L48_#{item}"
      $resol = item
      $sstid = "control"
      mknc_stspctrum
    }
#  }

}


def mknc_stspctrum

  $ncfile_path = "/work11/ape/NetCDF/#{$resol}/"
  $groupid = "AGUforAPE-03a"
  $expID = $sstid
   
  $t= ape_new "/GFD_Dennou_Work2/yukiko/eva01/work11/ape/NetCDF/#{$resol}/AGUforAPE-03a_TR_#{$sstid}.nc"    
#  $t= ape_new "/work11/ape/NetCDF/#{$resol}/AGUforAPE-03a_TR_#{$sstid}.nc"    
#  $t= ape_new "/work11/ape/NetCDF/#{$resol}/AGUforAPE-03a_TR_#{$sstid}_half1.nc"    
#   $t = Ape.new("/GFD_Dennou_Work2/yukiko/eva01/work11/ape/yukiko/data/T159L48_eml_expID01/T159L48_eml_expID01_0013.nc")

  file = NetCDF.create("#{$resol}_#{$sstid}_wk_40smooth.nc")
  
  puts "#{$resol}_#{$sstid}_wk_40smooth.nc create start"

  $t.netcdf_open.var_names.each { |name| 
    if name =~ /lw_toa|mslp|500|250|ps|tppn/
#    if name =~ /tppn/
#    if name =~ /OMG|V|U|mslp|tppn/
    
      puts name
      gphys = $t.go(name)
#      gphys = __tr_T159_marge(name)
#      gphys = __tr_T159_eml_marge(name)
#      gphys = __tr_T319_marge(name)
      

      gphys_sym, gphys_asym, gphys_sym_wk, gphys_asym_wk, gphys_bg = spct_wk1999(gphys)
      dim1 = gphys_sym.coord(0).shape.to_s.to_i
      dim2 = gphys_sym.coord(1).shape.to_s.to_i
      
      # netCDF 初期化
      GPhys::NetCDF_IO.write(file, gphys_sym_wk )
      GPhys::NetCDF_IO.write(file, gphys_asym_wk )
      GPhys::NetCDF_IO.write(file, gphys_bg)
      GPhys::NetCDF_IO.write(file, gphys_sym)
      GPhys::NetCDF_IO.write(file, gphys_asym)
      
    end
  }      

      # 大域属性
      file.put_att("Conventions", "CF-1.0")
      file.put_att("title","Aqua Planet: #{$resol}_#{$sstid} Experiment, Wheeler and Kiladis, 1999 Plot")
      file.put_att("history", "Original data produced: #{Time.now} from TR netCDF")
      file.put_att("institution", "AGU for APE")
      file.put_att("source","afes_115_ape_#{$cumulus}.LM")
      file.close
  
end

# --------------------------------------------------------------------
  
def __tr_T159_marge(name)
  
  file1 = Ape.new("/work11/ape/NetCDF/#{$resol}/AGUforAPE-03a_TR_#{$sstid}_half1.nc")
  file2 = Ape.new("/work11/ape/NetCDF/#{$resol}/AGUforAPE-03a_TR_#{$sstid}_half2.nc")

  # data
  data = NArray.sfloat(480,40,1440).indgen * 0.0
  
  gphys = file1.gphys_open(name)
  gphys = gphys.cut(true,-15..15,true)
  data[true,true,0..719] = gphys.data.val
  
  gphys = file2.gphys_open(name)
  gphys = gphys.cut(true,-15..15,true)
  data[true,true,720..-1] = gphys.data.val
  
  data = VArray.new(data).rename("#{name}").
    put_att("units",gphys.data.get_att("units")).
    put_att("ape_name",gphys.data.get_att("ape_name"))


  # 軸
#  x_size = NArray.sfloat(480).indgen*360.0/480
  t_size = NArray.sfloat(1440).indgen*1.0/4.0 + 0.25
  y_size = gphys.grid_copy.axis(1).pos.val
  x_size = gphys.grid_copy.axis(0).pos.val
  
  x = VArray.new(x_size).rename("lon").put_att("units","degrees_east")
  y = VArray.new(y_size).rename("lat").put_att("units","degrees_north")
  t = VArray.new(t_size).rename("time").put_att("units","days since 0000-01-01")
  x = Axis.new().set_pos(x)
  y = Axis.new().set_pos(y)
  t = Axis.new().set_pos(t)
  grid = Grid.new(x,y,t)

  data = GPhys.new( grid, data)
  
  return data
end
  
def __tr_T159_eml_marge(name)


    file1 = Ape.new("/GFD_Dennou_Work2/yukiko/eva01/work11/ape/yukiko/data/T159L48_eml_expID01/T159L48_eml_expID01_0011.nc")
    file2 = Ape.new("/GFD_Dennou_Work2/yukiko/eva01/work11/ape/yukiko/data/T159L48_eml_expID01/T159L48_eml_expID01_0012.nc")

    file3 = Ape.new("/GFD_Dennou_Work2/yukiko/eva01/work11/ape/yukiko/data/T159L48_eml_expID01/T159L48_eml_expID01_0013.nc")
    file4 = Ape.new("/GFD_Dennou_Work2/yukiko/eva01/work11/ape/yukiko/data/T159L48_eml_expID01/T159L48_eml_expID01_0014.nc")

  # data
  data = NArray.sfloat(480,40,1440).indgen * 0.0
  
  gphys = file1.gphys_open(name)
  gphys = gphys.cut(true,-15..15,true)
  data[true,true,0..359] = gphys.data.val
  
  gphys = file2.gphys_open(name)
  gphys = gphys.cut(true,-15..15,true)
  data[true,true,360..719] = gphys.data.val
  
  gphys = file3.gphys_open(name)
  gphys = gphys.cut(true,-15..15,true)
  data[true,true,720..1079] = gphys.data.val
  
  gphys = file4.gphys_open(name)
  gphys = gphys.cut(true,-15..15,true)
  data[true,true,1080..-1] = gphys.data.val

  data = VArray.new(data).rename("#{name}").
    put_att("units",gphys.data.get_att("units")).
    put_att("ape_name",gphys.data.get_att("ape_name"))


  # 軸
#  x_size = NArray.sfloat(480).indgen*360.0/480
  t_size = NArray.sfloat(1440).indgen*1.0/4.0 + 0.25
  y_size = gphys.grid_copy.axis(1).pos.val
  x_size = gphys.grid_copy.axis(0).pos.val
  
  x = VArray.new(x_size).rename("lon").put_att("units","degrees_east")
  y = VArray.new(y_size).rename("lat").put_att("units","degrees_north")
  t = VArray.new(t_size).rename("time").put_att("units","days since 0000-01-01")
  x = Axis.new().set_pos(x)
  y = Axis.new().set_pos(y)
  t = Axis.new().set_pos(t)
  grid = Grid.new(x,y,t)

  data = GPhys.new( grid, data)
  
  return data
end

# --------------------------------------------------------------------
  
def __tr_T319_marge(name)
   

  file1 = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_0031.nc")
  file2 = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_0032.nc")
  file3 = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_0033.nc")
  file4 = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_0034.nc")
  file5 = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_0035.nc")    
  file6 = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_0036.nc")
  file7 = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_0037.nc")
  file8 = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_0038.nc")
  file9 = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_0039.nc")
  file10 = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_0040.nc")
  file11 = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_0041.nc")
  file12 = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_0042.nc")
  
    
  # data
  data = NArray.sfloat(960,80,1440).indgen * 0.0
  
  
  gphys = file1.gphys_open(name).cut(true,-15..15,true)
  data[true,true,0..119] = gphys.data.val
  gphys = file2.gphys_open(name).cut(true,-15..15,true)
  data[true,true,120..239] = gphys.data.val
  gphys = file3.gphys_open(name).cut(true,-15..15,true)
    data[true,true,240..359] = gphys.data.val
  gphys = file4.gphys_open(name).cut(true,-15..15,true)
  data[true,true,360..479] = gphys.data.val
  gphys = file5.gphys_open(name).cut(true,-15..15,true)
  data[true,true,480..599] = gphys.data.val
  gphys = file6.gphys_open(name).cut(true,-15..15,true)
  data[true,true,600..719] = gphys.data.val
  gphys = file7.gphys_open(name).cut(true,-15..15,true)
  data[true,true,720..839] = gphys.data.val
  gphys = file8.gphys_open(name).cut(true,-15..15,true)
  data[true,true,840..959] = gphys.data.val
  gphys = file9.gphys_open(name).cut(true,-15..15,true)
  data[true,true,960..1079] = gphys.data.val
  gphys = file10.gphys_open(name).cut(true,-15..15,true)
  data[true,true,1080..1199] = gphys.data.val
  gphys = file11.gphys_open(name).cut(true,-15..15,true)
  data[true,true,1200..1319] = gphys.data.val
  gphys = file12.gphys_open(name).cut(true,-15..15,true)
  data[true,true,1320..1439] = gphys.data.val
  
  data = VArray.new(data).rename("#{name}").
    put_att("units",gphys.data.get_att("units")).
    put_att("ape_name",gphys.data.get_att("ape_name"))
  
  # 軸
  t_size = NArray.sfloat(1440).indgen*1.0/4.0 + 0.25
  y_size = gphys.grid_copy.axis(1).pos.val
  x_size = gphys.grid_copy.axis(0).pos.val
  
  x = VArray.new(x_size).rename("lon").put_att("units","degrees_east")
  y = VArray.new(y_size).rename("lat").put_att("units","degrees_north")
  t = VArray.new(t_size).rename("time").put_att("units","days since 0000-01-01")
  x = Axis.new().set_pos(x)
  y = Axis.new().set_pos(y)
  t = Axis.new().set_pos(t)
  grid = Grid.new(x,y,t)

  data = GPhys.new( grid, data)
  
  return data
end


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

def spct_wk1999(gphys)
  
  gphys1 = gphys.cut(true,-15..0,true)
  gphys2 = gphys.cut(true,0..15,true)

#  gphys1 = gphys.cut(true,-2..0,true)
#  gphys2 = gphys.cut(true,0..2,true)

  dim = gphys1.coord(1).shape.to_s.to_i  
  work1 = gphys1[true,0,0..359].stspct_fft("spct") * 0.0
  work2 = gphys2[true,0,0..359].stspct_fft("spct") * 0.0

  dim0 = gphys1.coord(0).shape.to_s.to_i  # 緯度
  dim1 = gphys1.coord(1).shape.to_s.to_i  # 経度

=begin
  dim.times { |dim_num|
    10.times { |num|
      work1 = gphys1[true,dim_num,(num*120)..359+(num*120)].
	stspct_fft("spct") + work1
      work2 = gphys2[true,dim_num,(num*120)..359+(num*120)].
	stspct_fft("spct") + work2
    }
  }
=end


  # 最小二乗法 + hamming taper + 時空間スペクトル足しあわせ
  10.times { |num|

    print "#{num}\n"

    g1 = gphys1[true,true,(num*120)..359+(num*120)].copy
    g2 = gphys2[true,true,(num*120)..359+(num*120)].copy
    	
    x   = NArray.sfloat(dim0,dim1).fill!(0.0)
    xx  = NArray.sfloat(dim0,dim1).fill!(0.0)
    xy1 = NArray.sfloat(dim0,dim1).fill!(0.0)
    xy2 = NArray.sfloat(dim0,dim1).fill!(0.0)
    y1 =  NArray.sfloat(dim0,dim1).fill!(0.0)
    y2 =  NArray.sfloat(dim0,dim1).fill!(0.0)
    
    time =  NArray.sfloat(360).indgen!
    
    time.each{ |time_num|

      xy1 = time_num * g1.data.val[true,true,time_num] + xy1
      xy2 = time_num * g2.data.val[true,true,time_num] + xy2
      x  = time_num + x  
      y1  = g1.data.val[true,true,time_num] + y1
      y2  = g2.data.val[true,true,time_num] + y2
      xx = time_num + xx
      
    }

    a1 = (360 * xy1 - x * y1  ) / (360 * xx - x**2 )
    a2 = (360 * xy2 - x * y2  ) / (360 * xx - x**2 )
    b1 = ( xx * y1  - x * xy1 ) / (360 * xx - x**2 )
    b2 = ( xx * y2  - x * xy2 ) / (360 * xx - x**2 )

    time.each{ |time_num|
      g1[true,true,time_num] = ( g1.data.val[true,true,time_num] - ( a1 * time_num + b1 ) ) * 
	( 1.0 - cos(2*Math::PI*(time_num - 0.5)/360.0) ) / 2.0 
      g2[true,true,time_num] = ( g2.data.val[true,true,time_num] - ( a2 * time_num + b2 ) ) * 
	( 1.0 - cos(2*Math::PI*(time_num - 0.5)/360.0) ) / 2.0 
    }

    dim1.times { |dim1_num|
      work1 = g1[true,dim1_num,true].stspct_fft("spct") + work1
      work2 = g2[true,dim1_num,true].stspct_fft("spct") + work2
    }

  }

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

  gphys_asym = (( work1 - work2 )/2 ).
    set_att("ape_name", "#{gphys.data.get_att("ape_name")} S-T power spectrum").
    set_att("units", "1" ). 
    rename("#{gphys.name}_asym_org_spct")
  gphys_sym  = (( work1 + work2 )/2 ).
    set_att("ape_name", "#{gphys.data.get_att("ape_name")} S-T power spectrum").
    set_att("units", "1" ). 
    rename("#{gphys.name}_sym_org_spct")

  gphys = ( gphys_sym + gphys_asym ).rename("#{gphys.name}")
  gphys_bg  = spct_121mean(gphys)

  gphys_sym_wk = spct_divide(gphys_sym,gphys_bg).
    set_att("ape_name", gphys.data.get_att("ape_name")).
    set_att("units","1"). 
    rename("#{gphys.name}_sym_spct")
  gphys_asym_wk = spct_divide(gphys_asym,gphys_bg).
    set_att("ape_name", gphys.data.get_att("ape_name") ).
    set_att("units","1"). 
    rename("#{gphys.name}_asym_spct")

  gphys_bg = gphys_bg.rename("#{gphys.name}_bg_spct").
    set_att("ape_name", gphys.data.get_att("ape_name")).
    set_att("units", "1" )

  return gphys_sym, gphys_asym, gphys_sym_wk, gphys_asym_wk, gphys_bg
end


def spct_121mean(gphys)
  
  dim1 = gphys.coord(0).shape.to_s.to_i
  dim2 = gphys.coord(1).shape.to_s.to_i
  
  gphys_bg = gphys.copy
  work     = gphys.copy

  40.times{ |num|
#  10.times{ |num|
    dim1.times { |dim1_num|
      if dim1_num == 0
	work[dim1_num,true] = 
	  ( gphys_bg[dim1_num,true] + gphys_bg[dim1_num+1,true] )/2 
      elsif dim1_num == (dim1-1)
	work[dim1_num,true] = 
	  ( gphys_bg[dim1_num-1,true] + gphys_bg[dim1_num,true] )/2 
      else
	work[dim1_num,true] = 
	  ( gphys_bg[dim1_num-1,true]+ gphys_bg[dim1_num,true] + gphys_bg[dim1_num,true] + gphys_bg[dim1_num+1,true] )/4 
      end
    }
    
    dim2.times { |dim2_num|
      if dim2_num == 0
	gphys_bg[true, dim2_num] = 
	  ( work[true,dim2_num] + work[true,dim2_num+1] )/2 
      elsif dim2_num == (dim2-1)
	gphys_bg[true, dim2_num] = 
	  ( work[true,dim2_num-1] + work[true,dim2_num] )/2 
      else
	gphys_bg[true, dim2_num] = 
	  ( work[true,dim2_num-1]+ work[true,dim2_num] + work[true,dim2_num] +  work[true,dim2_num+1] )/4 
      end
    }
  }    


  return gphys_bg

end

def spct_divide(gphys,gphys_bg)  

  dim1 = gphys.coord(0).shape.to_s.to_i
  dim2 = gphys.coord(1).shape.to_s.to_i

  dim1.times{ |dim1_num|
    dim2.times{ |dim2_num|
      unless gphys_bg[dim1_num,dim2_num].data.val == 0.0
	gphys[dim1_num,dim2_num]=
	  gphys[dim1_num,dim2_num].data.val/gphys_bg[dim1_num,dim2_num].data.val
      else 
	gphys[dim1_num,dim2_num] = 0.0 
      end
    }
  }

  return gphys
end




=begin
    $t.gropn(1)

    $t.plot(gphys_sym[((dim1+1)/2-31)..((dim1+1)/2+29),0..80])
    file_name = "#{$file_label}_#{gphys_sym.name}.gif"
    p $file_name

    # ダンプする場合
    GGraph.tone( gphys_sym[((dim1+1)/2-31)..((dim1+1)/2+29),0..80])
    `for file in *.xwd ; do xwdtopnm ${file} | pnmcut 2 2 904 654 > tmp.pnm; done`
    `ppmtogif tmp.pnm > #{$fig_path}#{file_name}`
    `rm dcl_* tmp.pnm` 

    $t.plot(gphys_asym[((dim1+1)/2-31)..((dim1+1)/2+29),0..80])
    file_name = "#{$file_label}_#{gphys_asym.name}.gif"
    p $file_name
    
    # ダンプする場合
    GGraph.tone( gphys_asym[((dim1+1)/2-31)..((dim1+1)/2+29),0..80])
   `for file in *.xwd ; do xwdtopnm ${file} | pnmcut 2 2 904 654 > tmp.pnm; done`
    `ppmtogif tmp.pnm > #{$fig_path}#{file_name}`
    `rm dcl_* tmp.pnm` 

  }
  
  # おしまい
  $t.grcls
  `rm tmp.pnm dcl_*`
=end

=begin
  dim.times { |dim_num|
    5.times { |num|
      work1 = gphys1[true,dim_num,(num*240)..359+(num*240)].
	stspct_fft("spct") + work1
      work2 = gphys2[true,dim_num,(num*240)..359+(num*240)].
	stspct_fft("spct") + work2
    }
  }
=end