close all; clear all; clc

%% define grid
V0     = 1;
ax     = [-3 3 0 2];        % plot axis (and control volume)
nx     = 51;                % no. of points in direction x
ny     = 21;                % no. of points in direction y

% define vector x and y
x = linspace(ax(1),ax(2),nx);
y = linspace(ax(3),ax(4),ny);

% create mesh
[xx,yy] = meshgrid(x,y);

% choose the departing position
x0 = 0.11;  % departing position in x
y0 = 1.8;   % departing position in y

% define time vector
dt=0.01; tf=4;
t=0:dt:tf;
nt=length(t);

n = input('Enter 1 for steady(L=1) and 2 for unsteady case: ');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

switch n
    %% A. Steady Flow, L=1
    case 1
        L      = 1;
        
        % 2: flow field
        v = -V0*yy./cosh(xx/L).^2;
        u = V0*L*tanh(xx/L);
        quiver(xx,yy,u,v,0.8,'k');
        
        % 3.Flux computation
        vtop   = -V0*y(end)./cosh(x/L).^2;
        uleft  = ones(size(y))*V0*L*tanh(x(1)/L);
        uright = ones(size(y))*V0*L*tanh(x(end)/L);
        
        fluxTop   = trapz(x,vtop);
        fluxLeft  = -trapz(y,uleft);
        fluxRight = trapz(y,uright);
        flucTot   = fluxTop+fluxLeft+fluxRight;
        
        % 5. Plotting streamlines, streaklines and trajectory
        
        vf = @(x,y,t) -V0*y/cosh(x/L)^2;
        uf = @(x,y,t)  V0*L*tanh(x/L);
        
        % xt,yt: trajectory points
        xt = zeros(nt,1);
        yt = zeros(nt,1);
        xt(1) = x0;
        yt(1) = y0;
        
        % xs,ys streamline points
        xs=zeros(nt,(nt-1)/10);
        ys=zeros(nt,(nt-1)/10);
        xs(1:10:end-1,:)=x0;
        ys(1:10:end-1,:)=y0;
        is=0;
        
        figure
        
        for ii=1:nt-1
            
            if mod(ii-1,10)==0 
                is=is+1;
                xs(ii,is)=x0;
                ys(ii,is)=y0;
                v = -V0*yy./cosh(xx/L).^2;
                u = V0*L*tanh(xx/L);
                
                hold off;
                quiver(xx,yy,u,v,0.8,'k');
                hold on;
                
                % plotting trajectory
                plot(xt(1:ii),yt(1:ii),'ro')
                plot(xt(ii),yt(ii),'ro','MarkerFaceColor','r');
                
                % plotting streakline
                plot(xs(ii,1:is),ys(ii,1:is),'bo','MarkerFaceColor','b');
                axis([-3 3 0 2]);
                
                % plotting streamline
                xlc= -3:0.05:3;
                ylc=yt(ii)*abs(tanh(xt(ii)/L)./tanh(xlc/L));
                plot(xlc,ylc,'k')
                
                %         pause;
                nfilm=['film',num2str(is)];
                print('-dpng',nfilm);
            end;
            
            % new position of the particle at t(ii+1)
            xt(ii+1)=xt(ii)+dt*uf(xt(ii),yt(ii),t(ii));
            yt(ii+1)=yt(ii)+dt*vf(xt(ii),yt(ii),t(ii));
            
            % streakline at t(ii+1)
            for ij=1:(nt-1)/10
                xs(ii+1,ij)=xs(ii,ij)+dt*uf(xs(ii,ij),ys(ii,ij),t(ii));
                ys(ii+1,ij)=ys(ii,ij)+dt*vf(xs(ii,ij),ys(ii,ij),t(ii));
            end;
            
            
        end
        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
    case 2
        %% B. Unsteady Flow L=1+0.5sin(2t)/2
        
        vf = @(x,y,Lnew) -V0*y/cosh(x/Lnew)^2;
        uf = @(x,y,Lnew) V0*Lnew*tanh(x/Lnew);
        
        % xt,yt: trajectory points
        xt=zeros(nt,1);
        yt=zeros(nt,1);
        xt(1) = x0;
        yt(1) = y0;
        
        % xs,ys streamline points
        xs=zeros(nt,(nt-1)/10);
        ys=zeros(nt,(nt-1)/10);
        xs(1:10:end-1,:)=x0;
        ys(1:10:end-1,:)=y0;
        is=0;
        
        figure
        
        for ii=1:nt-1
            
            L(ii) = 1+0.5*sin(2*t(ii));
            
            if mod(ii-1,10)==0
                is=is+1;
                xs(ii,is)=x0;
                ys(ii,is)=y0;
                v = -V0*yy./cosh(xx/L(ii)).^2;
                u = V0*L(ii)*tanh(xx/L(ii));
                
                hold off;
                quiver(xx,yy,u,v,0.8,'k');
                hold on;
                
                % plotting trajectory
                plot(xt(1:ii),yt(1:ii),'ro')
                plot(xt(ii),yt(ii),'ro','MarkerFaceColor','r');
                
                % plotting streakline
                plot(xs(ii,1:is),ys(ii,1:is),'bo','MarkerFaceColor','b');
                axis([-3 3 0 2]);
                
                % plotting streamline
                xlc= -3:0.05:3;
                ylc=yt(ii)*abs(tanh(xt(ii)/L(ii))./tanh(xlc/L(ii)));
                plot(xlc,ylc,'k')
                
                %         pause;
                nfilm=['film',num2str(is)];
                print('-dpng',nfilm);
            end;
            
            % new position of the particle at t(ii+1)
            xt(ii+1)=xt(ii)+dt*uf(xt(ii),yt(ii),L(ii));
            yt(ii+1)=yt(ii)+dt*vf(xt(ii),yt(ii),L(ii));
            
            % streakline at t(ii+1)
            for ij=1:(nt-1)/10
                xs(ii+1,ij)=xs(ii,ij)+dt*uf(xs(ii,ij),ys(ii,ij),L(ii));
                ys(ii+1,ij)=ys(ii,ij)+dt*vf(xs(ii,ij),ys(ii,ij),L(ii));
            end;
            
            
        end
        
        
end