% Ritzsches Verfahren: Biegeschwingungen gerader Traeger mit Polynom-Ansatzfunktionen,
% Aufgabe 32-11.

clear all ; close all ;

% 11111111111 ANPASSEN AN DAS AKTUELLE PROBLEM VON HIER ... 111111111111111111111111111
% Parameter:
L1    =    0.4 ;
L2    =    0.3 ;
L3    =    0.5 ;
L4    =    0.6 ;
EI1   =   3000 ;
EI2   =   4000 ;
rhoA1 =    2.4 ;
rhoA2 =    3.2 ;
m1    =      8 ;
m2    =     10 ;
c1    = 120000 ;
c2    =  70000 ;
cT    =   6000 ;
Lges  = L1+L2+L3+L4 ;

% Ansatzfunktionen und deren 1. und 2. Ableitungen:

m  = 5 ; tl = L1+L2+L3 ;          % Anzahl der Ansatzfuntionen (wenn m < 5 gesetzt ...
P1 = [1/tl^3 -1/tl^2 0 0] ;     % ... wird, werden nur die ersten ...
P2 = [1/tl^4 -1/tl^3 0 0 0] ;   % ... m Ansatzfunktionen verwendet)
P3 = [1/tl^5 -1/tl^4 0 0 0 0] ;
P4 = [1/tl^6 -1/tl^5 0 0 0 0 0] ;
P5 = [1/tl^7 -1/tl^6 0 0 0 0 0 0] ;
% 11111111111 ... BIS HIER 111111111111111111111111111111111111111111111111111111111111

P1D = polyder (P1);
P2D = polyder (P2);
P3D = polyder (P3);
P4D = polyder (P4);
P5D = polyder (P5);

P1DD = polyder (P1D) ;
P2DD = polyder (P2D) ;
P3DD = polyder (P3D) ;
P4DD = polyder (P4D) ;
P5DD = polyder (P5D) ;

% ------------------------------------------------------------------------------------
% Vorbereitung der numerischen Integration und der Ergebnisauswertung:
n  = 180 ;                % Anzahl der Abschnitte für numerische Integration
dz = Lges / n ;           % Breite eines Integrationsintervalls
nS = n + 1 ;              % Anzahl der Stützstellen
zS = 0 : dz : Lges ;      % Koordinaten der Stützpunkte

EIS   = zeros (nS , 1) ;
rhoAS = zeros (nS , 1) ;

% 22222222222222 ANPASSEN AN DAS AKTUELLE PROBLEM VON HIER ... 222222222222222222222222
ic1 = round (n*L1/Lges+1) ;
im1 = round (n*(L1+L2)/Lges+1) ;
icT = round (n*(L1+L2+L3)/Lges+1) ;
im2 = n+1 ;

EIS(1  :im1)   = EI1   ;     
EIS(im1:im2)   = EI2   ;     
rhoAS(  1:im1) = rhoA1 ;     
rhoAS(im1:im2) = rhoA2 ;     
% 22222222222222 ... BIS HIER 222222222222222222222222222222222222222222222222222222222

% Funktionswerte und 1. und 2. Ableitung der Ansatzfunktionen an den Stützstellen
vS  (:,1) = polyval (P1   , zS)' ;
vS  (:,2) = polyval (P2   , zS)' ;
vS  (:,3) = polyval (P3   , zS)' ;
vS  (:,4) = polyval (P4   , zS)' ;
vS  (:,5) = polyval (P5   , zS)' ;
vdS (:,1) = polyval (P1D  , zS)' ;
vdS (:,2) = polyval (P2D  , zS)' ;
vdS (:,3) = polyval (P3D  , zS)' ;
vdS (:,4) = polyval (P4D  , zS)' ;
vdS (:,5) = polyval (P5D  , zS)' ;
vddS(:,1) = polyval (P1DD , zS)' ;
vddS(:,2) = polyval (P2DD , zS)' ;
vddS(:,3) = polyval (P3DD , zS)' ;
vddS(:,4) = polyval (P4DD , zS)' ;
vddS(:,5) = polyval (P5DD , zS)' ;

plot(zS , vS(:,1:m)) ; grid ; title('Ansatzfunktionen:') ;

% ------------------------------------------------------------------------------------
% Aufbau des allgemeinen Matrizen-Eigenwertproblems:
K = zeros (m , m) ;
M = zeros (m , m) ;

for ii = 1:m                                          % Schleife über alle Gleichungen
  for jj = 1:m                                        % ii-te Zeile (Koeffizienten kij)
    Summe1 = EIS  (1)*vddS(1,ii)*vddS(1,jj) + EIS  (nS)*vddS(nS,ii)*vddS(nS,jj) ;
    Summe2 = rhoAS(1)*vS  (1,ii)*vS  (1,jj) + rhoAS(nS)*vS  (nS,ii)*vS  (nS,jj) ;
    faktor = 4 ;                                      % Numerische Integration ...
    for k = 2:n                                       % ... nach Simpsonscher Regel ... 
      Summe1 = Summe1 + EIS  (k) * vddS(k,ii) * vddS (k,jj) * faktor ;
      Summe2 = Summe2 + rhoAS(k) * vS  (k,ii) * vS   (k,jj) * faktor ;
      if   (faktor == 4) faktor = 2 ;                 
      else               faktor = 4 ;                 
      end ;
   end
% 33333333333333 ANPASSEN AN DAS AKTUELLE PROBLEM VON HIER ... 3333333333333333333333333
    K(ii,jj) = Summe1*dz/3 + 2*c1*vS(ic1,ii)*vS(ic1,jj) + 2*c2*vS(im2,ii)*vS(im2,jj) ...
                           + cT*vdS(icT,ii)*vdS(icT,jj); 
    M(ii,jj) = Summe2*dz/3 + m1*vS(im1,ii)*vS(im1,jj) + m2*vS(im2,ii)*vS(im2,jj) ;
% 33333333333333 ... BIS HIER ... 333333333333333333333333333333333333333333333333333333
  end   
end

% Lösen des allgemeinen Eigenwertproblems:
[V D] = eig (K , M) ;

for i=1:m
   f(i) = sqrt(D(i,i)) / (2*pi) ;
end
[f , Index] = sort(f) ;
if (m > 3) mm = 3 ; else mm = m ; end  ;  % Maximal 3 Eigenfrequenzen ausgeben ...
Eigenfrequenzen = f(1:mm)

figure ;
z = 0 : Lges/n : Lges ;
for i = 1:mm                              % ... bzw. Eigenschwingungsformen zeichnen
    v = vS(:,1:m) * V(:,Index(i)) ;
    subplot (mm,1,i) ; plot (z , v) , grid on , ... 
             title (strcat(num2str(i),'. Eigenschwingungsform:')) ;
end