% Finding a single vortex steady state using radial ODE and BVP solver:
%       i u_t = -(1/2) (u_xx+u_yy) + myalpha*(|u|^2) u + V(x) u
% myalpha = 1 is defocusing
% myalpha =-1 is focusing
%
% Examples:
%
% Ex1: Vortex inside MT
% S=1;mu=1;Omega=0.05;
% [r,u,up] = SteadyStateBVP(mu,S,Omega);
%
% Ex2: Vortex on a homogeneous background
% S=1;mu=1;Omega=0;
% [r,u,up] = SteadyStateBVP(mu,S,Omega);
%

function [r,u,up] = SteadyStateBVP(mu,S,Omega)

format compact
plotsteadystates=1;
plotderivs=0;

%mu=132;
%S = 1; % vortex charge;
%Omega=0*0.2;
myalpha=+1; 		% (+1:defocusing. -1:focusing)

r0 = 1e-8;
if(Omega)
 Rtf=sqrt(2*(mu))/Omega;
 rmax =2.0*Rtf;
else
 rmax = 35;
end

Nr = 400;
rini=linspace(r0,rmax,Nr);
drini=rini(2)-rini(1);

fprintf('S    =\t%f\n',S);
fprintf('myalpha=\t%f\n',myalpha);
fprintf('Omega=\t%f\n',Omega);
fprintf('mu   =\t%f\n',mu);

VMT = 0.5*Omega^2*rini.^2;
V=VMT;

utf = max(0,real(sqrt((mu-V)./myalpha)));

if(Omega==0)
 u0=sqrt(mu)+rini*0;
else
 u0=utf;
end

if(Omega)
 smooth=120;
else
 smooth=0;
end
for ii=1:smooth;
 u0p=[u0(1),u0(1:end-1)];
 u0m=[u0(2:end),0];
 u0=(u0p+u0+u0m)./3;
end

if(S)
 widthtanh=0.8*(mu/myalpha)^0.5/S;
 u0=u0.*tanh(widthtanh*rini).^sqrt(abs(S));
end

%u0pold =[0,(u0(3:end)-u0(1:end-2))/(2*drini),0];
%4th order:
u0p =[0,(u0(3)-u0(1))/(2*drini),...
    (-u0(5:end)+8*u0(4:end-1)-8*u0(2:end-3)+u0(1:end-4))./(12*drini),...
    (u0(end)-u0(end-2))/(2*drini),0];
if(S)
 u0p(1)=widthtanh*sqrt(mu/myalpha);
else
 u0p(1)=0;
 u0(1)=sqrt(mu);
end


solinit.x=rini;
solinit.y(1,:) = u0;
solinit.y(2,:) = u0p;

yprime0 = vortexode(rini,[u0;u0p],mu,myalpha,Omega,S);

%figure(10)
%subplot(2,1,1)
%plot(rini,u0,'o-',rini,u0p,'o-')
%subplot(2,1,2)
%plot(rini,yprime0,'o-')
%caca

ms=1.5;

if(plotsteadystates)
%figure(11), clf 
%plot(r,real(u0),r,imag(u0),r,real(u0p),r,imag(u0p)); 
%plot(rini,u0,':',rini,u0p,':'); 

subplot(2,1,1)
plot(rini,utf.^2,'-',rini,abs(u0).^2,'s--',rini,V,'r-','MarkerSize',ms); 
ax=[0 rmax -0.1*mu 1.10*mu];
axis(ax);
xlabel('r');
ylabel('$u_{0}^2$, $|u|^2$')
hold on
title(['$\alpha=$',num2str(myalpha),', $\Omega=$',num2str(Omega),', $\mu$=',num2str(mu),', $S=$',num2str(S)])
legend('$u_{TF}^2$','$u_{0}^2$','$V$')
drawnow
end %plotsteadystates

%options=bvpset;
options=bvpset('AbsTol',1e-10,'RelTol',1e-5,'Nmax',1000);
%options=bvpset('AbsTol',1e-12,'RelTol',1e-6,'Nmax',1000);
%options=bvpset('AbsTol',1e-15,'RelTol',1e-7,'Nmax',2000);
%options=bvpset('Nmax',1000);
%options=[];
%options

