function Biegno
%BIEGNO für Matlab
%
%Normalspannungen bei Biegung mit Normalkraftbeanspruchung beliebiger, aus
%Teilflächen zusammengesetzter Querschnitte
%
%Autor: Gerhard Knappstein
%Matlab-Bearbeitung: Philipp Fleiger
%
%1987 - 2007
%


men_contr=0;    %Menue-Kontrollvariable
beenden=0;

Flaechen=[]; %Daten Teilflächen
Lasten=[];  %Lastdaten [N, M_y, M_z]
Punkte=[];
Gesamtflaeche=[];


%--------------------------------------------------------------------------
%Textmenue
%--------------------------------------------------------------------------

while beenden == 0

    clc
    
    disp('----------------------------------------------------------------- ');
    disp('**************************');
    disp('* B I E G N O für Matlab *');
    disp('**************************');
    disp(' ');
    disp('Universität Siegen    Lehrstuhl für numerische Mechanik ');
    disp('G. Knappstein, P. Fleiger  2007 ');
    disp('----------------------------------------------------------------- ');
    disp(' ');
    disp(' ');

    switch men_contr
        
        case 0
            disp('( 1 ) Erste Teilfläche erzeugen');
            %disp('( 2 ) Datensatz laden');
            disp('( 2 ) Beenden');
            
        case 1
            disp('( 1 ) Weitere Teilfläche hinzufügen');
            disp('( 2 ) Mit Berechung fortfahren (Flächeneingabe abschließen)');
            disp('( 3 ) Ausgabe der aktuellen Konfiguration');
            disp('( 4 ) Beenden');
            
        case 2
            disp('( 1 ) Kräfte und Momente eingeben');
            disp('( 2 ) Punkt hinzufügen, in dem die Spannung berechnet werden soll');
            disp('( 3 ) Eingaben überprüfen');
            disp('( 4 ) Berechnung durchführen');
            disp('( 5 ) Beenden');
            
        case 3
            disp('( 1 ) Kräfte und Momente ändern');
            disp('( 2 ) Punkt hinzufügen, in dem die Spannung berechnet werden soll');
            disp('( 3 ) Eingaben überprüfen');
            disp('( 4 ) Berechnung erneut durchführen');
            disp('( 5 ) Ergebnisse anzeigen');
            disp('( 6 ) Ausgabedatei erzeugen');
            disp('( 7 ) Daten speichern');
            disp('( 8 ) Beenden');
    end  %switch
    
    disp(' ');
    disp(' ');
    Auswahl=input('Bitte wählen Sie eine Option aus:   ');
    
    switch men_contr
        
        case 0
            
            switch Auswahl
                
                case 1
                    clc
                    erste_flaeche;
                    men_contr=1;
                    
