function z4_12()
%
% Kopf für jeden Übungszettel
%
ZETTEL = '4';
AUFGABE = '12';

fprintf('\n\n###\tNumerik I SS 2010\t\t###\n###\tZettel %s, Aufgabe %s', ZETTEL, AUFGABE);
fprintf('\t\t###\n###\t\t\t\t\t###\n###\tvon\tFelix Nyenhuis\t\t###\n');
fprintf('###\t\tLouisa Marie Puls\t###\n###\t\tReidar Janssen\t\t###\n\n');
%
%
%


    global BKM;
    BKM = pascal(1001); % _B_inomial_K_oeffizienten_M_atrix
                        % global ausgelagert wg. Performance

    approxBezierAndPlot(10, 1);
    approxBezierAndPlot(100, 2);
    approxBezierAndPlot(1000, 3);

    clear BKM;
    clear all;


function y = f(x)
    y = 1/10 * ( sin(32*pi*x) * ((8*pi*x - 4*pi)^2) );

    
function approxBezierAndPlot(m, fig)
    % Fuehrt die Aufgabenstellung aus fuer m und
    % fig gibt nur eine interne Zahl fuers Plotten an
    fprintf('Berechnung fuer m = %d\n', m);
    TIME = clock;
    
    b = zeros(m + 1, 2); % Datenpaare b = (t, s)

    for i = 0:1:m
        b(i + 1, 1) = i / m; % Das ist t_i
        b(i + 1, 2) = f( b(i + 1, 1) ); % Das ist s_i
    end;
 
    u = 0:0.01:1;
    f_u = zeros(length(u), 1);
    bezier = zeros(length(u), 2);
    fehler = zeros(length(u), 1);
    for i = 1:1:length(u)
        f_u(i) = f(u(i));
        bezier(i, :) = calcBezier(u(i), b);
        fehler(i) = abs( f_u(i) - bezier(i, 2) );
    end;
    
    figure(fig);
    subplot(2, 1, 1);
    %plot(u, f_u, 'm-', bezier(:, 1), bezier(:, 2), 'b-');
    plot(u, f_u, 'm-', u, bezier(:, 2), 'b-');
    xlabel('x'); ylabel('y');
    title(sprintf('f(t) und b(t) fuer m = %d', m));
    legend('f(t)', 'b(t)');
    
    subplot(2, 1, 2);
    plot(u, fehler, 'r-');
    xlabel('x'); ylabel('Fehler');
    title('Fehler zwischen f und b');
    
    ETIME = etime(clock, TIME);
    fprintf('Fertig mit m = %d. Gebrauchte Zeit: %2.2fs\n', m, ETIME);

function [y] = calcBezier(t, b)
    % Wertet die Bezier-Kurve mit Hilfe von Bernsteinpolynomen
    % an der Stelle t aus, wobei b die Kontrollpunkte der Form
    % b_i = (t_i, s_i) mit t_i, s_i reell ist.
    % vgl. Skript Formel (4.8)
    y = zeros(1, 2);
    k = length(b) - 1;
    
    for i = 0:1:k
        y = y + (b(i + 1, :) * BernsteinPolynom(t, k, i + 1));
    end;
    
    
function y = BernsteinPolynom(t, k, i)
    % Wertet das Bernsteinpolynom B_i^k an der Stelle t aus
    % vgl. Skript Formel (4.2)
    y = 0;
    if(~(0 <= i && i <= k)) return; end;
    
    y = binomialKoeffEffizient(k, i) * (t^i) * ((1-t)^(k-i));

    
% Folgende Funktion ist nicht notwendig in diesem Programm
function b = binomialKoeff(n, k)
    % Binomialkoeffizient n ueber k fuer grosse Zahlen
    bK = pascal(max(n - k + 1, k + 1));
    b = bK(n - k + 1, k + 1);

    
function b = binomialKoeffEffizient(n, k)
    % Binomialkoeffizient n ueber k fuer grosse Zahlen
    % in effizienter Art und Weise.
    % Dazu wird das globale Pascalsche Dreieck BKM genutzt
    % Falls diese zu klein oder nicht verfuegbar ist, 
    % wird normal ueber die Funktion binomialKoeff() bestimmt
    global BKM;
    
    if(max(n - k + 1, k + 1) > length(BKM(:, 1)))
        b = binomialKoeff(n, k);
        return;
    end;
    b = BKM(n - k + 1, k + 1);
    
