function IP = InterPol(fg,kn)
% Berechnet die Interpolationspolynome zum 
% Freiheitsgrad fg 1,2 oder 3
% Knotenzahl kn >= 2
% Aufruf IP = InterPol(2,5)
% fg = 1: Lagrange, 2: Hermite, 3: Erweiteter Hermite

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


if nargin~=2  
   error('Zahl der Argumente muss 2 sein')
end
knt = rem(kn,1);
if knt ~=0 || kn<2
    error('Unzulaessiger Knotenwert')
end

switch fg
    case 1
        C = LagPoly(kn);
    case 2
        C = HermPoly(kn);
    case 3
        C = ErwHermPoly(kn);
    otherwise
        error('fg unerlaubter Wert, muss 1, 2 oder 3 sein')
end

IP.Freiheitsgrad = fg;
IP.Knotenzahl = kn;
IP.Polynom = C;
IP.PolyWie = 'C Spaltenweise auswerten: C(:,k) k-tes Polynom';


function C = LagPoly(Or)
Intervall = [0,1];   % Intervall vorgeben
dx = (Intervall(2)-Intervall(1))/Or;
A = zeros(Or+1);
lauf= 0:Or;
for n = lauf
    A(n+1,:)=(Intervall(1)+n*dx).^lauf;
end
C = eye(Or+1)/A;
% Polynomordnung
C = flipud(C);   % Spaltenweise

function C = HermPoly(Or)
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 

function C = ErwHermPoly(Or)

Or = Or-1;
%Or = 2;
dx = 1/Or;
lauf = 0:(3*(Or+1)-1);
A = zeros(3*(Or+1));
for n=0:Or
    A(3*n+1,:) = (n*dx).^lauf;
    A(3*n+2,:) = lauf.*(n*dx).^(lauf-1);
    A(3*n+3,:) = (lauf-1).*lauf.*(n*dx).^(lauf-2);
end
A(isnan(A))=0;
%format rat
C = eye(3*(Or+1))/A;
% Polynomordnung
C = flipud(C);