require "staggered_deriv.rb"
require "numru/ggraph"
include NumRu

def d_x(gp) #将来的には四角いmodel用moduleかclassで定義される。
  bc = gp.axis(0).get_aux("bc").val[0]
  GPhys::Stg_deriv::cderiv(gp, 0, 1, bc)
end

Grav = UNumeric.new(9.8, "m.s-2")
Depth = UNumeric.new(1.0, "m")
Dt = UNumeric.new(10.0,"s")
Dx = 100.0
Nx = 100
Nloop = 10000
XMAX = Dx * Nx

#<< make GPhys>>
pos_bounds = NArray.new("float", Nx+1).indgen(0.0, Dx)
pos_bounds = VArray.new(pos_bounds, {"long_name"=>"X coordinate","units"=>"m"},
"x_b")
pos_center = NArray.new("float", Nx).indgen(Dx/2.0, Dx)
pos_center = VArray.new(pos_center, {"long_name"=>"X coordinate","units"=>"m"}, "x")

axis_x = Axis.new(cell=true).set_cell( pos_center, pos_bounds )
#h_x_coord = Axis.new(cell=true).set_cell( pos_center, pos_bounds ).set_pos_to_center
#u_x_coord = Axis.new(cell=true).set_cell( pos_center, pos_bounds ).set_pos_to_bounds

h_x_coord = axis_x.copy.set_pos_to_center
u_x_coord  = axis_x.copy.set_pos_to_bounds

#<< initial wave >>
init_h = NArray.sfloat(Nx)
for i in 0..Nx-1
  init_h[i] = 0.5 * ( Math.cos((4.0 * i / (Nx-1) - 1.0) * Math::PI) + 1.0 )
end
init_h[Nx/2..-1] = 0

data_h = VArray.new( init_h, 
                   {"long_name"=>"depth perturbation", "units"=>"m"},
                   "h" )

data_u = VArray.new( NArray.sfloat(Nx+1),
                   {"long_name"=>"X velocity", "units"=>"m/s"},
                   "u" )

gp_h = GPhys.new(Grid.new(h_x_coord), data_h)
gp_u = GPhys.new(Grid.new(u_x_coord), data_u)


# Axis.set_bcをつくる？？
na = NArray.object(1)
na[0] = "asym"
na[0] = "cyclic"
u_bc = VArray.new( na, nil, "bc_u" )
na = NArray.object(1)
na[0] = "sym"
na[0] = "cyclic"
h_bc = VArray.new( na, nil, "bc_h" )

gp_u.axis(0).set_aux("bc", u_bc)
gp_h.axis(0).set_aux("bc", h_bc)
# ここまで 

gp_u_a = gp_u.copy
gp_u_b = gp_u.copy

gp_h_a = gp_h.copy
gp_h_b = gp_h.copy


#<<intial integrate>>

gp_u = gp_u_b - Grav * d_x(gp_h_b) * Dt
gp_h = gp_h_b - Depth * d_x(gp_u_b) * Dt

DCL.swpset('lalt',true)   
DCL.gropn(1)
DCL.sldiv('y',1,2) 
DCL.sgpset('lcntl', false)   # 制御文字を解釈しない
DCL.sgpset('lfull',true)     # 全画面表示
DCL.sgpset('lcorner',false)   
DCL.uzfact(0.6)             # 座標軸の文字列サイズを 0.75 倍
DCL.sgpset('lfprop',true)    # プロポーショナルフォントを使う
GGraph.set_fig('viewport'=>[0.1,0.95,0.1,0.25])  

#<< main loop >>
for itr in 1..Nloop
  time = Dt * itr

  gp_u_a = gp_u_b - Grav * d_x(gp_h) * Dt * 2.0
  gp_h_a = gp_h_b - Depth * d_x(gp_u) * Dt * 2.0

  GGraph.set_axes( 'xside'=>'tb', 'yside'=>'lr', 'xtitle'=>"")
  GGraph.set_fig("window"=>[0.0, XMAX, -1.0, 1.0])
  GGraph.line(gp_u_a)
  GGraph.line(gp_h_a)

  gp_u_b = gp_u 
  gp_u   = gp_u_a

  gp_h_b = gp_h
  gp_h   = gp_h_a
  
end

DCL.grcls