function stream505  

v0 = 1;
u0 = 1;
R0 = 1;
t0 = -5;
tf = 5;
dt = .02;
A = ones( (tf-t0)/dt+1,2);
B = A;
C = A;
i = 1;
j = 0;

[x,y] = meshgrid(-17:.33:8,-3:.33:5);

xA = u0*t0;
yA = -R0/2;
xB = u0*t0;
yB = -R0;
xC = u0*t0;
yC = -3*R0/2;


%Initial frames


%1. Rankine vortex profile

[o,r] = cart2pol(x,y);
uoin = (r <= R0);
uout = (r > R0);
uo = uoin+uout;
uo(uoin) =  v0*r(uoin)/R0;
uo(uout) =  v0*R0./r(uout);

u = -uo.*sin(o);
v = uo.*cos(o);	

hold off;
axis([-7,7,-2,2]);
q = axis;
quiver(x,y,u,v);
hold on;
c = circle([0,0],R0,100,'-');
set(c, 'Color', 'black')
axis equal
axis(q);
hold off
pause(3);
imwrite(frame2im(getframe),['startimg1.png']);
%print('-dpng', '-r0', 'startimg1.png]);

%2. Vort+Mean


    [o,r] = cart2pol(x,y);

    ur = u0*cos(o);

    uoin = (r <= R0);
    uout = (r > R0);
    uo = uoin+uout;
    uo(uoin) =  v0*r(uoin)/R0 - u0*sin(o(uoin));
    uo(uout) =  v0*R0./r(uout) - u0*sin(o(uout)) ;

    u = ur.*cos(o) - uo.*sin(o);
    v = ur.*sin(o) + uo.*cos(o);


hold off;
axis([-7,7,-2,2]);
q = axis;
quiver(x,y,u,v);
hold on;
c = circle([0,0],R0,100,'-');
set(c, 'Color', 'black')
axis equal
axis(q);
hold off
pause(3);

imwrite(frame2im(getframe),['startimg2.png']);

%3. Init quiver + streamlines




for t = t0:dt:tf

    [o,r] = cart2pol(x-u0*t,y);

    ur = u0*cos(o);

    uoin = (r <= R0);
    uout = (r > R0);
    uo = uoin+uout;
    uo(uoin) =  v0*r(uoin)/R0 - u0*sin(o(uoin));
    uo(uout) =  v0*R0./r(uout) - u0*sin(o(uout)) ;

    u = ur.*cos(o) - uo.*sin(o);
    v = ur.*sin(o) + uo.*cos(o);

%    trajectory plotting
     
%    need to find velocities at each point;
%    using 2nd order Runge-Kutta scheme

     [oA, rA] = cart2pol(xA-u0*t,yA);
     uA =  -v0*rA/R0*sin(oA)+u0;
     vA =  v0*rA/R0*cos(oA);
     xnp1 = xA+uA*dt;
     ynp1 = yA+vA*dt;

     [oAA, rAA] = cart2pol(xnp1-u0*t,ynp1);
     uAA =  -v0*rAA/R0*sin(oAA)+u0;
     vAA =  v0*rAA/R0*cos(oAA);     
     xA = xA+dt/2*(uA+uAA);
     yA = yA+dt/2*(vA+vAA);

     A(i,:) = [xA, yA];

     [oB, rB] = cart2pol(xB-u0*t,yB);
     uB =  -v0*R0/rB*sin(oB)+u0;
     vB =  v0*R0/rB*cos(oB);
     xnp1 = xB+uB*dt;
     ynp1 = yB+vB*dt;

     [oBB, rBB] = cart2pol(xnp1-u0*t,ynp1);
     uBB =  -v0*rBB/R0*sin(oBB)+u0;
     vBB =  v0*rBB/R0*cos(oBB);     
     xB = xB+dt/2*(uB+uBB);
     yB = yB+dt/2*(vB+vBB);

     B(i,:) = [xB, yB];

     [oC, rC] = cart2pol(xC-u0*t,yC);
     uC =  -v0*R0/rC*sin(oC)+u0;
     vC =  v0*R0/rC*cos(oC);
     xnp1 = xC+uC*dt;
     ynp1 = yC+vC*dt;

     [oCC, rCC] = cart2pol(xnp1-u0*t,ynp1);
     uCC =  -v0*rCC/R0*sin(oCC)+u0;
     vCC =  v0*rCC/R0*cos(oCC);     
     xC = xC+dt/2*(uC+uCC);
     yC = yC+dt/2*(vC+vCC);

     C(i,:) = [xC, yC];

     if rem(i-1,10) == 0

	axis([-7,7,-2,2]);
	q = axis;
	quiver(0,0,0,0);
	hold on;
	c = circle([u0*t,0],R0,100,'-');
	set(c, 'Color', 'black')
	h = streamline(x,y,u,v,(t0-tf-2)*ones(1,10)+t,-2:.70:4.8);
	set(h, 'Color', 'red');

	if i == 1

	   quiver(x,y,u,v);
	   axis equal;
	   axis(q);
	   pause(3);
	   imwrite(frame2im(getframe),['startimg3.png']);

	end

	plot(xA,yA,'.k'); 
	plot(xB,yB,'.c'); 
	plot(xC,yC,'.g'); 
	axis equal;
	axis(q);

	plot(A(1:i,1),A(1:i,2), '-b', 'LineWidth',2)
	plot(B(1:i,1),B(1:i,2), '-c', 'LineWidth',2)
	plot(C(1:i,1),C(1:i,2), '-g', 'LineWidth',2)

	
	arrow([xA(end),yA(end)],[xA(end)+uA/4,yA(end)+vA/4], ...
	      'Length',10,'Width',2);
	arrow([xB(end),yB(end)],[xB(end)+uB/4,yB(end)+vB/4], ...
	      'Length',10,'Width',2);
	arrow([xC(end),yC(end)],[xC(end)+uC/4,yC(end)+vC/4], ...
	      'Length',10,'Width',2);

	hold off;
	getframe;

	j = j+1;
	imwrite(frame2im(getframe),['image', int2str(j) , '.png']);

     end

    %getframe;
    %m(i) = getframe;


%    hold off;

    i = i+1;
end

%movie(m,(tf-t0/dt+1))
%clear;
