clear all;

%smethod = 'AB1';
%smethod = 'AB2';
%smethod = 'AM1';
%smethod = 'AM2';
%smethod = 'RK2';
smethod = 'RK4';

dReta = 0.0025;
Reta  = 0.0;
rpmax = 0.0;

while ((rpmax < 1.0) & (Reta < 1.0)),
  
  % Initialize parameters

  part_init;
  Reta = Reta + dReta;
  R0 = Reta*Ro;
  
  
  % Initial conditions
  
  u0 = 0; 
  v0 = 0;
  w0 = Wg;
  x0 = R0;
  y0 = 0;
  z0 = 0;
  
  % Set time length of integration, and number of steps
  
  Tmax  = L/Wg;
  N     = 10;
  dt    = Tmax/N;
  
  % Integrate particle position

  if (smethod == 'AB1'),
    part_ab1
  elseif (smethod == 'AB2'),
    part_ab2;
  elseif (smethod == 'AM1'),
    part_am1;
  elseif (smethod == 'AM2'),
    part_am2;
  elseif (smethod == 'RK2'),
    part_rk2;
  elseif (smethod == 'RK4'),
    part_rk4;
  else,
    fprintf('ERROR: unknown smethod value = %s\n',smethod);
    exit;
  end
  
  % Post-process to find rpmax = max(rp)/Rs
  xp = v(4,:);
  yp = v(5,:);
  rp = sqrt(xp.^2 + yp.^2);
  
  rpmax = max(rp)/Rs;
  
end

fprintf('%s: dt = %4.2e Reta = %5.3f  eta = %5.3f\n',smethod,dt,Reta,1-Reta^2);