format compact

myplot=1;  % plot trajectories?
myanim=1;  % animane lattice?
fs=18;

% params;
%A = 0.5;
%Ns = length(Y0);
params = [A,Ns,L];


%%%%%%%%%%%% If there is a trap, undo the following:
% Omega=0.05;
% params= [A,Omega];

%%%%%%%%%%%% Initial Conditions %%%%%%%%%%%%

%%% Two Bright solitons
% Xi0 = [-5 5];
% Vel0 = zeros(1,Ns);

%%% Four Bright solitons
% Xi0 = [-38 -1 15 45];
% Vel0 = zeros(1,Ns);

% Ns BSs
Xi0 = Y0;
% Vel0 = vini;
% Vel0 = zeros(1,Ns);
% Xi0 = -L+0.5*2*L/Ns + [0:Ns-1]*2*L/Ns;

IC = [Xi0,Vel0];

TLperiod=2*L/(velTL*r0);
tf=1000;
tf = 2*TLperiod;
tspan = linspace(0,tf,1000);

[T,Y] = ode45(@Ns_sol_ode_rcg,tspan,IC,[],params);

if(myplot)
 figure(100)
 subplot(2,1,1)
 plot(T,Y(:,1:Ns),'k-')
 xlabel('t');ylabel('${\xi}_n(t)$','Interpreter','Latex','FontSize',fs);
 axis tight
 subplot(2,1,2)
 %plot(T,diff(Y(:,1:Ns)'),'k-')
 DX=kron(linspace(1,Ns,Ns)*max(max(Y(:,[1:Ns]+Ns))),ones(length(T),1));
 plot(T,Y(:,[1:Ns]+Ns)+DX,'k-')
 xlabel('t');ylabel('$\dot{\xi}_n(t) + DX$','Interpreter','Latex','FontSize',fs);
 axis tight
end %(myplot)

if(myanim)
 figure(101)
 Nshow=100; % number of snapshots to show
 skip=round(length(tspan)/Nshow);
 myind=1:Ns;
 myind2=2:Ns;
 mm=min(min(Y(:,myind)));
 MM=max(max(Y(:,myind)));
 mm2=min(min(diff(Y(:,myind)')))*0;
 MM2=max(max(diff(Y(:,myind)')))*1.1;
 %mm=-L;MM=+L; % to be used if plotting mod

 it=1;
 subplot(2,1,1)
 plot(Y(it,myind),Y(it,myind)*0,'o')
 ylabel('${\xi}_n(t)$','Interpreter','Latex','FontSize',fs);axis tight
 axis([mm MM -1 1])
 set(gca,'FontSize',fs)
 text(MM-0.4*(MM-mm),0.6,['$t=',num2str(T(it),3),'$'],'Interpreter','Latex','FontSize',fs);
 subplot(2,1,2)
 plot(Y(it,myind2),diff(Y(it,myind)),'o-')
 xlabel('$x$','Interpreter','Latex','FontSize',fs);
 ylabel('$r_n(t)$','Interpreter','Latex','FontSize',fs);axis tight
 axis([mm MM mm2 MM2])
 set(gca,'FontSize',fs)
 drawnow
 
 fprintf('Press any key to see animation of TL\n');pause
 for it=1:skip:length(T);
   subplot(2,1,1)
   %plot(mod(Y(it,myind)+L,2*L)-L,Y(it,myind)*0,'o')
   plot(Y(it,myind),Y(it,myind)*0,'o')
   ylabel('${\xi}_n(t)$','Interpreter','Latex','FontSize',fs);axis tight
   set(gca,'FontSize',fs)
   axis([mm MM -1 1])
   text(MM-0.4*(MM-mm),0.6,['$t=',num2str(T(it),3),'$'],'Interpreter','Latex','FontSize',fs);
   subplot(2,1,2)
   %plot(mod(Y(it,myind2)+L,2*L)-L,diff(Y(it,myind)),'o-')
   plot(Y(it,myind2),diff(Y(it,myind)),'o-')
   xlabel('$x$','Interpreter','Latex','FontSize',fs);
   ylabel('$r_n(t)$','Interpreter','Latex','FontSize',fs);axis tight
   axis([mm MM mm2 MM2])
   set(gca,'FontSize',fs)
   drawnow
 end
end %(myanim)