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

%------------------------------
  h=0.1; % 画面表示する t の間隔
  nstep=1000; % 時間発展するステップ数
              % ( t = nstep * h まで時間発展されることになる)
  m=1; % h をルンゲクッタ法においてどれだけ分割して計算するか
       % (m が大きい程精度は上がるが, m に比例して計算時間がかかる)
%------------------------------  
  z(1,1)=0;  % 1番目の渦の x座標の初期値
  z(2,1)=1;  % 1番目の渦の y座標の初期値
  g(1)=1;    % 1番目の渦の循環Γの値
  
  z(1,2)=0;  % 2番目の渦の x座標の初期値
  z(2,2)=-1; % 2番目の渦の y座標の初期値
  g(2)=1;    % 2番目の渦の循環Γの値
  
  % 渦の位置や循環の値を変えたければ, 上記の値を変えればよい.
  % 3番目, 4番目の渦を追加したければ, 続けて 
  % z(1,3)=.... などとして増していけばよい.
  
  %------------------------------
  
%  hold on % もし軌跡を描きたいときは, このコメントを外せばよい.

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