#!/usr/bin/env ruby
# -*- coding: utf-8 -*-
# 
# 1 次元定常移流方程式
# \frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = 0
#
# 時間方向に 1 次の前進差分, 空間方向に 2 次の中央差分を取った場合.
# 
# = 履歴 (編集した際は上に追記すること)
# * 2012/07/25 納多哲史  解析解を一緒に表示するよう変更
# * 2012/06/12 納多哲史  作成
# 

require "numru/ggraph"
include NumRu
include GGraph

### パラメータ設定
nt     = 9000      # 全ステップ数
nt_out = 300       # 出力間隔

m  = 128           # 格子点数
l  = 2 * Math::PI  # 系の大きさ
c  = 1             # 移流速度
dx = l / m         # 格子の間隔
dt = 0.001         # 時間ステップ

### 座標の設定
x = dx * NArray.float(m).indgen

### 初期値の設定
u0 = NArray.float(m).fill(0.0)
n1 = m / 8
n2 = m / 4
u0[n1..n2] = 1.0

########################################
# メソッド定義
########################################

class NArray
  ### 周期境界条件で i (整数) 個シフトさせた配列を返す
  def sft(i)
    na = NArray.float(self.size)
    na[0..(-i-1)] = self[i..-1]
    na[-i..-1] = self[0..(i-1)]
    na
  end
  ### 周期境界条件で xsft (浮動小数点) 個シフトさせた配列を返す.
  ### 解析解の計算用.
  def sft_float(xsft, l, m, dx)
    na = NArray.sfloat(self.size)
    sft = (xsft % l) / dx
    dm = sft.divmod(1.0)
    isft = dm[0]
    sft_mod = dm[1]
    na = (1.0 - sft_mod) * self.sft(isft) + sft_mod * self.sft((isft+1) % m) 
    na
  end
end

### 1 次元の NArray からおまかせで GPhys を作成
def na2gphys(nary, x)
  # option
  axis_opt = {'long_name' => 'X-axis', 'units' => '1', 'name' => 'x'}
  data_opt = {'long_name' => 'u', 'units' => '1', 'name' => 'u'}

  # 軸情報の設定
  va_axis_x = VArray.new( x,
                   {"long_name"=>axis_opt['long_name'],
                     "units"=>axis_opt['units']},
                   axis_opt['name'] )
  axis_x = Axis.new.set_pos(va_axis_x)

  # データとその付属情報の設定
  va_data = VArray.new( nary,
                   {"long_name"=>data_opt['long_name'],
                     "units"=>data_opt['units']},
                   data_opt['name'] )

  # GPhys オブジェクトの作成
  gphys = GPhys.new( Grid.new(axis_x), va_data )
end

### 解析解の計算
def u_ref(u0, xsft, l, m, dx)
  u_ref = u0.sft_float(xsft, l, m, dx)
  u_ref
end

### 描画
def draw(u, x, t, u0, l, m, dx, c)
  # gphys 形式へ変換 (NetCDF 形式で出力するときに便利なように)
  gphys = na2gphys(u, x)
  gphys_ref = na2gphys(u_ref(u0, -c*t, l, m, dx), x)
  # 折れ線グラフとして可視化
  # 解析解の描画
  line(gphys_ref, true, 'max'=>2.0, 'min'=>-1.0, 'index'=>15, 'title'=>"t="+t.to_s)
  # 数値計算によって得られた値の描画
  line(gphys, false, 'max'=>2.0, 'min'=>-1.0, 'index'=>25)
end

########################################
# メイン
########################################

### ウィンドウを開く
DCL.gropn(4)

### 初期値の代入と描画
u = u0
draw(u, x, 0.0, u0, l, m, dx, c)

### メインループ開始
(1..nt).each{|it|
  # du / dx (2 次の中心差分)
  dudx = (u.sft(1) - ??????) / ??????
  dudt = ??????
  # 時間発展を解く (一次の前進差分)
  u = u + ??????

  # 可視化
  if it % nt_out == 0
    t = it * dt
    draw(u, x, t, u0, l, m, dx, c)
  end
### メインループ終了
}

### ウィンドウを閉じる
DCL.grcls