sol = bvp4c(@vortexode,@vortexbc,solinit,options,mu,myalpha,Omega,S);

r  = sol.x;
u  = sol.y(1,:);
up = sol.yp(1,:);
u2  = sol.y(2,:);
up2 = sol.yp(2,:);
u(end)=u(end-1);
up(end)=up(end-1);
u0n   = spline(rini,u0,r);
N=length(r);

computemass=1;
if(computemass)
 hbar = 1.05457148e-34;          %[Kg m^2/ s]
 amu=1.660538e-27;               %[Kg]
 mass = 87*amu;                  %[Kg]
 %ascatter = 55.4e-10;           %[m] (old value)
 ascatter = 54.5e-10;         %[m]

 fz=101.2;                       %[Hz]
 fr= 35.8;                       %[Hz]
 wz = 2*pi*fz;                   %[rad/s]
 wr = 2*pi*fr;                   %[rad/s]

 xi = sqrt(hbar/(mass*wz));      %[m] 
 az =  xi;                       %[m] 
 g3 = 4*pi*hbar^2*ascatter/mass;
 g2 = g3/(sqrt(2*pi)*az);
 rescM = hbar^2/(mass*g2);

 %dr=diff(r);
 %M2=sum(u(2:end).^2.*r(2:end).*2.*pi.*dr)
 %M3=rescM*M2
end

if(plotsteadystates)
%figure(11)
%plot(r,u,r,up,r,u2); 
%plot(r,up1-u2); 
%plot(r,real(u),r,imag(u),r,abs(u),'k-'); 

subplot(2,1,1)
plot(r,abs(u).^2,'ko-','MarkerSize',ms); 
%plot(r,r*0,'k-'); 
legend('$u_{TF}^2$','$u_{0}^2$','$V$','$|u|^2$')

subplot(2,1,2)
plot(r,real(u),'o-',r,imag(u),'o-',...
     r,real(u0n),'--',r,imag(u0n),'--','MarkerSize',ms); 
xlabel('r');
ylabel('Re(u), Im(u)')
ax=[0 rmax -0.1*sqrt(mu) 1.10*sqrt(mu)];
axis(ax);
legend('Re(u)','Im(u)','Re(u0)','Im(u0)')

%subplot(3,1,3)
%plot(r,unwrap(angle(u)),'o-',...
%     r,r*0,'k-','MarkerSize',ms); 
%xlabel('r');
%ylabel('Arg(u)')

hold off
drawnow
end %plotsteadystates

if(plotderivs)

% Residue:
%dr2 = r(3:end)-r(1:end-2);
%u0np  = [0,(u0n(3:end)-u0n(1:end-2))./(2*dr),0];
%4th order:
u0np =[0,(u0n(3)-u0n(1))/(r(3)-r(1)),...
    (-u0n(5:end)+8*u0n(4:end-1)-8*u0n(2:end-3)+u0n(1:end-4))./...
    (6*(r(4:end-1)-r(2:end-3))),...
    (u0n(end)-u0n(end-2))/(r(end)-r(end-2)),0];

%u0npp=[0,(u0n(3)-2*u0n(2)-u0n(1))/((r(2)-r(1))*(r(3)-r(2))),...
%    (-u0n(5:end)+16*u0n(4:end-1)-30*u0n(3:end-2)+16*u0n(2:end-3)-u0n(1:end-4))./...
%    (12*(r(3:end-2)-r(2:end-3)).*(r(4:end-1)-r(3:end-2))),...
%    (u0n(end)-2*u0n(end-1)+u0n(end-2))/((r(end-1)-r(end-2))*(r(end)-r(end-1))),0];
u0npp =[0,(u0np(3)-u0np(1))/(r(3)-r(1)),...
    (-u0np(5:end)+8*u0np(4:end-1)-8*u0np(2:end-3)+u0np(1:end-4))./...
    (6*(r(4:end-1)-r(2:end-3))),...
    (u0np(end)-u0np(end-2))/(r(end)-r(end-2)),0];
