%
% Driver to apply Newton iteration to find a COMPLEX
% steady state solution for |u|^2 for the GPE:
%       i u_t = -(1/2) (u_xx+u_yy) + alpha*|u|^2 u + V(x) u
%  alpha = 1 is defocusing
%  alpha =-1 is focusing
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% setup
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear

 format compact
 N=100;
 it=0			% Newton iterations
 factsize=1.5;		% how larger is domain wrt Rtf
 
 fs=16;
 mu=1;

 use_old_IC=0;

 Omega=0.05;        % Modulationally unstable
 Omega=0.30;        % Modulationally almost stable (dt=0.05)
 Rtf=sqrt(2*mu)/Omega;
 L=Rtf*factsize;

 xx=1:N^2;
 xxxx=1:2*N^2;
 dx=2*L/N;
 x=linspace(-L,L-dx,N);
 y=linspace(-L,L-dx,N);
 Lx=L;Ly=L;
 [X,Y] = meshgrid(x,y);

 onehalf = 0.5;         % NLS: (1/2)*u_xx or (1)*u_xx
 onehalfoodx2=onehalf/dx^2;
 allAUmf = [];

 plottop=1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GENERATING ICs AND POTENTIAL
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 alpha=1; 		% (+1:defocusing. -1:focusing)

 %MT:
 xV0=0*Rtf/2; yV0=0*Rtf/2;
 Omegax=1*Omega;
 Omegay=1*Omega;
 VMT = Omegax^2*(X-xV0).^2/2+Omegay^2*(Y-yV0).^2/2;

 V = VMT;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GENERATING ICs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(use_old_IC==0)
 utf = sqrt(max(1e-10,mu-V));
 Maxu = max(utf);
 u=utf;

% it=5;newton2d_complex;it=0;u=utf;

%%adding a vortex:
% sigma =+1;
% x0=1*Rtf/4; y0=-1*Rtf/4;
% myangle = atan2((Y-y0),(X-x0));
% u = utf.*tanh(sqrt((Y-y0).^2+(X-x0).^2));
% u = u.*exp(i*sigma*(myangle));

 else
  load('unst_N100.mat','allu','u00');
  u = squeeze(allu(end,:,:));
  clear allu;
 end

%Genratinf a Dark Soliton stripe
 u = u.*tanh(X);
 it=5;newton2d_complex;

 Nx=N;Ny=N;Lx=L;Ly=L;fit=0;saveframes=200;t0=0;disp=50;
 dt=0.1;
 maxtime=100;

% bec_2d
 %film