% Differenzenverfahren, Aufgabe 18-13
clear all

L1  = 200 ;
L2  = 300 ;
L3  = 100 ;
L   = L1+L2+L3 ;
q1  = 2 ;
q2  = 3 ;
E   = 210000 ;
I1  =    500 ;
I2  =    800 ;
I3  =    200 ;
F1  =    200 ;
F2  =    600 ;
M1  =  40000 ;
M2  =  50000 ;
c1  =     40 ;
c2  =     30 ;
cT  = 120000 ;

EIB = E*I1   ;       % Bezugs-Traegheitsmoment
nA  = 3000   ;       % Anzahl der Abschnitte
n   = nA + 5 ;       % Anzahl der Gleichungen
h   = L / nA ;

% Markante Punkte:
iF2 = round (nA*L1/L      + 3) ;
iM2 = round (nA*(L1+L2)/L + 3) ;
iRechts = n - 2 ;

A  = zeros (n,7) ;      % Rechteckmatrix fuer Aufnahme der Bandmatrix
b  = zeros (n,1) ;      % Nullvektor ("rechte Seite")
qi = zeros (n,1) ;      % Linienlastintensitaeten an Stuetzpunkten
mi = zeros (n,1) ;      % EIi = mi * EIB
ki = zeros (n,1) ;      % Bettungszahlen an Stuetzpunkten

% Biegesteifigkeit (mue-Werte):
mi(1:iF2-1) = E*I1/EIB ;
mi(iF2)     = 0 ;                                %Gelenk
mi(iF2+1:iM2-1) = E*I2/EIB ;
mi(iM2)     = (E*I2+E*I3)/(2*EIB) ;
mi(iM2+1:n) = E*I3/EIB ;

qi(iF2:iM2) = q1 : (q2-q1)/(iM2-iF2) : q2 ;      % Linienlast
qi(iF2)     = q1/2 ;
qi(iM2)     = q2/2 ;

qi(iF2)   = qi(iF2)   + F2/h ;                   % Kraft F2 ("verschmiert")
qi(iM2-1) = qi(iM2-1) + M2/(2*h^2) ;             % Moment M2, ersetzt durch Kraeftepaar,
qi(iM2+1) = qi(iM2+1) - M2/(2*h^2) ;             % dieses dann "verschmiert"
           
ki(iF2)   = c2/h ;                               % Feder c2

%Randbedingungen:
A(1:2,:) = [0 0  0 0 mi(3)+cT*h/(2*EIB) -2*mi(3) mi(3)-cT*h/(2*EIB) ;
            0 0  mi(2) -2*mi(2) mi(2)-mi(4)-2*c1*h^3/(EIB) 2*mi(4) -mi(4)] ;
b(1) = -M1*h^2/(EIB)   ;                        % Randbedingungen links: Feder c1, 
b(2) = -2*F1*h^3/(EIB) ;                        % Drehfeder, Kraft F1, Moment M1
                
A(n-1:n,:) = [ 0 0 1 0 0 0 0  ;
              -1 0 1 0 0 0 0] ;                  % Randbedingungen rechts: Einspannung

for i = 3:n-2    % Matrix A:                          
    A(i,:)= [0 mi(i-1) -2*(mi(i-1)+mi(i)) mi(i-1)+4*mi(i)+mi(i+1)+ki(i)*h^4/EIB ...
                       -2*(mi(i)+mi(i+1)) mi(i+1) 0] ;    
    b(i) = qi(i)*h^4/EIB ;                       % Standardgleichungen
end

A(iM2,:)= [0 0 0 1 0 0 0] ;                      % Zwischenstuetze
b(iM2)  = 0 ;

t1 = cputime ;
v = gabamp (A , b) ;                             % Berechnung der Durchsenkung
ZeitGlSyst = cputime - t1 

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

for i=3:n-2                       
    Mb(i) = -mi(i)*EIB * (v(i-1)-2*v(i)+v(i+1))/h^2 ;                 % Biegemoment
%    FQ(i) = 0.5*EIB * (mi(i-1)*v(i-2)-2*mi(i-1)*v(i-1)+(mi(i-1)-mi(i+1))*v(i) ...
%                        +2*mi(i+1)*v(i+1)-mi(i+1)*v(i+2))/h^3 ;       % Querkraft
end
    
subplot (2,1,2) ; plot (z , Mb(3:n-2)) , grid on , title ('Biegemoment') 
%subplot (3,1,3) ; plot (z , FQ(3:n-2)) , grid on , title ('Querkraft') 
   
AbsenkungLinks  = v(3)        % Ausgabe in das Command window:
AbsenkungF2     = v(iF2)      % Zahlenwerte der Verschiebungen ...          
MbEinspann      = Mb(iRechts) % ... und das Einspannmoment rechts