% This code is for set the initial condition of BSs and its Toda Lattice,
% It comes from the formula in Paper "Non-linear Wave propagation in the one dimensional Toda lattice"
% By Itai Afek.

format compact
fs=18;

global Y0 Vel0;
% % The amplitude of the chain of solitons and how many solitons are there in
% % the chain.
% A=input('The equal amplitude of solitons =');
% Ns=input('The number of solitons construct the TL =');
A=1;
Ns=20;

% Build the matrix to store the initial data
nspan=1:Ns;
nspan2=2:Ns; % to plot rn's
Y0=zeros(1,Ns);
Vel0=zeros(1,Ns);

% % Arbitrary constant:
% y0=input('You wanna put the first soliton here = ');
% v0=input('You wanna give the first soliton velocity as =');
% alphaTL=input('You wanna set the 1/width of Toda lattice to be =');
y0=0;
v0=0;
r0=10; %separation of homogeneous state;
alphaTL=1.0;

% Other pamameters needed:
B=4*A*A*A;b=A;  % for BRIGHT soliton
%mu=1;B=8*mu^(3/2);b=sqrt(mu);  % for DARK soliton
Bbar=B*exp(-b*r0);
a=Bbar;
betaTL=sqrt(a*b).*sinh(alphaTL);
velTL=betaTL/alphaTL
rmin=-(1/b)*log(betaTL^2/(a*b)+1)

% Positions:
No2=floor(Ns/2);
Y0(No2)=y0;
for ii=No2:Ns-1 % forward positions from No2
    Z=alphaTL.*(ii+1-No2);
    rn0p1 = -(1/b).*log((betaTL.^2/(a*b)).*sech(Z).^2+1);
    Y0(ii+1)=Y0(ii)+r0+rn0p1;
end 
for ii=No2:-1:2 % backward positions from No2
    Z=alphaTL.*(ii+0-No2);
    rn0   = -(1/b).*log((betaTL.^2/(a*b)).*sech(Z).^2+1);
    Y0(ii-1)=Y0(ii)-r0-rn0;
end 

% Velocity:
Vel0(No2)=v0;
for ii=No2:Ns-1 % forward velocities from No2
    Z=alphaTL.*(ii+1-No2);
    rn0p1 = -(1/b).*log((betaTL.^2/(a*b)).*sech(Z).^2+1);
    Sn0p1 = (betaTL^2/(a*b))*sech(Z).^2;
    Tn0p1 = -((2*betaTL)/(b)).*tanh(Z);
    dvel  = Sn0p1.*Tn0p1.*exp(b*rn0p1);
    Vel0(ii+1)=Vel0(ii)+dvel;
end 
for ii=No2:-1:2 % backward velocities from No2
    Z=alphaTL.*(ii+0-No2);
    rn0 = -(1/b).*log((betaTL.^2/(a*b)).*sech(Z).^2+1);
    Sn0 = (betaTL^2/(a*b))*sech(Z).^2;
    Tn0 = -((2*betaTL)/(b)).*tanh(Z);
    dvel  = Sn0.*Tn0.*exp(b*rn0);
    Vel0(ii-1)=Vel0(ii)-dvel;
end 

%substract background veolcity:
Vel0=Vel0-Vel0(end);
%Vel0=Vel0-mean(Vel0);
%Vel0=Vel0-velTL*r0-mean(Vel0);

%domain for periodic lattice:
m=Y0(1)-(r0/2);
M=Y0(end)+(r0/2);
L=(M-m)/2;
Y0=Y0-m-L;


rn  = diff(Y0);
rnp = diff(Vel0);

figure(1);clf;
subplot(2,2,1)
plot(nspan,Y0,'o-')
xlabel('$n$','Interpreter','Latex','FontSize',fs);
ylabel('$\xi_n(0)$','Interpreter','Latex','FontSize',fs);axis tight
set(gca,'FontSize',fs)
subplot(2,2,3)
plot(nspan2,rn-r0,'o-')
xlabel('$n$','Interpreter','Latex','FontSize',fs);
ylabel('$r_n(0)$','Interpreter','Latex','FontSize',fs);axis tight
set(gca,'FontSize',fs)

subplot(2,2,2)
plot(nspan,Vel0,'o-')
xlabel('$n$','Interpreter','Latex','FontSize',fs);
ylabel('$\dot{\xi}_n(0)$','Interpreter','Latex','FontSize',fs);axis tight
set(gca,'FontSize',fs)
subplot(2,2,4)
plot(nspan2,rnp,'o-')
xlabel('$n$','Interpreter','Latex','FontSize',fs);
ylabel('$\dot{r}_n(0)$','Interpreter','Latex','FontSize',fs);axis tight
set(gca,'FontSize',fs)