%% FE Berechnung der Eigenwerte des Wasserstoffatoms
%  Verwendung von Lagrange-Interpolationspolynomen
%  Hier 6 Knoten ( kann in Zeile 14 geaendert werden)
%  1000 FE Elemente s. Zeile 70
%  Integrationsgrenze rmax = 100000 s. Zeile 72

% Vers. 3.3 09.2016 (Vers. 1 03/2000) 
% Copyright W.Schweizer
% Simulation physikalischer Systeme - Computational Physics mit MATLAB 
% 2016

clear, close all, clc
%% Berechnung der Lagrangeschen Interpoaltionspolynome basierend auf [0,1] 
n=6;% 4                  Zahl der Knoten
x = linspace(0,1,n);   % ?quidistant auf [0,1] (prinzipiell freie Wahl)  
xx = repmat(x,n-1,1)';
xy = ones(n,1);
xx = [xy,xx];
%%
xc = cumprod(xx,2);    % Berechnung der Matrix A zu A*C=1 
%format rat
C = inv(xc);
%format short
%% Lagrangepolynom
% Darstellung der ersten drei Polynome
xv = linspace(0,1);
Cli = fliplr(C');  % Cli(n,:) n-tes Lagrangepolynom
plot(xv,polyval(Cli(1,:),xv),xv,polyval(Cli(2,:),xv),xv,polyval(Cli(3,:),xv)),legend('0','1','2')
grid on, shg
%
%% Ableitungen
% Da es sich um Polynome handelt, ist die Berechnung der Ableitungen einfach 
abl = ones(size(Cli));
abl(:,1)=0;
Cliabl=Cli.*fliplr(cumsum(abl,2));
Cliabl(:,n)=[];
%Cli
%% Hilfsintegrale
for ni=1:n
    for nj=1:ni
        zwi = polyint(conv(Cli(ni,:),Cli(nj,:)));                 % x^0*
        Clint0(ni,nj) = polyval(zwi,1)-polyval(zwi,0);
        Clint0(nj,ni) = Clint0(ni,nj);
        %
        zwi = polyint(conv([1,0],conv(Cli(ni,:),Cli(nj,:))));     % x^1*
        Clint1(ni,nj) = polyval(zwi,1)-polyval(zwi,0);
        Clint1(nj,ni) = Clint1(ni,nj);
        %
        zwi = polyint(conv([1,0,0],conv(Cli(ni,:),Cli(nj,:))));   % x^2*
        Clint2(ni,nj) = polyval(zwi,1)-polyval(zwi,0);
        Clint2(nj,ni) = Clint2(ni,nj);
        %
        %                                                                Ableitungen
        zwi = polyint(conv(Cliabl(ni,:),Cliabl(nj,:)));                  % x^0
        Clabint0(ni,nj) = polyval(zwi,1)-polyval(zwi,0);
        Clabint0(nj,ni) = Clabint0(ni,nj);
        %
        zwi = polyint(conv([1,0],conv(Cliabl(ni,:),Cliabl(nj,:))));      % x^1
        Clabint1(ni,nj) = polyval(zwi,1)-polyval(zwi,0);
        Clabint1(nj,ni) = Clabint1(ni,nj);
        %
        zwi = polyint(conv([1,0,0],conv(Cliabl(ni,:),Cliabl(nj,:))));    % x^2
        Clabint2(ni,nj) = polyval(zwi,1)-polyval(zwi,0);
        Clabint2(nj,ni) = Clabint2(ni,nj);
        %
    end
end
%% FE
% quadratische Aufweitung
m = 1000; %3000              % Anzahl der finiten Elemente
ml = 1:m;                    % Laufvariable FE
rmax = 100000;%0;%100;              % Integrationsgrenze (rechts) 
h0 = rmax/(m-1)^2;           % Schrittweite
h = (2*ml-1)*h0;
r0 = (ml-1).^2*h0;
%% Berechnung der Integrale
% Vgl. Gleichungen (4.33 ff)
%
S = h(1)*(r0(1)^2*Clint0 + 2*h(1)*r0(1)*Clint1 + h(1)^2*Clint2);
S = S(:);
%
T = (r0(1)^2*Clabint0 + 2*h(1)*r0(1)*Clabint1 + h(1)^2*Clabint2)/h(1);
T = T(:);
%
V = -2*h(1)*(r0(1)*Clint0 + h(1)*Clint1);
V = V(:);
%
for k=2:m
    Sneu = h(k)*(r0(k)^2*Clint0 + 2*h(k)*r0(k)*Clint1 + h(k)^2*Clint2);
    Sneu = Sneu(:);
    S(end) = S(end)+Sneu(1);
    Sneu(1) = [];
    S = [S;Sneu];
    %
    Tneu = (r0(k)^2*Clabint0 + 2*h(k)*r0(k)*Clabint1 + h(k)^2*Clabint2)/h(k);
    Tneu = Tneu(:);
    T(end) = T(end)+Tneu(1);
    Tneu(1) = [];
    T = [T;Tneu];
    %
    Vneu = -2*h(k)*(r0(k)*Clint0 + h(k)*Clint1);
    Vneu = Vneu(:);
    V(end) = V(end)+Vneu(1);
    Vneu(1) = [];
    V = [V;Vneu];
end
%% Eigenwertgleichung
H = T + V;
%% Indexgymnastik
%n-1           % Ordnung Interpolationspolynom
p = 0;         % Interpolationstyp  p=fg
%m             % Zahl der FE
blockdim = n*(p+1);
indlin = 1:blockdim^2;
[iz,is]=ind2sub([blockdim,blockdim],indlin);
izneu = iz;
isneu = is;
for k= 2:m  
    izneu = izneu+(blockdim-1);%2;   % 2 blockdim-1 f??r Lagrange  
    isneu = isneu+(blockdim-1);
    switch p
        case 0
            iz = [iz(1:(end-(p+1))),izneu];  % -1 Abh??ngig von ??berlappung p+1
            is = [is(1:(end-(p+1))),isneu];
    end
end
%
S = sparse(iz,is,S);
H = sparse(iz,is,H);
%figure, spy(S),shg
figure, spy(H), shg
%% Eigenwertproblem loesen
E = eigs(H,S,25,-1);%'sa')%100
E = E/2
    
    
 

