%% FE Berechnung der Eigenwerte des Yukawapotentials
%  Verwendung von Lagrange-Interpolationspolynomen
%  Hier 3 Knoten (kann in Zeile 15 geaendert werden)
%  FE Elemente s. Zeile 72ff

% Vers. 1.4 08.2016 
% Copyright W.Schweizer
% Simulation physikalischer Systeme - Computational Physics mit MATLAB 
% 2016

%%
clear, close all, clc
%% Berechnung der Lagrangeschen Interpolationspolynome basierend auf [0,1] 
n=3;%6;% 4                  Zahl der Knoten
x = linspace(0,1,n);   % aequidistant auf [0,1] (prinzipiell freie Wahl)  
xx = repmat(x,n-1,1)';
xy = ones(n,1);
xx = [xy,xx];
%% Yukawa-Parameter
rya = 0.01;%1.0;
V0ya = 1;
%%
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 + aequidistante Aufweitung
m = 500;%3000;%500;%1000;%1250;%1500; %3000              % Anzahl der finiten Elemente
m1= 70;%1000;%70;                                        % Anzahl fuer quadratische Aufweitung
ml = 1:m1;                    % Laufvariable FE
rmax = 100;%0;%100;           % Quadratische Grenze (rechts) 
h0 = rmax/(m1-1)^2;           % Schrittweite Start
h = (2*ml-1)*h0;
h = [h,h(end)*ones(1,m-m1)];
%r0 = (ml-1).^2*h0;
r0 = [0,cumsum(h)];
r0(end) = [];
%% Berechnung der Integrale
% Vgl. Gleichungen
%
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(:);
%
% Hilfsintegrale fuer das Potential

nlauf = 2*n;
alpha = -rya*h(1);
cint(nlauf) = 1/alpha*(exp(alpha)-1);
cvor = 1/alpha*exp(alpha);
for k=1:nlauf-1
    cint(nlauf-k) = cvor - k/alpha*cint(nlauf-k+1);
end
cint0 = cint(2:end);
for ni=1:n
    for nj=1:ni
        zwi=(conv(Cli(ni,:),Cli(nj,:)));                   % x^0*
        ClintV0(ni,nj)= zwi*cint0';    
        ClintV0(nj,ni) = ClintV0(ni,nj);
        %
        zwi = (conv([1,0],conv(Cli(ni,:),Cli(nj,:))));     % x^1*
        ClintV1(ni,nj) = zwi*cint';
        ClintV1(nj,ni) = ClintV1(ni,nj);
    end
end
V = -2*h(1)*V0ya*exp(alpha)*(r0(1)*ClintV0 + h(1)*ClintV1);
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];
    %
    
    nlauf = 2*n;
    alpha = -rya*h(k);
    cint(nlauf) = 1/alpha*(exp(alpha)-1);
    cvor = 1/alpha*exp(alpha);
        for kv=1:nlauf-1
            cint(nlauf-kv) = cvor - kv/alpha*cint(nlauf-kv+1);
        end
    cint0 = cint(2:end);
    
    for ni=1:n
    for nj=1:ni
        zwi=(conv(Cli(ni,:),Cli(nj,:)));                   % x^0*
        ClintV0(ni,nj)= zwi*cint0';    
        ClintV0(nj,ni) = ClintV0(ni,nj);
        %
        zwi = (conv([1,0],conv(Cli(ni,:),Cli(nj,:))));     % x^1*
        ClintV1(ni,nj) = zwi*cint';
        ClintV1(nj,ni) = ClintV1(ni,nj);
    end
    end
    Vneu = -2*h(k)*V0ya*exp(alpha)*(r0(k)*ClintV0 + h(k)*ClintV1);
    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-1
%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 fuer Lagrange  
    isneu = isneu+(blockdim-1);
    switch p
        case 0
            iz = [iz(1:(end-(p+1))),izneu];  % -1 Abhaengig von Ueberlappung p+1
            is = [is(1:(end-(p+1))),isneu];
    end
end
%
S = sparse(iz,is,S);
H = sparse(iz,is,H);
figure, spy(H), shg
%% Eigenwertproblem l?sen
E = eigs(H,S,25,-3);%'sm')%-1)%'sa')%100
E = E/2