%%    UMD AOSC 615  
%%    Point Vortex Model

%% startup
clear all
close all
cpu0=cputime;
disp(sprintf('\n point vortex model with tracer: %s',datestr(now)))

%% global
global ode_prmtr

%% input
% Number of variables
Nv=2; Nt=1;
% initial 
%- initial condition
%-- vortex
xv0=zeros(2*Nv,1); Gv0=zeros(Nv,1);
n=1; n2=2*n; n1=n2-1; xv0(n1:n2)=[ 1; 0];  Gv0(n)=1;   
n=2; n2=2*n; n1=n2-1; xv0(n1:n2)=[-1; 0];  Gv0(n)=1;
% n=3; n2=2*n; n1=n2-1; xv0(n1:n2)=[ 0; 0.5];  Gv0(n)=1;
%-- tracer
xt0=zeros(2*Nt,1); Gt0=zeros(Nt,1);
n=1; n2=2*n; n1=n2-1; xt0(n1:n2)=[0.5;0];  
% integration time
Tot=300;       %  integration time
dTot=0.1;   %  output
% plot
ixn=[1 3];
c_vt=['r','b'];
lw_vt=[2 1];
axs_vt=[-2 2 -2 2];
c_tx=[c_vt(1), c_vt(2)];
lw_tx=[lw_vt(1), lw_vt(2)];
axs_tx=[0 Tot/2 -2 2];

%% setup
% dimension (two dimnsional position and circulation)
N=3*(Nv+Nt);
% inital condition (concatination of the position & circulation)
x0=[xv0; xt0; Gv0; Gt0]; 
% output 
tot=[0:dTot:Tot];
KTot=length(tot);
xot=zeros(N,KTot);
% model parameter for ode
ode_prmtr=[Nv Nt];    

%% model forward integration
xot(:,1)=x0;
for k=2:KTot;  kn1=k-1; xk0=xot(:,kn1);
    [tk,xk]=ode45(@pvt_dxdt,[tot(kn1) tot(k)],xk0);
    xot(:,k)=xk(end,:)'; 
end

%% plot
% t and selected xk
figure
for n=1:length(ixn); n2=2*ixn(n); n1=n2-1;
    plot(tot,xot(n1,:),c_tx(n),'LineWidth',lw_tx(n)); if(n==1); hold on; end
    plot(tot,xot(n2,:),[c_tx(n) '-.'],'LineWidth',lw_tx(n)); 
end
legend('x_1','x_2','x_5','x_6','Location','NorthEast')
axis(axs_tx)
xlabel('time'); ylabel('x_n')

% xk
figure
for n=1:Nv; n2=2*n; n1=n2-1;
    plot(xot(n1,:),xot(n2,:),c_vt(1),'LineWidth',lw_vt(1));  if(n==1); hold on; end
    plot(xot(n1,1),xot(n2,1),'o','MarkerSize',8,'MarkerFaceColor','w','MarkerEdgeColor',c_vt(1))
end
for n=1:Nt; n2=2*(Nv+n); n1=n2-1;
    plot(xot(n1,:),xot(n2,:),c_vt(2),'LineWidth',lw_vt(2)); 
    plot(xot(n1,1),xot(n2,1),'^','MarkerSize',8,'MarkerFaceColor','w','MarkerEdgeColor',c_vt(2))    
end    
axis(axs_vt); axis square;
xlabel('x_2_n_-_1'); ylabel('x_2_n')


%% finish
disp(sprintf(' ... end of point vortex model (cpu=%s sec)\n',num2str(cputime-cpu0)))
