

clear ; close;
disp(' ')
disp('                    CURL OF VECTOR FIELD    ')

disp(' ')
disp('      Four different vector fields are presented  ')
disp('  The scalar quantity  v_x - u_y is indicated by the colorbar on ')
disp('  the right. You can then click on the figure window 4 times  ')
disp('  to see the the direction of rotation ')

disp(' ')
disp('                        Shear Flow  ')
disp(' ')
disp('            F = (u(x,y), v(x,y))  = ( 1-y^2, 0)  ')
    x = linspace(-1.2,1.2, 51);
    [X,Y] = meshgrid(x,x);
    crl = inline('2*y', 'x', 'y');
    pcolor(X,Y, crl(X,Y)); shading interp; colorbar;
    
    hold on
    xx = linspace(-1,1, 11);
    [XX,YY] = meshgrid(xx,xx);
    quiver(XX,YY,1-YY.^2, 0*YY)
    
    
    for j =1:4
    drawnow
    [p1,p2] = ginput(1);
 
 s = linspace(0,2*pi, 101);
 plot(p1 + .15*cos(s), p2+.15 *sin(s), 'k')
 if crl(p1, p2) > 0
     theta = .5*pi + .2;
     tip = [p1 + .15*cos(theta), p2+ .15*sin(theta)];
     plot([.07, 0, .07] +tip(1), [.04, 0, -.04]+tip(2), 'k')
 elseif crl(p1,p2) < 0     
     theta = .5*pi -.2;
     tip = [p1+.15*cos(theta), p2+.15*sin(theta)]; 
     plot([-.07 0 -.07]+ tip(1), [.04, 0, -.04] + tip(2),'k')

 end
end 
    hold off
   pause
   clf
    

disp(' ')
disp('  Counter clockwise rotational flow with positive curl  ')
disp('   F(x,y) = ( -y/(r^2+1), x/(r^2+1)),   r = sqrt(x^2 +y^2)  ')
x = linspace(-2,2, 51);
[X,Y] = meshgrid(x,x);
clear crl
crl = inline('2./(x.^2 + y.^2+1).^2', 'x', 'y');
  pcolor(X,Y,crl(X,Y)); shading interp; 
  map = jet; myjet = map(33:64, :); colormap(myjet);
  colorbar;
  hold on
  
  xx = linspace(-1.7, 1.7 ,11);
  [XX,YY] = meshgrid(xx,xx);
  U = -YY./(XX.^2+YY.^2 +1);
  V = XX./(XX.^2 +YY.^2+1);
  quiver(XX,YY,U,V)
  
  for j =1:4
 drawnow
 [p1,p2] = ginput(1);
 
 s = linspace(0,2*pi, 101);
 plot(p1 + .15*cos(s), p2+.15 *sin(s), 'k')
 if crl(p1, p2) > 0
     theta = .5*pi + .2;
     tip = [p1 + .15*cos(theta), p2+ .15*sin(theta)];
     plot([.07, 0, .07] +tip(1), [.04, 0, -.04]+tip(2), 'k')
 elseif crl(p1,p2) < 0     
     theta = .5*pi -.2;
     tip = [p1+.15*cos(theta), p2+.15*sin(theta)]; 
     plot([-.07 0 -.07]+ tip(1), [.04, 0, -.04] + tip(2),'k')

 end
end 
  
  
  
  hold off
  
  
  pause
  clf
 
  
 disp(' ')
 disp(' Counter clockwise rotational Flow with positive and negative curl ')
 disp('        F(x,y) = (-yg(r), xg(r) ),  g(r) = (rsin(r) + cos(r) -1)/r^2  ')
 
 clear Z U V 
 R  = sqrt(X.^2 +Y.^2 + eps);
 crl = inline('cos(sqrt(x.^2 + y.^2))', 'x', 'y');;
 
 
 pcolor(X,Y,crl(X,Y)); shading interp; colorbar;colormap(jet)
 hold on
 contour(X,Y,crl(X,Y), [0,0],'k') 
 
 RR= sqrt(XX.^2+YY.^2+eps);
 G = (RR.*sin(RR) + cos(RR) -1)./RR.^2;
 
 
 
 U = -YY.*G;
 V = XX.*G;
 quiver(XX,YY,U,V)
 
 for j =1:4
 drawnow
 [p1,p2] = ginput(1);
 
 s = linspace(0,2*pi, 101);
 plot(p1 + .15*cos(s), p2+.15 *sin(s), 'k')
 if crl(p1, p2) > 0
     theta = .5*pi + .2;
     tip = [p1 + .15*cos(theta), p2+ .15*sin(theta)];
     plot([.07, 0, .07] +tip(1), [.04, 0, -.04]+tip(2), 'k')
 elseif crl(p1,p2) < 0     
     theta = .5*pi -.2;
     tip = [p1+.15*cos(theta), p2+.15*sin(theta)]; 
     plot([-.07 0 -.07]+ tip(1), [.04, 0, -.04] + tip(2),'k')

 end
end 
hold off
 pause
 clf

 
 disp(' ')
 disp('  General Vector field  ')
 disp('   F = (u(x,y), v(x,y)) = .5*( y^3, x^2 +2x)  ')
 
 clear Z U V R
 crl = inline(' x+1 -1.5*y.^2', 'x', 'y')
 pcolor(X,Y,crl(X,Y)); shading interp; colorbar;
 hold on
 contour(X,Y,crl(X,Y), [0,0], 'k')
 
 U = .5*YY.^3;
 V = .5*XX.^2 + XX;
 
 quiver(XX,YY, U,V)
 for j =1:4
 drawnow
 [p1,p2] = ginput(1);
 
 s = linspace(0,2*pi, 101);
 plot(p1 + .15*cos(s), p2+.15 *sin(s), 'k')
 if crl(p1, p2) > 0
     theta = .5*pi + .2;
     tip = [p1 + .15*cos(theta), p2+ .15*sin(theta)];
     plot([.07, 0, .07] +tip(1), [.04, 0, -.04]+tip(2), 'k')
 elseif crl(p1,p2) < 0     
     theta = .5*pi -.2;
     tip = [p1+.15*cos(theta), p2+.15*sin(theta)]; 
     plot([-.07 0 -.07]+ tip(1), [.04, 0, -.04] + tip(2),'k')

 end
end 
 
 pause
 close


