% % code02.m % % save file as "code02.m" (make sure that it is on the matlab path!) % to run, type "code02" at the matlab prompt ">>" % % initialize workspace & set constants(try B=0,1,2,3) clear; close all L = 4.0; ds = 2*L/45; B = 1.0 % set x & y-vectors & grid matrices ("help meshgrid") x = -L:ds:+L; y = -L:ds:+L; [xx,yy] = meshgrid(x,y); rr2 = xx.^2 + yy.^2; % calculate psi psi = yy .* (1 - 1./(rr2)) + (B/2)*log(rr2); % calculate velocities (u = psi_y, v = -psi_x) uu = (1 - 1./(rr2)) + (2*yy.*yy ./ (rr2.^2) ) + (B*yy ./ rr2 ); vv = - (2*xx.*yy ./ (rr2.^2) ) - (B*xx ./ rr2 ); % plot vector field ("help quiver") outside unit circle ("help gt") clf quiver(x,y, uu .* (rr2>=1), vv .* (rr2>=1),0.3) hold on quiver(x,y,-uu .* (rr2>=1),-vv .* (rr2>=1),0.3,'.') axis([-L +L -L +L]); axis square; % plot unit circle th = 0:pi/50:2*pi; plot(cos(th),sin(th),'k') % plot level curves of psi (huh?) contour(x,y,psi .* (rr2>=1),10) % label plot title(['\bf my first vector plot (B=' num2str(B) '\bf)']) xlabel('\bf x-axis') ylabel('\bf y-axis') % plot a point xs = 0; ys = 0; plot(xs,ys,'r*')