%DEMO2  Demo of a "vorticity meter" in a 2D streamflow
%   Demonstrates the difference between a rotational flow
%   and a irrotational one using a "vorticity meter".

%   bugs: if your initial conditions have no complex component
%   then strange things can happen...

% orig by: cbm (24/01/2002)
% hacked on by: cbm, djm
% 'complex'ified: cbm (10/02/2002)

% ****

global DOMAIN FLOW_TITLE

flo = [];  
while (isempty(flo)==1)
  disp(' ')
  disp('  flow choices:')
  disp('    1) line vortex')
  disp('    2) solid body rotation')
  disp('    3) shear')
  disp('    4) flow past cylinder')
  disp('    5) user-specified problem')
  disp(' ')

  flo = input(' choice:  ');
end


% plot the vorticity meter a total of N times
N = 7;

% take skipn steps between each drawing of the vorticity meter
% (each step is .01)
skipn = 10;


% the true initial size of the VM
vmradius = .01;
% how many times should we scale up the VM for visualization?
scaling = 10;


% ****
% Changes past this point should not typically be necessary



% define global parameters by making an extra call to demo_ufunc
demo_ufunc(0,1,flo);


% initial figure setup
figure(2);
clf;
axis([DOMAIN.left DOMAIN.right DOMAIN.bottom DOMAIN.top]);
axis square;
hold on;


% plot vector field
VF_DENSITY = 10;
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 = demo_ufunc(0, xg + i*yg, flo);
ug = real(uvg);
vg = imag(uvg);

quiver(xx, yy, ug, vg, .3)
quiver(xx, yy, -ug, -vg, .3, '.')

if (flo == 4)
  % draw a unit circle
  th = 0:pi/50:2*pi;
  plot(cos(th),sin(th),'k')
end



% graphical IVs
disp('  clik for IV');  disp(' ')
IVinput = ginput(1);

old_midpt = complex(IVinput(1),IVinput(2));

% generate the corner points of the VM
theta = 0;
offsets = [0; pi/2; pi; 3*pi/2];
old_corners = old_midpt + vmradius*exp(i*offsets);

% plot the initial VM
scal_corners = old_midpt + scaling*vmradius*exp(i*(theta+offsets));
plot(scal_corners(1), 'k*');
plot([scal_corners(1), scal_corners(3)], 'k-');
plot([scal_corners(2), scal_corners(4)], 'k-');
pause(0);


tstep = .01;
t0 = 0;
for n = 1:(N*skipn)
  [T,MIDPT] = ode45(@demo_ufunc, [t0 (t0+tstep)], old_midpt, [], flo);
  [T,CORNERS] = ode45(@demo_ufunc, [t0 (t0+tstep)], old_corners, [], flo);

  % trace the streamline
  plot(MIDPT, 'r-');

  % the output from ode45 contains more then just the "answer"
  new_midpt = MIDPT(end);
  new_corners = transpose(CORNERS(end,:));

  % the average change in the 4 angles formed by the corners and
  % the midpoint
  dtheta = mean(angle((new_corners-new_midpt)./(old_corners-old_midpt)));
  theta = theta + dtheta;
    
  if (mod(n,skipn) == 0)
    scal_corners = new_midpt + scaling*vmradius*exp(i*(theta+offsets));
    plot(scal_corners(1), 'ko');
    plot([scal_corners(1), scal_corners(3)], 'k-');
    plot([scal_corners(2), scal_corners(4)], 'k-');
  end
  
  % force a redraw
  pause(0);  
 
  % prepare for the next iteration (basically reset the VM to its
  % original position parallel to the x,y axes)
  old_corners = new_midpt + vmradius*exp(i*offsets);
  old_midpt = new_midpt;
  t0 = t0 + tstep;
end

title(['\bf vorticity meter:  ' FLOW_TITLE]);
xlabel('\bf x-axis');
ylabel('\bf y-axis');