#!/usr/bin/env ruby
require "numru/gphys"
include NumRu

dir = "."
varname = "sd"
# number of application of 1-2-1 filter
smoothing_ntimes0 = 40
smoothing_ntimes1 = 10

def smoothing_in_freq(na0,smoothing_ntimes)
  na1 = NArray.sfloat(*na0.shape)
  n = na0.shape[1]
  smoothing_ntimes.times do |i|
    # 1-2-1 filter: はしっこは 2-1 filter, 1-2 filterにする
    na1[true, 0] = (2*na0[true, 0] + na0[true, 1])/3.0
    (1..n-2).each do |j|
      na1[true, j] = (na0[true, j-1] + 2*na0[true, j] + na0[true, j+1])/4.0
    end
    na1[true, n-1] = (na0[true, n-2]+ 2*na0[true, n-1])/3.0
    na0 = na1
  end
  na0
end

def smoothing_in_k(na0,smoothing_ntimes)
  na1 = NArray.sfloat(*na0.shape)
  n = na0.shape[0]
  smoothing_ntimes.times do |i|
    # 1-2-1 filter: はしっこは 2-1 filter, 1-2 filterにする
    na1[0,true] = (2*na0[0,true] + na0[1,true])/3.0
    (1..n-2).each do |j|
      na1[j,true] = (na0[j-1,true] + 2*na0[j,true] + na0[j+1,true])/4.0
    end
    na1[n-1,true] = (na0[n-2,true]+ 2*na0[n-1,true])/3.0
    na0 = na1
  end
  na0
end

gp_sp_s, gp_sp_a = ["sym", "antisym"].map do |t|
  file = "sp-#{t}.nc"
  GPhys::IO.open_gturl("#{dir}/#{file}@#{varname}").copy
end
gp_sp = gp_sp_s + gp_sp_a

na0 = (gp_sp.val.class == NArrayMiss ? gp_sp.val.to_na(0) : gp_sp.val )
na0 = smoothing_in_freq(na0, smoothing_ntimes1)
na0 = smoothing_in_k(na0, smoothing_ntimes0)
gp_sp.val = na0

outfile = NetCDF.create("sp-bg.nc")
GPhys::IO.write(outfile, gp_sp)
outfile.close