#!/usr/bin/env ruby
# -*- coding: utf-8 -*-
#
#*******************************************************************************
# 表題 地球流体力学 計算機実験集
#
#     粒子拡散：
#     ブラウン運動による多数粒子群の拡散
#
# 履歴 1999/04/24 金尾 学
#      1999/05/11 余田 成男
#      2012/03/05 村上 真也 (Ruby/DCL版)
#*******************************************************************************

NMAX=3000  # 粒子数
ITMAX=200  # 繰り返し計算回数
DT=0.1      # 時間刻み
JBOX=40    # ヒストグラムのボックス数
XMIN, XMAX, YMIN, YMAX = -5.0, 5.0, -5.0, 5.0  # ウインドウの範囲

####
require 'narray'
require 'numru/dcl'
include NumRu

#*******************************************************************************
def random()
  
  # 正規乱数作成サブルーチン
  #     水田 亮 作
  #     参考文献: Numerical Recipes
  
  begin
    v1 = 2.0*rand() - 1.0
    v2 = 2.0*rand() - 1.0
    r  = v1**2 + v2**2
  end until r <= 1.0
  
  factor = Math::sqrt( -2.0*Math::log(r)/r )
  
  return v1*factor, v2*factor
end
#*******************************************************************************
def initia(x, y, nmax, iseed)

  # 各粒子の初期位置は, 乱数を用いて半径rの円内に一様等方的に配置
  r = 0.2
  (0..nmax-1).each do |n|
    begin
      px = r*(2.0*rand() - 1.0)
      py = r*(2.0*rand() - 1.0)
      dd = px**2 + py**2
    end until dd < r**2
    x[0,n] = px
    y[0,n] = py
  end
end
#*******************************************************************************
def timein(x, y, a0, a1, b, iseed, nmax, itmax, dt)

  # 各粒子毎に時間発展
  (0..nmax-1).each do |n|
    # 各粒子の初期速度はゼロ
    u = 0.0
    v = 0.0
    (1..itmax).each do |i|
      srand1, srand2 = random()
      du =      - a1*u + b*srand1
      dv = - a0 - a1*v + b*srand2
      u  = u + du*dt
      v  = v + dv*dt
      x[i,n] = x[i-1,n] + du*dt
      y[i,n] = y[i-1,n] + dv*dt
    end
  end
end
#*******************************************************************************
def xypdf(x, y, pdfx, pdfy, istep, nmax, itmax, jbox, xmin, xmax, ymin, ymax)

  # 分布関数の計算
  delx = (xmax-xmin)/jbox
  dely = (ymax-ymin)/jbox
  
  DCL::rset(pdfx, jbox*(itmax+1), 1, 0.0)
  DCL::rset(pdfy, jbox*(itmax+1), 1, 0.0)
  
  (0..itmax/istep).each do |ii|
    i = ii*istep
    (0..nmax-1).each do |n|
      j=(0..jbox-1).each do |j|
        xl = xmin + delx*j
        xm = xl    + delx
        if (x[i,n] >= xl and x[i,n] <= xm)
          pdfx[j,i] = pdfx[j,i] + 100.0/nmax # %で表示
        end
          
        yl = ymin + dely*j
        ym = yl    + dely
        if (y[i,n] >= yl and y[i,n] <= ym)
          pdfy[j,i] = pdfy[j,i] + 100.0/nmax # %で表示
        end
      end
    end
  end
