
disp('   Another way of interpreting divergence  ')
disp('   We use the vector field     ')
disp('    F(x,y) = (u(x,y), v(x,y))  =  (x+y, cos(y))  ')
disp('    div(F) = 1 - sin(y)  ')

 xmin = -3; xmax = 1;
 width = xmax - xmin;
 ymin = -1.5; ymax = ymin+.8*width; 

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

 U = X +Y; V = cos(Y);

 quiver(X,Y,U,V)
 axis([xmin,xmax, ymin, ymax])
 hold on
 
 times = [.4, .8, 1.2, 1.6];

 %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')

disp(' Hit return to see where the flow carries the dot ')

pause

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

f = inline('[w(1)+w(2); cos(w(2))]', 't', 'w');

for j = 1:n+1

   w0 = [x0(j), y0(j)];
   [t,w] = ode45(f, 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

pause
 
close

