% Kreisringplatte, konstante Flaechenlast, am Innenrand eingespannt
clear all

ri = 1000  ; 
ra = 3000  ;
t  =   30  ;
E  = 2.1e5 ;
ny = 0.3   ;
p0 = 0.004 ;
K  = E*t^3/(12*(1-ny^2)) ;

nA   = 1000    ;       % Anzahl der Abschnitte
n    = nA + 5  ;       % Anzahl der Gleichungen
h    = (ra-ri)/nA ;    % Schrittweite

A  = zeros (n,7) ;     % Rechteckmatrix fuer Aufnahme der Bandmatrix
b  = zeros (n,1) ;     % Nullvektor ("rechte Seite")
pi = zeros (n,1) ;     % Flaechenlastintensitaeten an Stuetzpunkten
r  = zeros (n,1) ;     % Radien bis zu Stuetzpunkte


iLinks  = 3 ;                             % Nummern der Randpunkte
iRechts = n - 2 ;

pi(1:n) = p0 ;                            % Belastung
r(iLinks:iRechts) = ri:h:ra ;             % Radien

%Randbedingungen:
A(1:2,:)   = [0  0  0   0  0  1  0 ;
              0  0  0  -1  0  1  0 ] ;    % Starr eingespannter Rand links und freier Rand rechts
A(n-1:n,:) = [    -0.5      1+h/ra+h^2/(2*ra^2)   -2*h/ra     -1+h/ra-h^2/(2*ra^2)  0.5  0  0 ;
              1-ny*h/(2*ra)        -2           1+ny*h/(2*ra)           0            0   0  0 ] ; 

for i = 3:n-2    % Matrix A:
   hdr = h/r(i) ;
   A(i,:)= [0  1-hdr  -4+2*hdr-hdr^2-0.5*hdr^3  6+2*hdr^2  -4-2*hdr-hdr^2+0.5*hdr^3  1+hdr  0] ;    
   b(i) = pi(i)*h^4/K ;                                        % Standardgleichungen
end 

w = gabamp (A , b) ;                                           % Berechnung der Durchsenkung 

clf;                                                           % Graphische Ausgabe der Biegelinie:
subplot (3,1,1) ; plot (r(3:n-2) , w(3:n-2)) , axis ij , grid on , title ('w(r)') , xlim([0 ra]) 

for i=3:n-2                      
   hdr = h/r(i) ;
   mr(i) = -K*((1-0.5*ny*hdr)*w(i-1) - 2*w(i) + (1+0.5*ny*hdr)*w(i+1))/h^2 ;  % Biegemoment
   qr(i) = -K*(-0.5*w(i-2) + (1+hdr+0.5*hdr^2)*w(i-1) - 2*hdr*w(i) + ...
                            (-1+hdr-0.5*hdr^2)*w(i+1) + 0.5*w(i+2))/h^3 ;     % Querkraft
end
   
subplot (3,1,2) ; plot (r(3:n-2) , mr(3:n-2)) , grid on , title ('m_r(r)') , xlim([0 ra]) 
subplot (3,1,3) ; plot (r(3:n-2) , qr(3:n-2)) , grid on , title ('q_r(r)')   , xlim([0 ra])   

disp(['Maximale Durchbiegung: w_max  = ', num2str(max(w))]) ;
disp(['Maximales Biegemoment: mr_max = ', num2str(max(abs(mr)))]) ;
disp(['Maximale Querkraft:    qr_max = ', num2str(max(abs(qr)))]) ;   