
clear all
close all
 e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

clear e
global sigma 
global mu

sigma=1;                                  %1 for focussing case
                                           %-1 for defocussing case 
solution_number=5;
%%%%%%%%% Setting Parameters
%sphi                 decay rate for the solution Phi
%s                    decay rate for the eigenvalue/eigenvectors < sphi
%mphi                 finite dimensional parameter for the solution Phi
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%% IMPORT NUMERICAL DATA
cd num_data
if sigma==1
     cd focusing
     if solution_number==1
          alphas=load('alphas_9_0.txt') ;
          Matckdk=importdata('CsDs_9_0.txt') ;
          sphi=4; 
          s=3; 
          mphi=18; 
     elseif solution_number==2
          alphas=load('alphas_99_0.txt') ;
          Matckdk=importdata('CsDs_99_0.txt') ;
          sphi=4; 
          s=3; 
          mphi=24; 
     elseif solution_number==3
          alphas=load('alphas_999_0.txt');
          Matckdk=importdata('CsDs_999_0.txt');
          sphi=3.5; 
          s=3; 
          mphi=30; 
     elseif solution_number==4
          alphas=load('alphas_778_1.txt') ;
          Matckdk=importdata('CsDs_778_1.txt') ;
          Matckdkcomplex=importdata('CsDs_778_1_complex.txt') ;
          sphi=4; 
          s=3; 
          mphi=30; 
     elseif solution_number==5
          alphas=load('alphas_882_1.txt') ;
          Matckdk=importdata('CsDs_882_1.txt');
          Matckdkcomplex=importdata('CsDs_882_1_complex.txt');
          sphi=3.5; 
          s=3; 
          mphi=30; 
     elseif solution_number==6
          alphas=load('alphas_944_1.txt') ;
          Matckdk=importdata('CsDs_944_1.txt');
          Matckdkcomplex=importdata('CsDs_944_1_complex.txt');
          sphi=3.2; 
          s=3; 
          mphi=30; 
     elseif solution_number==7
          alphas=load('alphas_741_2.txt');
          Matckdk=importdata('CsDs_741_2.txt');
          Matckdkcomplex=importdata('CsDs_741_2_complex.txt');
          sphi=3.1; 
          s=3; 
          mphi=44; 
     elseif solution_number==8
          alphas=load('alphas_806_2.txt') ;
          Matckdk=importdata('CsDs_806_2.txt');
          Matckdkcomplex=importdata('CsDs_806_2_complex.txt');
          sphi=3.1; 
          s=3; 
          mphi=58; 
     elseif solution_number==9
          alphas=load('alphas_862_2.txt') ;
          Matckdk=importdata('CsDs_862_2.txt');
          Matckdkcomplex=importdata('CsDs_862_2_complex.txt');
          sphi=3; 
          s=2.8; 
          mphi=80; 
     end
else
   cd defocusing
   if solution_number==1
          alphas=load('defoc_alphas_99_0.txt') ;
          Matckdk=importdata('defoc_CsDs_99_0.txt') ;
          sphi=4; 
          s=3; 
          mphi=20; 
     elseif solution_number==2
          alphas=load('defoc_alphas_999_0.txt') ;
          Matckdk=importdata('defoc_CsDs_999_0.txt') ;
          sphi=3.5; 
          s=3; 
          mphi=30; 
     elseif solution_number==3
          alphas=load('defoc_alphas_9999_0.txt');
          Matckdk=importdata('defoc_CsDs_9999_0.txt');
          sphi=3.1; 
          s=3; 
          mphi=36; 
     elseif solution_number==4
          alphas=load('defoc_alphas_741_1.txt') ;
          Matckdk=importdata('defoc_CsDs_741_1.txt') ;
          sphi=4; 
          s=3; 
          mphi=20; 
     elseif solution_number==5
          alphas=load('defoc_alphas_914_1.txt') ;
          Matckdk=importdata('defoc_CsDs_914_1.txt');
          sphi=3.1; 
          s=3; 
          mphi=30; 
     elseif solution_number==6
          alphas=load('defoc_alphas_972_1.txt') ;
          Matckdk=importdata('defoc_CsDs_972_1.txt');
          sphi=3.2; 
          s=3; 
          mphi=42; 
     elseif solution_number==7
          alphas=load('defoc_alphas_055_2.txt');
          Matckdk=importdata('defoc_CsDs_055_2.txt');
          sphi=4; 
          s=2.8; 
          mphi=16; 
     elseif solution_number==8
          alphas=load('defoc_alphas_640_2.txt') ;
          Matckdk=importdata('defoc_CsDs_640_2.txt');
          sphi=3.5; 
          s=2.8; 
          mphi=54; 
     elseif solution_number==9
          alphas=load('defoc_alphas_826_2.txt') ;
          Matckdk=importdata('defoc_CsDs_826_2.txt');
          sphi=3; 
          s=2.8; 
          mphi=54; 
     end
   
   
   
 %  alphas=load('defoc_alphas_9999_0.txt') ;
 %  Matckdk=importdata('defoc_CsDs_9999_0.txt') ; 
