clear all; clf;  close all;
k = 9;
rmax = 50;
% nbreak exponential brakepoints 
nbreak = 20; 
i = linspace(1,nbreak,nbreak);
gamma = 5.0 ;
breaks = rmax*(exp(gamma.*(i-1)./(nbreak-1))-1)/(exp(gamma)-1);
%breaks = linspace(0,rmax,nbreak);
%plot(breaks,zeros(nbreak,1),'bo')
% knot sequence, multiplicities k,1,1,1 ...
multip = ones(nbreak,1);
multip(1) = k;
multip(nbreak)=k;
n = sum(multip)-k;
% n+k knots
ii = 0;
for i = 1:nbreak
    for j = 1:multip(i)
        ii = ii + 1;
        t(ii) = breaks(i);
    end 
end 




%x=linspace(breaks(1),breaks(nbreak),1000);
%figure (1)
%plot(breaks,zeros(nbreak,1),'bo');hold on;
%for i = 1:n
%    for ix = 1:size(x,2)
%        y(ix)=Bi(x(ix),i,t,n,k);
%    end
%    plot(x,y,'Color',rand(1,3))
%end


% Gauss_legendre integrals
% 

nGL=10;  % # of GL points
S=zeros(n,n);
H1=zeros(n,n);
H2=zeros(n,n);
H3=zeros(n,n);
for seg = 1: nbreak-1  % t-segments
    a = breaks(seg);
    b = breaks(seg+1);
    %disp([num2str(a) ' -> ' num2str(b) ])
         
    [xi,wi]=gauleg(a,b,nGL); % abscissae and weights
    [B,Bp,Bpp,left]=Bspline(xi,t,n,k); 
    left= left-k+1; % shift because k first values are skipped in left
    for ir = 1: nGL
        %disp([num2str(xi(ir)) ' ' num2str(ir) ' ' num2str(left(ir)) ' ' num2str(left(ir)+k-1) ] )
        for i = left(ir): left(ir)+k-1
            if(i>n) continue; end
            for j = left(ir): left(ir)+k-1
                if (j>n) continue; end
                ii = i-left(ir)+1;
                jj = j-left(ir)+1;
                S(i,j) = S(i,j) + wi(ir)* B(ir,ii)*B(ir,jj);
                H1(i,j) = H1(i,j) -  wi(ir)* B(ir,ii)*Bpp(ir,jj)/2;
                H2(i,j) = H2(i,j) +  wi(ir)* B(ir,ii)*(-1/xi(ir))*B(ir,jj);
            end 
        end 
    end
end
S;
H1;
H2;

H = H1+H2;

% drop out first and last B_i(r) from solution 
% to satisfy boundary condition R(0)=0
[c,E]=eig(H(2:n-1,2:n-1),S(2:n-1,2:n-1));

E=diag(E);
[E, order] = sort(E);
cc = c;
for i = 1: n-2
    c(:,i) = cc(:,order(i)); 
end
%check that sorting went ok
%state = 1;
%[H(2:n-1,2:n-1)*c(:,state) E(state)*S(2:n-1,2:n-1)*c(:,state)]

for i=1:min(10,size(E))
    format long;
    exact = -1/2*1/i^2;
    disp(['i= ',sprintf('%5i',i) '   E = ',sprintf('%25.20f',E(i)), ...
         '   Exact = ',sprintf('%25.20f',exact),'  error = ',sprintf('%10.5e',E(i)-exact)])
end



 nn=1000;
 nstates=5;
 rmax
 rtab=linspace(0,rmax,nn);
 psi=zeros(nn-2,nstates);

 
 % ineffective way 
 % sums over all B_i(r), although most are zero
 %figure(1)
 %for state=1:nstates
 %    for ir = 1:nn
 %        r = rtab(ir);        
 %        tmp = 0;
 %        for i=2:n-2
 %            tmp = tmp  + c(i,state)*Bi(rtab(ir),i,t,n,k); 
 %        end
 %        psi(ir,state) = tmp;
 %    end
 %    
 %     plot(rtab,psi(:,state),'Color',rand(1,3));hold on;
 %end
 %figure(2)
 
 % Effective way
 % only k terms in the inner sum
 % all B-splines at once for rtab points
 [B,Bp,Bpp,left]=Bspline(rtab,t,n,k); 
 left = left-k+1;
 leg = [];

 for state=1:nstates
     for ir = 1:nn
         r = rtab(ir);        
         tmp = 0;
         for i=left(ir):left(ir)+k-1  % only non-zero B's
             if(i>n-2) continue;end
             ii = i-left(ir)+1;
             tmp = tmp  + c(i,state)*B(ir,ii);
         end
         psi(ir,state) = tmp;
     end
     plot(rtab,psi(:,state),'Color',rand(1,3));hold on;
     leg = [leg; num2str(state) '. state'];
 end
 legend(leg)



    