end
#*******************************************************************************
def graph(x, y, pdfx, pdfy, istep, dt, nmax, itmax, xmin, xmax, ymin, ymax)
  
  # 絵を描く
  
  pmax=12.0
  
  rundef = DCL::glrget('RUNDEF')
  DCL::swiset('iwidth',  540)
  DCL::swiset('iheight', 540)
  DCL::swlset('lalt',   true)
  DCL::swlset('lwait',  false)
  DCL::sglset('lcorner',false)
  DCL::sglset('lclip',  true)
  
  DCL::gropn(4)
  
  # istep 毎に絵を描く
  (0..itmax/istep).each do |ii|
    i = ii*istep
    DCL::grfrm
    DCL::grswnd(xmin, xmax, ymin, ymax)
    DCL::grsvpt(0.15, 0.75, 0.15, 0.75)
    DCL::grstrn(1)
    DCL::grstrf
    
    DCL::uxaxdv('b', 0.5, 2.0)
    DCL::uxaxdv('t', 0.5, 2.0)
    DCL::uyaxdv('l', 0.5, 2.0)
    DCL::uyaxdv('r', 0.5, 2.0)
    DCL::uxmttl('b', 'X', 0.0)
    DCL::uymttl('l', 'Y', 0.0)
    
    # 粒子群
    DCL::sgspmi(41)
    DCL::sgpmu(x[i,true], y[i,true])
    
    # x-分布関数
    DCL::grfig
    DCL::grswnd(XMIN, XMAX, 0.0, pmax)
    DCL::grsvpt(0.15, 0.75, 0.77, 0.95)
    DCL::grstrn(1)
    DCL::grstrf
    
    DCL::uuslni(22)
    DCL::uusidv(XMIN, XMAX)
    DCL::uvbxl(rundef, pdfx[true,i])
    
    DCL::uzlset('labelxb', false)
    DCL::uxaxdv('b', 0.5, 2.0)
    DCL::uzlset('labelxb', true)
    DCL::uxaxdv('t', 0.5, 2.0)
    DCL::uyaxdv('l', 1.0, 5.0)
    DCL::uyaxdv('r', 1.0, 5.0)
    DCL::uysttl('l', 'PDF(%)', 0.0)
    
    # y-分布関数
    DCL::grfig
    DCL::grswnd(0.0, pmax, YMIN, YMAX)
    DCL::grsvpt(0.77, 0.95, 0.15, 0.75)
    DCL::grstrn(1)
    DCL::grstrf
    
    DCL::uusidv(YMIN, YMAX)
    DCL::uhbxl(pdfy[true,i], rundef)
    
    DCL::uxaxdv('b', 1.0, 5.0)
    DCL::uxaxdv('t', 1.0, 5.0)
    DCL::uzlset('labelyl', false)
    DCL::uyaxdv('l', 0.5, 2.0)
    DCL::uzlset('labelyl', true)
    DCL::uyaxdv('r', 0.5, 2.0)
    DCL::uxsttl('b', 'PDF(%)', 0.0)
    
    # 時刻
    DCL::grfig
    DCL::grswnd(0.0 , 1.0, 0.0 , 1.0)
    DCL::grsvpt(0.75, 1.0, 0.75, 1.0)
    DCL::grstrn(1)
    DCL::grstrf
    
    DCL::sgstxi(2)
    DCL::sgstxs(0.03)
    DCL::sgtxv(0.88, 0.88, sprintf("T=%4.1f", i*dt))
  end
  
  DCL::grcls
end
#*******************************************************************************

x    = NArray.float(ITMAX+1,NMAX)
y    = NArray.float(ITMAX+1,NMAX)
pdfx = NArray.float(JBOX,ITMAX+1)
pdfy = NArray.float(JBOX,ITMAX+1)

# 実験パラメータの読み込み
file = 'expar2.data'
if RUBY_VERSION >= '1.9.0'
  io = open(file, 'r+', {:external_encoding=>"euc-jp"} )
else
  io = open(file, 'r+')
end
l = io.readline.split(' ') # 1行目のみ読みだす
a0    = l[0].to_f   # 定常外力の大きさ
a1    = l[1].to_f   # 線型外力の大きさ(レーリー摩擦の係数)
b     = l[2].to_f   # ストカスティック外力の大きさ
iseed = l[3].to_i   # 疑似乱数用の種
istep = l[4].to_i   # グラフ出力のステップ間隔
io.close

# 疑似乱数初期化
srand(iseed)

# 初期値設定
initia(x, y, NMAX, iseed)
# 時間発展
timein(x, y, a0, a1, b, iseed, NMAX, ITMAX, DT)
# 分布関数
xypdf(x, y, pdfx, pdfy, istep, NMAX, ITMAX, JBOX, XMIN, XMAX, YMIN, YMAX)
# グラフ
graph(x, y, pdfx, pdfy, istep, DT, NMAX, ITMAX, XMIN, XMAX, YMIN, YMAX)