end
cd ..
cd ..

%%%%%%%%% SETTING the data

number_nodes=alphas(end);                 % indicates the parity/ number of nodes
mu=alphas(end-1);                         % value of mu
alphas(end-1:end)=[];                     % Fourier coefficient of the num. solu
% the zero entries of alphas  due to symmetry  have been already deleted.
parity=mod(number_nodes,2);               % symmetry of the solution Phi
%%%%%%%%%%%%%%%%% Parameters

m=3*mphi;
M=m+4*mphi;
power=12;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp([' mu = ' num2str(mu)])
disp([' sigma = ' num2str(sigma)])



b=(alphas)/2;                     % vector b= bodd or beven 
b=[b;zeros(mphi-length(b),1)];    % make b to be mphi lenght by appending or removing items
b(mphi+1:end)=[];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(' Computation of Phi') %%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%
%%% only odd terms - ground state

if parity==0           % even j means odd coefficients
   disp('Refine numerical solution')
    ris=newton_phi_odd(b,power);
   
   disp('Compute enclosure of the bound state Phi')
   r = radii_polynomials_phi_odd(ris,sphi);
   maxb=max(abs(ris));
   ris_int=midrad([ris;zeros(M-length(ris),1)],r*[1:1:M].'.^(-sphi));

   bint=intval(zeros(2*size(ris_int,1),1));
   bint([1:2:size(bint,1)],1)=ris_int;                   % vector b_k for k>=1
else    %parity=1   odd solution = even coefficient
     disp('Refine numerical solution')
    ris=newton_phi_even(b,power);
     disp('Compute enclosure of the bound state Phi')
    r = radii_polynomials_phi_even(ris,sphi);
    maxb=max(abs(ris));
    ris_int=midrad([ris;zeros(M-length(ris),1)],r*[1:1:M].'.^(-sphi));

    bint=intval(zeros(2*size(ris_int,1),1));
    bint([2:2:size(bint,1)],1)=ris_int;               % vector b_k for k>=1

end
rphi=r;				  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Definition of the matrix F_nl
F=intval(zeros(m));
Bint=[0;bint];           % vector b_k for k>=0

Quadb=real(quadratic_sumFFT_int_sign(Bint,Bint,-1));
for l=1:m
    n=l:m;
    n(mod(n+l,2)==1)=[]; % if n+l is odd also n-l is odd and k,p in the sum are one odd and one even. Hence the sum is always zero
    F(n,l)=Quadb(n+l+1)-Quadb(n-l+1);
    F(l,n)=F(n,l);
end
F=intval(2)*F;
F=F+intval(8)*intval(2)^intval(sphi)*intval(r)*(intval(maxb)+intval(r))/((intval(sphi)-intval(1))*intval(mphi)^(intval(sphi)-intval(1))); % takes into account the tail of b_k

%%%%%%%%%%%%%%%%%%
%%% Enclosure eigenvalues/eigenvectors
disp('Enclosure of the Eigenvalues/eigenvectors')
Results=[];
Eigenvalue=[];
Eigenvector=[];
radius=[];

Matckdk=mat2cell(Matckdk,size(Matckdk,1),ones(size(Matckdk,2),1));
if exist('Matckdkcomplex')
    Matckdk=[Matckdk,mat2cell(Matckdkcomplex,size(Matckdkcomplex,1),ones(size(Matckdkcomplex,2),1))];
end
for number_eigenvalue=1:size(Matckdk,2)

ckdk=Matckdk{1,number_eigenvalue};
num_beta=ckdk(end-2);

N=(length(ckdk)-3)/2;
ck=[ckdk(1:N,1);zeros(m-N,1)];
dk=[ckdk(N+1:2*N,1);zeros(m-N,1)];
num_data=reshape([ck,dk]',[],1);


num_sol=[num_beta;num_data(1:end)];
power=12;

[maxc_k,index_max]=max(abs(ck));

x=newton_NLS(m,F.mid,num_sol,power,index_max);

display(['Eigenvalue=' num2str(x(1))]);
I2=NLS_control_int(m,M,mphi,F,x,index_max,Bint,rphi,s,sphi);


r2=I2(1);
Eigenvalue=[Eigenvalue,x(1)];
Eigenvector=[Eigenvector,x(2:end)];
radius=[radius,r2];
testGamma=sum((-1).^[1:1:m].'.*[1:1:m].'.*(ck(1:1:m)-dk(1:1:m)));
RadiusBeta=2*r2*(1+1/(s-2));
Results=[Results,[I2(1);testGamma;RadiusBeta]];
end

%%%% PLOT the bound state %%%%%%%
falphak=2*bint.mid;
X=linspace(0,1,1000);
phi=falphak(1)*sin(pi*1*X);
for k=2:mphi
    phi=phi+falphak(k)*sin(pi*k*X);
end
plot(X,phi);
title('Bound state')

%%%%%% PLOT eigenfunctions %%%%%
  %%% choose different number eigenvalue if needed 
% ckdk=Matckdk{1,number_eigenvalue};
% N=(length(ckdk)-3)/2;
% ck=[ckdk(1:N,1);zeros(m-N,1)];
% dk=[ckdk(N+1:2*N,1);zeros(m-N,1)];

wx=ck(1)*sin(pi*1*X);
zx=dk(1)*sin(pi*1*X);
for k=2:size(ck,1)
    wx=wx+ck(k)*sin(pi*k*X);
    zx=zx+dk(k)*sin(pi*k*X);
end
figure
plot(X,wx);
hold on
plot(X,zx,'r');
title('Function w(x) in blue and z(x) in red')



% SIMPLICITY of EIGENVALUE
% change 1 to any other eigenvalue number 
rad=radius(1);
V=Eigenvector(:,1);
Beta=midrad(Eigenvalue(1),rad);
V=[V(1:2:end,:),V(2:2:end,:)];
Vint=midrad(V,r*[1:1:size(V,1);1:1:size(V,1)].'.^(-s));
m=size(Vint,1);
M=max(m+1,4*mphi+m+1);
V=reshape(Vint',[],1);
Xbar=[0.01;0;zeros(2*m-1,1)];
INTsimpl=checksimpleig(Beta,V,m,M,Xbar,mphi,F,Bint,rphi,s,sphi);


%%%%%%%%% DISPLAY RESULTS %%%%
format short e
disp(' ')
disp('Some results')
disp(' ')

disp(' The bound state is computed with')
disp([' decay rate parameter s=' num2str(sphi)])
disp([' enclosure radius r=' num2str(rphi)])
disp(' ')
disp(['The parameter mu is ' num2str(mu)]);
disp(' ')
display('Eigenvalues ')
display(num2str(Eigenvalue));
disp(['are computed with decay rate parameter s=' num2str(s)])
display('and  enclosure radius  ')
display(num2str(Results(1,:)));
disp(' ')
display('Value of Gamma')
display(num2str(Results(2,:)));
display('enclosure of Gamma ')
display(num2str(Results(3,:)));