%                 case 2
%                     clc
%                     laden;
%                     beenden=0;
%                     men_contr=1;
                    
                case 2
                    beenden=1;
                    
                otherwise
                    disp(' ');
                    disp('Ungültige Eingabe! Bitte versuchen Sie es erneut.')
                    pause
            end
            
        case 1
            
            switch Auswahl
                
                case 1
                    clc
                    naechste_flaeche;
                    
                case 2
                    clc
                    zus_flaeche=flaech_berech (Flaechen);
                    men_contr=2;
                    
                case 3
                    clc
                    ausgeben_fl (Flaechen, Gesamtflaeche);                    
                  
                case 4
                    beenden=1;
                    
                otherwise
                    disp(' ');
                    disp('Ungültige Eingabe! Bitte Taste drücken und erneut versuchen.')
                    pause
            end
            
        case 2
            
            switch Auswahl
                
                case 1
                    clc
                    Lasten=kraefte_eingeben;
                    
                case 2
                    clc
                    Punkte=punkt_eingeben(Punkte);
                    
                case 3
                    clc
                    ausgeben(Lasten, Punkte);
                    
                case 4
                    [sigma, m, n]=berechnung(Lasten, Punkte, zus_flaeche);
                    men_contr=3;
                    
                case 5
                    beenden=1;
                    
                otherwise
                    disp(' ');
                    disp('Ungültige Eingabe! Bitte Taste drücken und erneut versuchen.')
                    pause
            end
                    
        case 3
            
            switch Auswahl
                
                case 1
                    clc
                    Lasten=kraefte_eingeben;
                    
                case 2
                    clc
                    Punkte=punkt_eingeben(Punkte);
                    
                case 3
                    clc
                    ausgeben(Lasten, Punkte);
                    
                case 4
                    sigma=berechnung(Lasten, Punkte, zus_flaeche);
                    
                case 5
                    clc
                    ergeb_aus(sigma, Lasten, Punkte, zus_flaeche, m, n);
                    
                case 6
                    clc
                    datei_aus(sigma, Lasten, Punkte, zus_flaeche);
                    
                case 7
                    clc
                    speichern(sigma, Lasten, Punkte, zus_flaeche);                    
                    
                case 8
                    beenden=1;
                    
                otherwise
                    disp(' ');
                    disp('Ungültige Eingabe! Bitte Taste drücken und erneut versuchen.')
                    pause
            end
    end  %switch men_contr
           



end %Menue


%--------------------------------------------------------------------------
%Unterfunktionen
%--------------------------------------------------------------------------


