%
%        Function mfile  flow2.m
%
%   Function mfile to compute the deformation of a disk following
%   the flow of a vector field.  The call is
%   flow2(corners, times) where as usual corners = [a,b,c,d]
%   is the vector of corner coordinates of the rectangle where the
%   vector field if displayed. times = [t1 t2 t3 t4] is a vector of
%   four times. The deformed disk is displayed at the times
%   t1, t2, t3,t4.    The vector field must be provided in
%   mfiles u.m and v.m.   The mfile wdot.m is also needed. 
%   t1,t2,t3,t4 are the times at which the deformed disk is 
%   observed. After the call, the program asks user to click on the
%   to place a disk. The image of the disk under the flow is then
%   calculated and plotted using ode45 on the boundary points of the
%   disk. The area of each image is calculated approximately.



function out = flow2(corners, times)
 xmin = corners(1); xmax = corners(2);
 width = xmax - xmin;
 ymin = corners(3); ymax = ymin+.8*width; 

 x = linspace(xmin,xmax,16); y = linspace(ymin,ymax ,16);
 [X,Y] = meshgrid(x,y);

 U = u(X,Y); V = v(X,Y);

quiver(X,Y,U,V)
axis([xmin,xmax, ymin, ymax])
hold on

%center = input('enter coords of center of circle [a,b]  ')
%a = center(1); b = center(2);

disp(' click where you want to the flow to begin   ')
disp('  ')
[a,b] = ginput(1)
r = width/20
n = 20; theta = 0:2*pi/n: 2*pi;

x0 = a + r*cos(theta); y0 = b + r*sin(theta) ;
fill(x0, y0, 'r')

x1 = zeros(1,n); y1 = x1;
x2 = x1; y2 = y1;
x3 = x2; y3 = y2;
x4 = x3; y4 = y3;
tspan = [0 times];

for j = 1:n+1

   w0 = [x0(j), y0(j)];
   [t,w] = ode45('wdot', tspan, w0);
   x1(j) = w(2,1); y1(j) = w(2,2);
   x2(j) = w(3,1); y2(j) = w(3,2);
   x3(j) = w(4,1); y3(j) = w(4,2);
   x4(j) = w(5,1); y4(j) = w(5,2);
end

area0 = 0;
for j = 1:n
    area0 = area0 +.5*(x0(j+1) +x0(j))*(y0(j+1) - y0(j));
end
area0


fill(x1,y1, 'r')
area1 = 0;
for j = 1:n
    area1 = area1 + .5*(x1(j+1) + x1(j))*(y1(j+1) - y1(j));
end
area1
    
fill(x2,y2, 'r')
area2 =0;
for j = 1:n
    area2 = area2 + .5*(x2(j+1) +x2(j))*(y2(j+1) - y2(j));
end
area2

fill(x3,y3, 'r')
area3 = 0;
for j = 1:n
    area3 = area3 + .5*(x3(j+1) + x3(j))*(y3(j+1) - y3(j));
end
area3

fill(x4,y4, 'r')
area4 = 0;
for j = 1:n
  area4 = area4 + .5*(x4(j+1) + x4(j))*(y4(j+1) - y4(j));
end
area4

hold off