%DEMO1  Demo of control volume deformation in a 2D streamflow
%   Simulates the evolution of a blob of fluid in a streamflow.
%   Uses the file demo_ufunc.m for its odefile.
%
%   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: jrg, cbm (08/02/2002)
%

% ****

global DOMAIN FLOW_TITLE;

% Customizable parameters

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 blob a total of N+1 times
N = 6;

% take skipn steps between each drawing of the blob
% (each step is .1)
skipn = 3;


% the initial size of the box
boxsize = .20;


% total number of points in blob (*must* be divisible by 4)
numpts = 20*4;

% ***
% 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(1);
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 sphere
  th = 0:pi/50:2*pi;
  plot(cos(th),sin(th),'k')
end




% graphical IVs
disp('  clik for IV');  disp(' ')
IVinput = ginput(1);
% generate complex initial conditions for midpoint of the blob
midpt = IVinput(1) + i * IVinput(2);

% generate complex initial conditions for boundaries of the blob
edgepts = zeros(numpts,1);
S = numpts/4;
w = boxsize;
h = boxsize;
m1 = real(midpt);
m2 = imag(midpt);
for j=0:(S-1)
  edgepts(j+1)     = (m1-w/2+j*w/S) + i * (m2+h/2);
  edgepts(j+S+1)   = (m1+w/2)       + i * (m2+h/2-j*h/S);
  edgepts(j+2*S+1) = (m1+w/2-j*w/S) + i * (m2-h/2);
  edgepts(j+3*S+1) = (m1-w/2)       + i * (m2-h/2+j*h/S);
end



plot(midpt, 'k*');
% note we need to plot the first point twice to close the shape
plot([edgepts; edgepts(1)], 'k-');

% force a redraw
pause(0);

tstep = .1;
t0 = 0;
for n = 1:(N*skipn)
  [T,MIDPT] = ode45(@demo_ufunc, [t0 (t0+tstep)], midpt, [], flo);
  [T,EDGEPTS] = ode45(@demo_ufunc, [t0 (t0+tstep)], edgepts, [], flo);

  % trace the streamline
  plot(MIDPT, 'r-');
    
  % the output from ode45 contains more then just the "answer"
  midpt = MIDPT(end);
  edgepts = transpose(EDGEPTS(end,:));
  

  if (mod(n,skipn) == 0)
    plot(midpt, 'ko');
    plot([edgepts; edgepts(1)], 'k-');
  end

  % force a redraw
  pause(0);  

  % get ready for the next step
  t0 = t0 + tstep;
end

title(['\bf control volume deformation:  ' FLOW_TITLE])
xlabel('\bf x-axis')
ylabel('\bf y-axis')