%Flächen
    function erste_flaeche
        disp(' ');
        disp('Bitte wählen Sie die Art der Teilfläche gemäß der Kennzahl aus Tablle 3!')
        disp(' ');
        Kennzahl=input('Kennzahl =  ');
        Flaechen=[Flaechen; eingabe_flaeche(Kennzahl)];
    end

    function naechste_flaeche
        disp(' ');
        disp('Bitte wählen Sie die Art der zuzufügenden Teilfläche gemäß der Kennzahl aus Tablle 3!')
        disp(' ');
        Kennzahl=input('Kennzahl =  ');
        Flaechen=[Flaechen; eingabe_flaeche(Kennzahl)];
    end

    function Geometrie=eingabe_flaeche (Kennzahl)
        %Eingabe der Geometriedaten
        clc
        disp(' ');
        disp('Bitte geben Sie die Parameterwerte für die neue Fläche ein!')
        disp(' ');
        
        Kz=abs(Kennzahl);
        fehler=0;
        
        switch Kz
            
            case 1
                disp(' ');
                disp('Rechteck:')
                disp(' ');
                b=input('Breite b:  ');
                h=input('Höhe h:  ');                
                y_quer=input('y-Koordinate des Schwerpunktes:  ');
                z_quer=input('z-Koordinate des Schwerpunktes:  ');
                %Berechnung der notwendigen Größen:
                A=b*h;
                I_y=(A*h^2)/12;
                I_z=(A*b^2)/12;
                I_yz=0;
            case 2
                disp(' ');
                disp('Profil oder jede andere Flächenform:')
                disp(' ');
                A=input('Fläche A:  ');                               
                y_quer=input('y-Koordinate des Schwerpunktes:  ');
                z_quer=input('z-Koordinate des Schwerpunktes:  ');
                I_y=input('axiales Flächenträgheitsmoment bezüglich der y-Achse I_y:  ');
                I_z=input('axiales Flächenträgheitsmoment bezüglich der z-Achse I_z:  ');
                I_yz=input('Zentrifugalmoment I_yz (Vorzeichen bitte Tabelle 1 entnehmen):  ');
            case 3
                disp(' ');
                disp('Zwei gleiche Profile oder zwei andere gleiche Flächenformen mit symmetricher Anordnung; Symmetrieachse y_i parallel zur y-Achse')
                disp(' ');
                A_T=input('eine Teilfläche A_T:  ');
                h_z=input('z-Abstand der Schwerpunkte der Teilflächen h_z:  ');
                y_quer=input('y-Koordinate des Gesamt-Schwerpunktes:  ');
                z_quer=input('z-Koordinate des Gesamt-Schwerpunktes:  ');
                I_yT=input('axiales Flächenträgheitsmoment bezüglich der y-Achse I_yT einer Fläche:  ');
                I_zT=input('axiales Flächenträgheitsmoment bezüglich der z-Achse I_zT einer Fläche:  ');
                %Berechnung der notwendigen Größen:
                A=2*A_T;
                I_y=2*I_yT+A_T*h_z^2/2;
                I_z=2*I_zT;
                I_yz=0;
            case 4
                disp(' ');
                disp('Zwei gleiche Profile oder zwei andere gleiche Flächenformen mit symmetricher Anordnung; Symmetrieachse z_i parallel zur z-Achse')
                disp(' ');
                A_T=input('eine Teilfläche A_T:  ');
                h_y=input('z-Abstand der Schwerpunkte der Teilflächen h_y:  ');
                y_quer=input('y-Koordinate des Gesamt-Schwerpunktes:  ');
                z_quer=input('z-Koordinate des Gesamt-Schwerpunktes:  ');
                I_yT=input('axiales Flächenträgheitsmoment bezüglich der y-Achse I_yT einer Fläche:  ');
                I_zT=input('axiales Flächenträgheitsmoment bezüglich der z-Achse I_zT einer Fläche:  ');
                %Berechnung der notwendigen Größen:
                A=2*A_T;
                I_z=2*I_zT+A_T*h_y^2/2;
                I_y=2*I_yT;
                I_yz=0;
            otherwise
                disp(' ');
                disp('Ungültige Eingabe! Bitte Taste drücken und erneut versuchen.')
                pause
                fehler=1;
        end

        if fehler==0    %Abfangen ungültiger Eingaben
            Geometrie=[A, y_quer, z_quer, I_y, I_z, I_yz, Kennzahl];

            %Berücksichtigung negativer Kennzahlen (Flächen werden abgezogen)
            if Kennzahl < 1
                Geometrie=Geometrie.*[-1, 1, 1, -1, -1, -1, 1];
            end
        else
            Geometrie=[];
        end
        
    end
               
    
    function Gesamtflaeche=flaech_berech (Flaechen)
        %Flaeche=[A, y_quer, z_quer, I_y, I_z, I_yz, Kennzahl]
        
        Gesamtflaeche=Flaechen(1,:);
        Flaechen(1,:)=[];
        
        while length(Flaechen) > 0
            F1=Gesamtflaeche;
            F2=Flaechen(1,:);
            Gesamtflaeche(1)=F1(1)+F2(1);
            Gesamtflaeche(2)=(F1(1)*F1(2)+F2(1)*F2(2))/Gesamtflaeche(1);
            Gesamtflaeche(3)=(F1(1)*F1(3)+F2(1)*F2(3))/Gesamtflaeche(1);
            Gesamtflaeche(4)=F1(4)+F1(1)*(F1(3)-Gesamtflaeche(3))^2+F2(4)+F2(1)*(F2(3)-Gesamtflaeche(3))^2;
            Gesamtflaeche(5)=F1(5)+F1(1)*(F1(2)-Gesamtflaeche(2))^2+F2(5)+F2(1)*(F2(2)-Gesamtflaeche(2))^2;
            Gesamtflaeche(6)=F1(6)-F1(1)*(F1(3)-Gesamtflaeche(3))*(F1(2)-Gesamtflaeche(2))+F2(6)-F2(1)*(F2(3)-Gesamtflaeche(3))*(F2(2)-Gesamtflaeche(2));
            Gesamtflaeche(7)=0;
            Flaechen(1,:)=[];            
        end

        %Berechnung phi_stern, I1 und I2
        
        if Gesamtflaeche(4) ~= Gesamtflaeche(5)
            
            phi_stern = 0.5*atand(2*Gesamtflaeche(6)/(Gesamtflaeche(4)-Gesamtflaeche(5)));
            
        else
            
            phi_stern = 0;
            
        end
        
        I1=(Gesamtflaeche(4)+Gesamtflaeche(5))*0.5+sqrt(((Gesamtflaeche(4)-Gesamtflaeche(5))*0.5)^2+Gesamtflaeche(6)^2);
        I2=(Gesamtflaeche(4)+Gesamtflaeche(5))*0.5-sqrt(((Gesamtflaeche(4)-Gesamtflaeche(5))*0.5)^2+Gesamtflaeche(6)^2);
        
        Gesamtflaeche=[Gesamtflaeche, I1, I2, phi_stern];
                
        if phi_stern == 0
            disp(' ');
            disp('Formel für phi_stern versagt!');
        else
            disp(' ');
            disp('phi_stern wurde berechnet!');
        end


        disp(' ');
        disp(' ');
        disp('A,   y_quer,   z_quer,   I_y,   I_z,   I_yz,   I1,   I2,   phi_stern');
        disp(' ');
        Gflaeche=Gesamtflaeche;
        Gflaeche(7)=[];
        disp(Gflaeche)
        pause
    end

    function ausgeben_fl (Flaechen, Gflaeche)
        disp(' ');
        disp('Teilflächen: ');
        disp(' ');
        disp('  A,    y_quer,    z_quer,    I_y,    I_z,    I_yz,    Kennzahl');
        Flaechen
        if isempty(Gflaeche)
            disp(' ');
            disp('Gesamtflaeche nocht nicht berechnet!');
            pause
        else
            disp(' ');
            disp('Gesamtfläche: ');
            disp(' ');
            disp('  A,    y_quer,    z_quer,    I_y,    I_z,    I_yz,    I1,    I2,    phi_stern');
            Gflaeche(7)=[];
            Gflaeche
            pause
        end
        
        pause
    end

    function laden
        clear all
        load('biegno_save');
        disp(' ');
        disp('Daten wurden erfoglreich geladen!');
        pause
    end


