%---------------------------------------------------------  
% 複数の質点の運動を計算するプログラム
%---------------------------------------------------------
function vortex
  
  global g;

%------------------------------
  h=0.1; % 画面表示する t の間隔
  nstep=1000; % 時間発展するステップ数
              % ( t = nstep * h まで時間発展されることになる)
  m=10; % h をルンゲクッタ法においてどれだけ分割して計算するか
       % (m が大きい程精度は上がるが, m に比例して計算時間がかかる)
%------------------------------  
  q(1,1)=0;    % 1番目の質点の x座標の初期値
  q(2,1)=0;    % 1番目の質点の y座標の初期値
  q(3,1)=0;    % 1番目の質点の z座標の初期値
  q(4,1)=0;    % 1番目の質点の x方向速度(u)の初期値
  q(5,1)=0;    % 1番目の質点の y方向速度(v)の初期値
  q(6,1)=0;    % 1番目の質点の z方向速度(w)の初期値
  g(1)=1;      % 1番目の質点の質量
  
  q(1,2)=3;    % 2番目の質点の x座標の初期値
  q(2,2)=0;    % 2番目の質点の y座標の初期値
  q(3,2)=0;    % 2番目の質点の z座標の初期値  
  q(4,2)=0.3;  % 2番目の質点の x方向速度(u)の初期値
  q(5,2)=0.2;  % 2番目の質点の y方向速度(v)の初期値
  q(6,2)=0;    % 2番目の質点の z方向速度(w)の初期値
  g(2)=0.001;  % 2番目の質点の質量
  
  % 質点の位置や質量の値を変えたければ, 上記の値を変えればよい.
  % 3番目, 4番目の渦を追加したければ, 続けて 
  % q(1,3)=.... などとして増していけばよい.
  
%------------------------------
% グラフィクスのための設定  
  gset parametric
  gset nokey
  axis([-1,4,-1,1.5,0,1],"equal") 
    % 画面の表示範囲の設定(x,y,z の最小, 最大を順に並べる)
  gset view 60,30
  % 視点の設定(1番目の数値が x軸まわり, 2番目の数値が z軸まわりの回転角)
  % gset view 0,0 % なら xy平面への投影を表示することになる.
  % hold on % もし軌跡を描きたいときは, このコメントを外せばよい.  
%------------------------------  

  t=0;
  graph3d(t,q); % 質点の位置の表示(初期値)
  
  for istep=1:nstep
    [t,q]=rk4u(m,h,t,q,'gfunc');
    graph3d(t,q); % 質点の位置の表示
    sleep(0.1); % 画面表示を何秒毎に更新するかの設定
  end