% Anhang B: Gerader mehrfeldriger Traeger
% Differenzenverfahren

clear all

L  = 2100  ;     % Länge
lb = 800   ;
q0 = 1     ;     % Linienlast
F1 = 2000  ;
F2 = 1000  ;
EI = 4.2e9 ;     % Biegesteifigkeit  

nA = 210    ;    % Anzahl der Abschnitte
k  = nA + 5 ;    % Knotenanzahl
h  = L / nA ;    % Schrittweite

A = zeros (k,k) ;      % Matrix A
b = zeros (k,1) ;      % Rechte Seite b

A(1:2,1:5)     = [0  0 1 0 0 ; ...
                  0 -1 0 1 0];         % Randbedingung: Einspannung links
A(k-1:k,k-4:k) = [0  1 -2 1  0 ; ...
                  1 -2  0 2 -1];     
b(k)           = 2*F2*h^3/EI   ;       % Randbedingung: Freier Rand rechts mit Einzelkraft

for i=3:k-2                            % Beschreibung der Matrix A in einer Schleife mit
    A(i,i-2:i+2)=[1 -4 6 -4 1];        % den Standardwerten  1 -4 6 -4 1
end

iB = 25 ;     % Punkt B (Lager)
A(iB,iB-2:iB+2) = [0  0 1  0 0] ;

iC = 105 ;    % Punkt C (Lager)
A(iC,iC-2:iC+2) = [0  0 1  0 0] ;

iD = 126 ;    % Punkt D (Gelenk)
A(iD-1,iD-3:iD+1) = [1 -4 5 -2 0] ;
A(iD  ,iD-2:iD+2) = [1 -2 2 -2 1] ;
A(iD+1,iD-1:iD+3) = [0 -2 5 -4 1] ;

iF = 182 ;    % Punkt E (Lager)
A(iF,iF-2:iF+2) = [0  0 1  0 0] ;

for i=iB+1:iC-1                          
    b(i) = q0*(3+4*(i-iB)*h/lb)*h^4/EI ;  % Linienlast
end

iE = 153 ;          % Punkt E
b(iE) = F1*h^3/EI ; % "Verschmierte" Einzelkraft F1 

v =A \ b ;          % Loesung des Gleichungssystems (Berechnung der Durchbiegung v)

clf;
x = 0 : h : L ;
subplot (3,1,1) ; plot (x , v(3:k-2)) , grid on , title ('Durchbiegung')            % Graphische Ausgabe der Biegelinie

for i=1:k                          
    EIi(i) = EI ;  
end

EIi(iD) = 0 ;

for i=3:k-2                       
    M (i) = -EIi(i)*(v(i-1)-2*v(i)+v(i+1))/h^2 ;  % Biegemoment
    FQ(i) = 0.5*(EIi(i-1)*v(i-2)-2*EIi(i-1)*v(i-1)+(EIi(i-1)-EIi(i+1))*v(i)+2*EIi(i+1)*v(i+1)-EIi(i+1)*v(i+2))/h^3 ;  % Querkraft
end
    
subplot (3,1,2) ; plot (x , M (3:k-2)) , grid on , title ('Biegemoment')            % Graphische Ausgabe des Momentenverlaufs
subplot (3,1,3) ; plot (x , FQ(3:k-2)) , grid on , title ('Querkraft')              % Graphische Ausgabe des Querkraftverlaufs
   
AbsenkungGelenk = v(iD)   % Ausgabe in das Command window:    
AbsenkungPunktE = v(iE)   % Zahlenwerte der Verschiebungen ...
AbsenkungPunktG = v(213)  % ... an drei speziellen Punkten ...
MbA = M(3)                % ... und das Einspannmoment am Punkt A