yprime0 = vortexode(r,[u0n;u0np],mu,myalpha,Omega,S);

upp1 =[0,(up(3)-up(1))/(r(3)-r(1)),...
    (-up(5:end)+8*up(4:end-1)-8*up(2:end-3)+up(1:end-4))./...
    (6*(r(4:end-1)-r(2:end-3))),...
    (up(end)-up(end-2))/(r(end)-r(end-2)),0];
%upp2=[0,(u(3)-2*u(2)-u(1))/((r(2)-r(1))*(r(3)-r(2))),...
%    (-u(5:end)+16*u(4:end-1)-30*u(3:end-2)+16*u(2:end-3)-u(1:end-4))./...
%    (12*(r(3:end-2)-r(2:end-3)).*(r(4:end-1)-r(3:end-2))),...
%    (u(end)-2*u(end-1)+u(end-2))/((r(end-1)-r(end-2))*(r(end)-r(end-1))),0];
%upp = (up(3:end)-up(1:end-2))./dr;upp = [0,upp,0];
yprime = vortexode(r,[u;up],mu,myalpha,Omega,S);

figure(13)

ind=52:N;

whos r u0npp yprime upp2 

subplot(4,1,1)
plot(r(ind),(u0npp(ind)),'-','MarkerSize',ms);

hold on
plot(r(ind),(yprime0(2,(ind))),'s','MarkerSize',ms);
plot(r,r*0,'k-'); 
hold off
xlabel('r');
ylabel('u_0'''' & rhs')
axis tight
legend('u_0''''','rhs',1)
title('Residue')

subplot(4,1,2)
plot(r(2:end),real(u0npp(2:end)-yprime0(2,2:end)),'s-',...
     r,r*0,'k-','MarkerSize',ms);
xlabel('r');
ylabel('u_0''''-rhs')
axis tight

subplot(4,1,3)
%plot(r(ind),up(ind),'o-',r(ind),yprime(2,(ind)),'o-','MarkerSize',ms);
%plot(r(ind),u(ind),'--',r(ind),up(ind),'--','MarkerSize',ms);
%hold on
%plot(r(ind),yprime(1,(ind)),'s-',r(ind),yprime(2,(ind)),'s-','MarkerSize',ms);
%hold off
plot(r(ind),(upp1(ind)),'-');
hold on
%plot(r(ind),real(upp2(ind)),'o',r,imag(upp2(ind)),'o','MarkerSize',ms);
plot(r(ind),(yprime(2,(ind))),'s','MarkerSize',ms);
plot(r,r*0,'k-'); 
hold off
xlabel('r');
ylabel('u'''' & rhs')
axis tight
%legend('Re(u_1'''')','Im(u_1'''')','Re(u_2'''')','Im(u_2'''')','Re(rhs)','Im(rhs)',1)
legend('u''''','rhs',1)
title('Residue')

subplot(4,1,4)
plot(r(ind),upp1(ind)-yprime(2,(ind)),'s-',r(ind),r(ind)*0,'k-','MarkerSize',ms);
xlabel('r');
ylabel('u''''-rhs')
axis tight

end %plotderivs


%---------Boundary Condition--------%
function res=vortexbc(ua,ub,mu,myalpha,Omega,S)

%ua:left, ub:right, ua(1);u, ua(2):u'

if(Omega==0)
 %res=[ua(1);ub(1)-sqrt(mu/myalpha)];  %BC on height=sqrt(mu/myalpha)
 res=[ua(1);ub(2)];                 % BC on u'(right)=0
else
 res=[ua(1);ub(1)];
end

if(S==0)
 res=[ua(2);ub(1)];
end

%---------ODE function-------------%
function yprime = vortexode(r,u,mu,myalpha,Omega,S)

VMT = 0.5*Omega^2*r.^2;
V=VMT;
u2=u(1,:).^2;
I=sqrt(-1);
onehalf=1/2;
oneoveronehalf=1/onehalf;

yprime=[
    u(2,:); ...
    -u(2,:)./r+S^2*u(1,:)./(r.^2)-oneoveronehalf*(mu-V-myalpha*(u2)).*u(1,:)...
    ];