% LVDEMO  demo of multiple line vorticies
%
%  2 vortices:  +1 at (-0.75,+0.4)  and
%               -1 at (-0.75,-0.4)
%
%  4 vortices:  +1 at (-0.75,+0.4),(0.0,+0.4)  and  
%               -1 at (-0.75,-0.4),(0.0,-0.4)
%
% cbm, 13/03/2002

global VORTEX_RAD_SQ;
% don't ask -- look at lvdemo_vel4quiv.m if you really must know
VORTEX_RAD_SQ = sqrt(eps);

clear;

% how dense the quiver field is
VF_DENSITY = 30;
% the size of the domain
DOMAIN.left = -1; DOMAIN.right = 1;
DOMAIN.bottom = -1; DOMAIN.top = 1;

% the colors (RGB triplets) to cycle through for plotting the LVs
colors = [1 0 0; 0 1 0; 0 0 1; 0 0 0; 1 0 1; 0 1 1];
ncolors = length(colors);

% start and end time
t0 = 0; tf = 5;
% small timesteps
tstep = .01;
% how often should we stop to replot the quivers
quiver_skip_steps = 5;

% end of customizable parameters

M = round(tf/tstep);
N = input('enter number of line vortices: ');
lvs = zeros(M, N);

for j = 1:N
  circs(j) = input(['enter the circulation for vortex ' num2str(j) ': ']);
end

figure(1);
clf;
axis equal;
axis([DOMAIN.left DOMAIN.right DOMAIN.bottom DOMAIN.top]);
hold on;

for j = 1:N
  disp(['click to place (*) vortex ' num2str(j) '...']);
  tempin = ginput(1);
  lvs(1,j) = tempin(1) + i*tempin(2);
  H = plot(lvs(1,j), '*');
  set(H, 'Color', colors(mod(j-1,ncolors)+1,:));  
end

t = t0;
for n = 1:M;
  if (mod(n-1, quiver_skip_steps) == 0)
    figure(1);
    clf;
    hold on;
    axis equal;
    axis([DOMAIN.left DOMAIN.right DOMAIN.bottom DOMAIN.top]);

    % draw the quiver field
    dx = (DOMAIN.right - DOMAIN.left) / VF_DENSITY;
    dy = (DOMAIN.top - DOMAIN.bottom) / VF_DENSITY;
    xx = (DOMAIN.left   + dx/2):dx:(DOMAIN.right - dx/2);
    yy = (DOMAIN.bottom + dy/2):dy:(DOMAIN.top   - dy/2);
    [xg, yg] = meshgrid(xx, yy);
    uvg = lvdemo_vel4quiv(0, xg + i*yg, lvs(n,:), circs);

    ug = real(uvg);
    vg = imag(uvg);
    % divide the quivers by their own lengths
    lens = sqrt(ug.^2 + vg.^2);
    ug = ug./lens;
    vg = vg./lens;

    H = quiver(xx, yy, ug, vg, .3);
    set(H, 'Color', [.6 .6 .6]);
    H = quiver(xx, yy, -ug, -vg, .3, '.');
    set(H, 'Color', [.6 .6 .6]);

    for j = 1:N
      H = plot(lvs(1,j), '*');
      set(H, 'Color', colors(mod(j-1,ncolors)+1,:));  
      H = plot(lvs(1:n,j));
      set(H, 'Color', colors(mod(j-1,ncolors)+1,:));
      H = plot(lvs(n,j), 'o');
      set(H, 'Color', colors(mod(j-1,ncolors)+1,:));
    end
    
    disp('Press enter to proceed');
    pause
  end
  
  
  [T,LVS] = ode45(@lvdemo_ode, [t t+tstep], lvs(n,:), [], circs);

  lvs(n+1,:) = LVS(end,:);
  for j = 1:N
    H = plot([lvs(n,j) lvs(n+1,j)]);
    set(H, 'Color', colors(mod(j-1,ncolors)+1,:));
  end
  pause(0);
  t = t + tstep;
end