%Lasten und Berechnung

    function Last=kraefte_eingeben
        clc
        
        disp(' ');
        disp('Bitte geben Sie Normalkraft und die Momente an!')
        disp(' ');
        N=input('Normalkraft N:  ');
        M_y=input('Moment um die y-Achse:  ');
        M_z=input('Moment um die z-Achse:  ');
        Last=[N, M_y, M_z]; 
    end

    function punkte=punkt_eingeben (punkte)
        clc
        
        disp(' ');
        disp('Bitte geben Sie einen Punkt an, in dem die Spannung berechnet werden soll!')
        disp(' ');
        Py=input('y_q-Koordinate:  ');
        Pz=input('z_q-Koordinate:  ');
        punkte=[punkte; Py, Pz];
    end

    function ausgeben (Lasten, Punkte)
        clc
        
        disp(' ');
        disp('Lasten:')
        disp(' ');
        disp(' ');
        disp('  N     M_y     M_z')
        Lasten
        disp(' ');
        disp(' ');
        disp('Punkte, in denen die Spannung berechnet werden soll: ');
        disp(' ');
        disp(' ');
        Punkte
        pause
    end

    function [sigma, m, n]=berechnung (Lasten, Punkte, zus_flaeche)
        sigma=[];
        clc
        
        if isempty(Punkte)
            disp(' ');
            disp('Es wurden keine Punkte angegeben. Bitte geben Sie jetzt einen Punkt an: ');
            disp(' ');
            Punkte=punkt_eingeben(Punkte);
        end
        
        if isempty(Lasten)
            disp(' ');
            disp('Es wurden keine Lasten angegeben. Bitte geben Sie jetzt die Lasten an: ');
            disp(' ');
            Lasten=kraefte_eingeben;
        end

         %A,    y_quer,    z_quer,    I_y,    I_z,    I_yz,    I1,    I2,    phi_stern

        %Punkte in y-z-Koordinatensystem transformieren
        punkte_t=[Punkte(:,1)-zus_flaeche(2), Punkte(:,1)-zus_flaeche(3)];
        
        %Geradengleichung der Spannungs-Null-Linie
        
        m=(Lasten(3)*zus_flaeche(4)-Lasten(2)*zus_flaeche(6))/(Lasten(2)*zus_flaeche(5)-Lasten(3)*zus_flaeche(6));
        n=-1*Lasten(1)/zus_flaeche(1)*(zus_flaeche(4)*zus_flaeche(5)-zus_flaeche(6)^2)/(Lasten(2)*zus_flaeche(5)-Lasten(3)*zus_flaeche(6));
        
        %Berechnung der Spannungen
        
        vor_y=(Lasten(3)*zus_flaeche(4)-Lasten(2)*zus_flaeche(6))/(zus_flaeche(4)*zus_flaeche(5)-zus_flaeche(6)^2);
        vor_z=(Lasten(2)*zus_flaeche(5)-Lasten(3)*zus_flaeche(6))/(zus_flaeche(4)*zus_flaeche(5)-zus_flaeche(6)^2);
        NA=Lasten(1)/zus_flaeche(1);
        
        s=size(Punkte);
        for i=1:s(1)
            temp=NA+vor_z*Punkte(i,2)-vor_y*Punkte(i,1);
            sigma=[sigma, temp];
        end
        disp(' ');
        disp('Berechnung abgeschlossen.');
        pause
        
      
    end
    
