%  function Multivar( a, F )
%
%    Analyse der beiden Kräfte F1 und F2, 
%    die auf einen Träger an den Positionen a einwirken,
%    mittels Multivariater Regression.
%    Die Funktion DurchbiegKoeff berechnet die 
%    Durchbiegung w an den Punkten z1.
%    Diese Werte werden um 10% verrauscht.
%    Mit den verrauschten w-Werten werden am Ende 
%    die Kräfte rekonstruiert
%
%    Beispielaufruf: Multivar( [0.3;0.8], [15;45] )
%      a = [a1;a2] : relative Position zwischen 0 und 1
%      F = [F1;F2] : vorgegebene Kräfte in Newton
%
%    Verwendetes Modell für den Träger:
%      Länge: 1 m
%      Flächentr.-Moment: I = 0.1e-9 m^4 (= 100 mm^4)
%      E-Modul: E = 0.2e12 N/m^2 (= 200 000 N/mm^2)
%
function Multivar( a, F )
    
    % rel. Messpunkte z1 zwischen 0 und 1
    d  = 0.01;          % Abstand der Messpunkte
    z1 = 0.0:d:1.0;
    anz = length( z1 ); 
    
    % Berechnung der beiden Vektoren mit den Durchbiege-Koeff.
    A1 = zeros(anz,1);
    A2 = zeros(anz,1);
    for( k = 1:anz )
       A1(k,1) = DurchbiegKoeff( z1(k), a(1) );
       A2(k,1) = DurchbiegKoeff( z1(k), a(2) );
    end
    
    % Matrix für beide Positionen 
    A = [ A1, A2 ];
        
    % Berechnung der einzelnen Biegelinien
    w1 = A1*F(1);
    w2 = A2*F(2);
    w  = A*F;
    
    % Daten plotten
    plotData( z1, w1, w2, w, F, a, d );
      
    % w-Daten verrauschen
    mult = 0.1;      % 10% Variation
    wr = zeros(anz,1);
    for( k=1:anz )
       wr(k,1) = w(k,1) * ( 1 + (rand-0.5)*mult );        
    end
    
    % Multivariate Regression
    FM = (A'*A)^-1 * A' * wr;
    
    % Ausgabe des Ergebnisses
    fprintf( 'Rekonstruierte Kräfte: \nFM1 = %.2f N, FM2 = %.2f N\n', ...
             FM(1), FM(2) );
    fprintf( 'Ursprüngliche Kräfte: \n F1 = %.2f N,  F2 = %.2f N\n', ...
             F(1), F(2) );
    
end


% Testfunktion für DurchbiegKoeff
function Test()
    % F = Kraft (m*g), m = 2 kg
    F = 2*9.81;
    A = DurchbiegKoeff( 0.5, 0.5 );
    % Durchbiegung sollte ca. -2 cm = -0.02 m sein
    w = A*F;
    fprintf( 'Test von DurchbiegKoeff: %g sollte ca. -0.02 sein \n', w );
end


function plotData( z1, w1, w2, w, F, a, d )
    % plot( z1, w1, 'b--', z1, w2, 'k-.', z1, w, 'r-' )
    % SW-Plot: 
    plot( z1, w1, 'k--', z1, w2, 'k-.', z1, w, 'k-' )
    
    xlabel( 'Mess-Position / m' );
    ylabel( 'Durchbiegung / m' );
    t = sprintf( 'Durchbiegung, gemessen im Abstand %.2f m', d );
    title( t )
    t1 = sprintf( 'F1 = %.0f N, a = %.1fm', F(1), a(1) );
    t2 = sprintf( 'F2 = %.0f N, a = %.1fm', F(2), a(2) );
    legend( t1, t2, 'F1 + F2' );
    
    % Pfeile anden Positionen der Kräfte
    w_a1 = DurchbiegKoeff(a(1),a(1))* F(1);
    w_a2 = DurchbiegKoeff(a(2),a(2))* F(2);

    % relative Positionen de Axes:  x,y,width,heigth
    pos = get( gca, 'Position' );
    xpos = pos(1);  ypos = pos(2);
    wpos = pos(3);  hpos = pos(4);
    YLim = get( gca, 'YLim' );
    yMin = YLim(1); 
    
    % Skalierung
    a1x = a(1)*wpos + xpos;
    a2x = a(2)*wpos + xpos;
    a1y = (1-w_a1/yMin)*hpos + ypos;
    a1y0 = a1y + hpos/6; 
    a2y = (1-w_a2/yMin)*hpos + ypos;
    a2y0 = a2y + hpos/6; 
    
    annotation( 'arrow', [a1x,a1x], [a1y0,a1y], ...
         'HeadStyle', 'cback2' );
    annotation( 'arrow', [a2x,a2x], [a2y0,a2y], ...
         'HeadStyle', 'cback2' );  
end


% function A = DurchbiegKoeff( z1, a )
%   berechnet die Durchbiegung
%     z1 = rel. Punkt zur Berechnung, vom linken Ende, zwischen [0,1]
%     a  = rel. Abstand der Kraft, vom linken Ende [0,1]
function A = DurchbiegKoeff( z1, a )

    L = 1;         % 1 m langer Träger
    E = 0.2e12;    % N/m^2 E-Modul 200 000 N/mm^2
    I = 0.1e-9;    % m^4 Flächenträgheitsmoment (100 mm^4)
    % I = 20e-9;   % m^4 Flächenträgheitsmoment (20 000 mm^4)
    
    z2 = z1 - a;   % rel. Abstand des Punktes von der Kraft
    b  = 1 - a;    % rel. Abstand der Kraft zum rechten Ende
    
    % Punkt z1 liegt links von der Krafteinwirkung a
    if( z1 <= a )
       A = - L^3*a*b^2/(6*E*I) * ...
           ( (1+L/b)*z1/L - z1^3/(a*b*L) ); 
    % Punkt liegt rechts von der Krafteinwirkung   
    else
       A = - L^3*a^2*b/(6*E*I) * ...         
           ( (1+L/a)*(b-z2)/L - (b-z2)^3/(a*b*L) ); 
    end

end