function ex1(profil, schrittweite, geschwindigkeit, m, l, k, d)
% Implementieren Sie das Modell mit MATLAB, die dazu notwendigen Gleichungen sind
% im Anhang beschrieben. Führen Sie einen Simulationslauf mit den angegebenen
% Parametern durch, plotten Sie die Verläufe sämtlicher Zustandsgrößen über der
% Zeit und interpretieren Sie die Ergebnisse. Die Fahrbahnunebenheiten sind dabei
% in tabellarischer Form mit Höhenangaben an charakteristischen Punkten entlang der
% Fahrbahn gegeben. Dazwischen soll der Verlauf linear interpoliert werden. Sie
% werden auch die Ableitung dieser so definierten Funktion benötigen.

% Wie kann diese in einfacher Weise aus den Tabellenwerten errechnet werden?
%  Ableitung = Steigung der Tangente -> (profil[i+1] - profil[i]) / (Distanz zw. i und i+1)

% Differenzialgleichungen: RM * z'' = FD + FE - RK * z
%    RM .. Massematrix
%    RK .. Steifigkeitsmatrix
%    FD .. Vektor mit den Dämpfkräften der Federn
%    FE .. Vektor mit den Erregerkräften

% Beispiel aus ACSL Programm:
%   >> m1v = 60.0; m1h = 80.0; m2v = 1200.0; m2h = 1200.0; m2k = 1200.0; lv = 1.0; lh = 1.5; lg = lv + lh;
%   >> k1h = 260.e3; k1v = 260.e3; k2h = 50.e3; k2v = 40.e3; d1h = 0.2e3; d1v = 0.2e3; d2h = 0.0; d2v = 0.0;
%   >> ex1(m1v, m1h, m2v, m2h, m2k, lv, lh, lg, k1h, k1v, k2h, k2v, d1h, d1v, d2h, d2v)

    m1v = m(1);
    m1h = m(2);
    m2v = m(3);
    m2h = m(4);
    m2k = m(5);

    lv = l(1);
    lh = l(2);
    lg = lv + lh;

    k1v = k(1, 1);
    k1h = k(1, 2);
    k2v = k(2, 1);
    k2h = k(2, 2);

    d1v = d(1, 1);
    d1h = d(1, 2);
    d2v = d(2, 1);
    d2h = d(2, 2);

    deltaT = schrittweite / geschwindigkeit;

    % Erstelle Vektor mit dem Fahrbanprofil in Abhängigkeit von der Schrittweite -> h
    x = [profil(1 , 1) profil(1 , 2)];
    y = [profil(2 , 1) profil(2 , 2)];
    xi = profil(1 , 1) : schrittweite : profil(1 , 2);
    h = interp1(x, y, xi);
    for i = 2 : (length(profil) - 1)
        x = [profil(1 , i) profil(1 , i + 1)];
        y = [profil(2 , i) profil(2 , i + 1)];
        xi = profil(1 , i) : schrittweite : profil(1 , i + 1);
        hi = interp1(x, y, xi);
        h = [h , hi(2 : length(hi))];
    end;

    clear x y xi;

    x = 0 : schrittweite : profil(1, length(profil)) + lg;
    hh = [zeros(size(0 : schrittweite : lg)) h];
    hv = [h zeros(size(0 : schrittweite : lg))];
    dhh = diff(hh);
    dhv = diff(hv);
    hv(length(hv)) = []; % Make all arrays equal length
    hh(length(hh)) = []; % Make all arrays equal length
    x(length(x)) = [];
    x(length(x)) = [];

    % Plotte Steigung der Fahrbahn (1. Ableitung)
%    plot(x, dhh, '-r', x, hh, '-b');

    % Definiere die Matritzen für die systembeschreibenden Diffenrentialgleichungen:  RM * z'' + RK * z = FD + FE
    RM = [m2v+m2k*(lh/lg)*(lh/lg) , m2k*(lv*lh/(lg*lg)) , 0 , 0 ; m2k*(lv*lh/(lg*lg)) , m2h+m2k*(lv/lg)*(lv/lg), 0 , 0 ; 0 , 0 , m1v , 0 ; 0 , 0 , 0 , m1h]
    RK = [k2v , 0 , -k2v , 0 ; 0 , k2h , 0 , -k2h ; -k2v , 0 , k2v , 0 ; 0 , -k2h , 0 , k2h];
    % Modell mit linearen Dämpfern
%    FD = [-d2v*(dz2v-dz1v) , -d2h*(dz2h-dz1h) , -d2v*(dz1v-dz2v)-d1v*(dz1v-dhv) , -d2h*(dz1h-dz2h)-d1h*(dz1h-dhh)]';
    % Modell mit konstanten Dämpfern
    FD = [-d2v , -d2h , -d2v-d1v , -d2h-d1h]';
    FE = [0 , 0 , k1v*hv , k1h*hh]';
%    FE = [0 , 0 , k1v*(z1v-hv) , k1h*(z1h-hh)]';

% Umformen: RM * z'' = -RK*z + FD + FE
% Obige Matrixschreibweise in das DGL 1. Ordnung auflösen und in die Funktion unten einfügen.
% u = dz/dt    RM * u' = -RK*z + FD + FE
% Danach mit einem ODE Solver lösen (ode23, ode45, ...)

% Löse das DGL mit Schrittweite durch x definiert
zd2v =0;
zd1v =0;
zd2h = 0;
zd1h = 0;
hdv =0;
hdh=0;
z1h=0;
z1v=0;

    tspan = [1 20];
    y0 = [0 ; 0 ; 0 ; 0];
    opt = odeset('Mass', RM);
    [T, Y] = ode45(@diffgl, tspan, y0, opt);

    hold on;
    plot(T,Y(:,1));
%    plot(T, Y(:,1), '-r', T, Y(:,2), '-g' T, Y(:,3), '-b' T, Y(:,4), '-');
    grid off;
    xlabel('t');
    ylabel('Auslenkung');
    legend('Aufhängung vorne', 'Aufhängung hinten', 'Radmittelpunkt vorne', 'Radmittelpunkt hinten');
    hold off;


    function dy = diffgl(t, y,flag)

                dy = zeros(4,1);%-RK+FD+FE
                dy(1) =  -k2v +k2v  -d2v*(zd2v-zd1v);
                dy(2) =  -k2h +k2h  -d2h*(zd2h-zd1h);
                dy(3) =  +k2v -k2v  -d2v*(zd1v-zd2v)-d1v*(zd1v-hdv )+k1v*(z1v-hv(t));
                dy(4) =   k2h -k2h  -d2h*(zd1h-zd1h)-d1h*(zd1v-hdh)+k1h*(z1h-hh(t));

    end

    end