%% Hermitsches Interpolationpolynome
%  Intervall von [0, 1] vgl. HerPoly

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

Or=input('Anzahl der Knoten  ')
Or = Or-1;
%Or = 2;
dx = 1/Or;
lauf = 0:(2*Or+1);
A = zeros(2*(Or+1));
for n=0:Or
    A(2*n+1,:) = (n*dx).^lauf;
    A(2*n+2,:) = lauf.*(n*dx).^(lauf-1);
end
format rat
A(isnan(A))=0;
%C = eye(2*(Or+1))/A;
C = inv(A)
C(abs(C) < 1e-12) = 0;
%% Polynomordnung
C = flipud(C);    % Spalten entsprechen MATLAB Polynom 
%% Visualisierung
x = linspace(0,1);
figure
subplot(3,1,1), hold on
for k=1:length(C)/2
    plot(x,polyval(C(:,2*k-1),x))               % ungerade
end
grid on
subplot(3,1,2), hold on
for k=1:length(C)/2
    plot(x,polyval(C(:,2*k),x))                 % gerade
end
grid on
subplot(3,1,3), hold on
for k=1:length(C)/2                             % Ableitungen gerade
    plot(x,polyval(polyder(C(:,2*k)),x))
end
grid on
%% Berechnung der Integrale
I0 = zeros(2*(Or+1));
I1 = I0;
I2 = I0;
Id0 = I0;
Id1 = I0;
Id2 = I0;
for kp1=1:2*(Or+1)
    for kp2=1:kp1
        p1 = C(kp1,:);
        p2 = C(kp2,:);
        %
        pi0 = conv(p1,p2);
        pi1 = conv([1,0],pi0);
        pi2 = conv([1,0],pi1);
        %
        pd1 = polyder(p1);
        pd2 = polyder(p2);
        %
        pid0 = conv(pd1,pd2);
        pid1 = conv([1,0],pid0);
        pid2 = conv([1,0],pid1);
        % Integrale
        I0(kp1,kp2) = diff(polyval(polyint(pi0),[0,1]));
        I0(kp2,kp1) = I0(kp1,kp2);
        I1(kp1,kp2) = diff(polyval(polyint(pi1),[0,1]));
        I1(kp2,kp1) = I1(kp1,kp2);
        I2(kp1,kp2) = diff(polyval(polyint(pi2),[0,1]));
        I2(kp2,kp1) = I2(kp1,kp2);
        %
        Id0(kp1,kp2) = diff(polyval(polyint(pid0),[0,1]));
        Id0(kp2,kp1) = Id0(kp1,kp2);
        Id1(kp1,kp2) = diff(polyval(polyint(pid1),[0,1]));
        Id1(kp2,kp1) = Id1(kp1,kp2);
        Id2(kp1,kp2) = diff(polyval(polyint(pid2),[0,1]));
        Id2(kp2,kp1) = Id2(kp1,kp2);
    end
end
%%


