format compact
set(groot,'defaulttextinterpreter','latex');
set(groot,'defaultAxesTickLabelInterpreter','latex');
set(groot,'defaultLegendInterpreter','latex');
set(groot,'defaultAxesFontSize',14);

mu=2
Omega=0.1
Rtf=sqrt(2*mu)/Omega;

S = 1;               % charge of vortex;
xx1=0.25*Rtf;yy1=0;  % position to seed vortex

Nx=251;Ny=251;Lx=Rtf*1.5;Ly=Lx;

% UINI:
x=linspace(-Lx,Lx,Nx);
y=linspace(-Ly,Ly,Ny);
skip=1;
[X,Y]=meshgrid(x,y);

%compute TF profile using BVP
fprintf('\nComputing TF with BVP\n');
[r,u,up] = SteadyStateBVP(mu,0,Omega);
r=r(1:skip:end); u=u(1:skip:end);
RTF=r;RTF=[RTF,linspace(r(end),1000,100)];
UTF=u;UTF=[UTF,zeros(1,100)];

%compute vortex profile at (0,0) using BVP
fprintf('\nComputing vortex with BVP\n');
[r,u,up] = SteadyStateBVP(mu,abs(S),0);
%u=u./max(u);
RS1=r;RS1=[RS1,linspace(r(end),1000,100)];
US1=u;US1=[US1,ones(1,100)*max(u)];
NN=100;
theta=linspace(0,2*pi,NN);
[TTTF,RRTF]=meshgrid(theta,RTF);
[TTS1,RRS1]=meshgrid(theta,RS1);
[XTF,YTF]=pol2cart(TTTF,RRTF);
[XS1,YS1]=pol2cart(TTS1,RRS1);
UUTF=kron(ones(1,NN),UTF');
UUS1=kron(ones(1,NN),US1');
uuTF=griddata(XTF,YTF,UUTF,X,Y);

clear TTTF RRTF TTS1 RRS1 XTF YTF UUTF RS1 US1 UTF RTF r u up theta NN

  %imprint vortex#1
  uuS1=griddata(XS1+xx1,YS1+yy1,UUS1,X,Y);
  ANG=atan2(Y-yy1,X-xx1);
  uuS1=uuS1.*exp(1i*ANG*S);
  uini=uuTF.*uuS1./(max(max(uuS1)));

% secondvortex=1;
% if(secondvortex) %imprint vortex#2
%  xx2=-xx1;yy2=0;
%  uuS2=griddata(XS1+xx2,YS1+yy2,UUS1,X,Y);
%
%  %removing NaN from uuS2
%  A2=1-(abs(uuS2).^2<1e9)*1;
%  [ii,jj]=find(A2);
%  for iii=1:length(ii)
%   ix=ii(iii);
%   iy=jj(iii);
%   xdisp=[-1,1,0,0]; ydisp=[0,0,-1,1];Asum=[];
%   for jjj=1:4
%    uudisp=uuS2(ix+xdisp(jjj),iy+ydisp(jjj));
%    if(uudisp<1e99) Asum=[Asum,uudisp];end
%   end
%   %uuS2(ix,iy)=mean([uuS2(ix-1,iy),uuS2(ix+1,iy),uuS2(ix,iy-1),uuS2(ix,iy+1)]);
%   uuS2(ix,iy)=mean(Asum);
%   %uuS2(ix,iy)
%  end
%
%  ANG=atan2(Y-yy2,X-xx2);
%  uuS2=uuS2.*exp(1i*ANG*S);
%  uini=uini.*uuS2./(max(max(uuS2)));
% end % if(secondvortex)

clear ANG A2 XS1 YS1 UUS1 ii iii jj ix iy xdisp uudisp secondvortex skip

figure(12)

subplot(2,2,1)
surf(X,Y,abs(uuS1).^2)
shading interp;view([0 90]);axis tight;
title('$|u_{vortex}|^2$');xlabel('x'); ylabel('y');

subplot(2,2,2)
surf(X,Y,atan2(imag(uuS1),real(uuS1)))
shading interp;view([0 90]);axis tight;
title('$phase(u_{vortex})$');xlabel('x'); ylabel('y');

subplot(2,2,3)
surf(X,Y,uuTF)
shading interp;view([0 90]);axis tight;
title('$u_{TF}$');xlabel('x'); ylabel('y');

subplot(2,2,4)
surf(X,Y,abs(uini).^2)
shading interp;view([0 90]);axis tight;
title('$|u_{TF} * u_{vortex}|^2$');xlabel('x'); ylabel('y');