%Ausgabe und Speichern

    function ergeb_aus(sigma, Lasten, Punkte, zus_flaeche, m, n)
        clc
        disp(' ');
        disp('Geometrie: ');
        disp(' ');
        %disp('     A,     y_quer,     z_quer,     I_y,     I_z,     I_yz,     I1,     I2,     phi_stern');
        out=['Fläche A:        ', num2str(zus_flaeche(1))];
        disp(out);
        out=['y-Koord. y_quer: ', num2str(zus_flaeche(2))];
        disp(out);
        out=['z-Koord. z-quer: ', num2str(zus_flaeche(3))];
        disp(out);
        disp(' ');
        disp('axiale Flächenträgheitsmomente:');
        out=['I_y:             ', num2str(zus_flaeche(4))];
        disp(out);
        out=['I_z:             ', num2str(zus_flaeche(5))];
        disp(out);
        disp(' ');
        disp('Zentrifugalmoment:');
        out=['I_yz:            ', num2str(zus_flaeche(6))];
        disp(out);
        disp(' ');
        disp('Hauptträgheitsmomente:');
        out=['I1:              ', num2str(zus_flaeche(7))];
        disp(out);
        out=['I2:              ', num2str(zus_flaeche(8))];
        disp(out);
        disp(' ');
        disp(' ');
        disp('Lasten: ');
        disp(' ');
        disp('Normalkraft: ');
        out=['N:              ', num2str(Lasten(1))];
        disp(out);
        disp('Biegemoment um die y-Achse: ');
        out=['M_y:            ', num2str(Lasten(2))];
        disp(out);
        disp('Biegemoment um die z-Achse: ');
        out=['M_z:            ', num2str(Lasten(3))];
        disp(out);
        disp(' ');
        disp(' ');
        disp('Berechnete Spannungen: ');
        disp(' ');
        for i=1:length(sigma)
            aus=['Punkt:  ( ', num2str(Punkte(i,1)), ' , ', num2str(Punkte(i,2)), ' )   Spannung:  ', num2str(sigma(i))];
            disp(aus);
        end
        
        disp(' ');
        disp(' ');
        disp('Gleichung der Spannungs-Null-Linie: ');
        disp(' ');
        glg=['z = ', num2str(m), ' * y + (', num2str(n), ')' ];
        disp(glg);
        
        pause
    end

    function datei_aus(sigma, Lasten, Punkte, zus_flaeche)
        clc
        fid = fopen('ausgabe_biegno.txt','wt');
        fprintf(fid,'%s \n',' ');
        fprintf(fid,'%s \n','Geometrie');
        fprintf(fid,'%s \n',' ');
        out=['Fläche A:        ', num2str(zus_flaeche(1))];
        fprintf(fid,'%s \n',out);
        out=['y-Koord. y_quer: ', num2str(zus_flaeche(2))];
        fprintf(fid,'%s \n',out);
        out=['z-Koord. z-quer: ', num2str(zus_flaeche(3))];
        fprintf(fid,'%s \n',out);
        fprintf(fid,'%s \n',' ');
        fprintf(fid,'%s \n','axiale Flächenträgheitsmomente: ');
        out=['I_y:             ', num2str(zus_flaeche(4))];
        fprintf(fid,'%s \n',out);
        out=['I_z:             ', num2str(zus_flaeche(5))];
        fprintf(fid,'%s \n',out);
        fprintf(fid,'%s \n',' ');
        fprintf(fid,'%s \n','Zentrifugalmoment: ');
        out=['I_yz:            ', num2str(zus_flaeche(6))];
        fprintf(fid,'%s \n',out);
        fprintf(fid,'%s \n',' ');
        fprintf(fid,'%s \n','Hauptträgheitsmomente: ');
        out=['I1:              ', num2str(zus_flaeche(7))];
        fprintf(fid,'%s \n',out);
        out=['I2:              ', num2str(zus_flaeche(8))];
        fprintf(fid,'%s \n',out);
        fprintf(fid,'%s \n',' ');
        fprintf(fid,'%s \n',' ');
        fprintf(fid,'%s \n','Lasten:  ');
        fprintf(fid,'%s \n',' ');
        fprintf(fid,'%s \n','Normalkraft:  ');
        out=['N:              ', num2str(Lasten(1))];
        fprintf(fid,'%s \n',' ');
        fprintf(fid,'%s \n','Biegemoment um die y-Achse:  ');
        out=['M_y:            ', num2str(Lasten(2))];
        fprintf(fid,'%s \n',out);
        fprintf(fid,'%s \n','Biegemoment um die z-Achse:  ');
        out=['M_z:            ', num2str(Lasten(3))];
        fprintf(fid,'%s \n',out);
        fprintf(fid,'%s \n',' ');
        fprintf(fid,'%s \n',' ');
        fprintf(fid,'%s \n','Berechnete Spannungen:  ');
        fprintf(fid,'%s \n',' ');
        for i=1:length(sigma)
            aus=['Punkt:  ( ', num2str(Punkte(i,1)), ' , ', num2str(Punkte(i,2)), ' )   Spannung:  ', num2str(sigma(i))];
            fprintf(fid,'%s \n',out);
        end
        
        fprintf(fid,'%s \n',' ');
        fprintf(fid,'%s \n',' ');
        fprintf(fid,'%s \n','Gleichung der Spannungs-Null-Linie:  ');
        fprintf(fid,'%s \n',' ');
        glg=['z = ', num2str(m), ' * y + (', num2str(n), ')' ];
        fprintf(fid,'%s \n',glg);
        
        pause
    end
    
    function speichern(sigma, Lasten, Punkte, zus_flaeche)
        clc
        disp(' ');
        disp('Geometrie- und Berechnungsdaten werden werden in der Datei biegno_save.mat gespeichert.');
        disp(' ');
        disp('ACHTUNG: Falls die Datei biegno_save.mat bereits existiert, wird diese überschrieben. ');
        pause
        save('biegno_save', 'sigma', 'Lasten', 'Punkte', 'zus_flaeche');
        clc
        disp(' ');
        disp('Daten gespeichert!');
        pause
    end

        
clc
end %Biegno