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

  %  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;
%   tspan = [0 20];
%     y0 = [0,0,0,0 ; 0,0,0,0 ; 0,0,0,0 ; 0,0,0,0];
%
%     options=odeset('mass','M)')
% %[t,y]=ode113('indmot_ode1',tspan,y0,options)
% [T, Y] = ode45('diffgl', tspan, y0, options);
%
%
%  function dy = diffgl(t, y, flag)
%         switch flag
%             case ''
%                 dy = zeros(4,1);%-RK+FD+FE
%                 dy(1) = -k2v*z(1) +k2v*z(3)   -d2v(zd2v-zd1v)+0;
%                 dy(2) =  - k2h*z(2) +k2h*z(4) -d2h(zd2h-zd1h)+0;
%                 dy(3) =  +k2v*z(1) -k2v*z(3)  -d2v(zd1v-zd2v)-d1v(zd1v-hdv )+k1v(z1v-hv);
%                 dy(4) =  k2h*z(2)  - k2h*z(4) -d2h(zd1h-zd1h)-d1h(zd1v-hdh)+k1h(z1h-hh);
%             case 'mass'
%           dy = RM;
%         end
%     end
% 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 = [0 20];
    y0 = [0,0,0,0 ; 0,0,0,0 ; 0,0,0,0 ; 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)

                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);
                dy(4) =   k2h -k2h  -d2h*(zd1h-zd1h)-d1h*(zd1v-hdh)+k1h*(z1h-hh);

    end



% function varargout=indmot_ode1(t,y,flag)
% switch flag
% case ''
% %no input flag
% varargout{1}=FF1(t,y);
% case 'mass'
% %flag of mass calls mass.m
% varargout{1}=MM1(t,y);
% otherwise
% error(['unknown flag ''' flag '''.']);
% end
%
%
%     function yp = FF1(t, y)
%                 dy = zeros(4,1);%-RK+FD+FE
%                 dy(1) = -k2v*z(1) +k2v*z(3)   -d2v(zd2v-zd1v)+0;
%                 dy(2) =  - k2h*z(2) +k2h*z(4) -d2h(zd2h-zd1h)+0;
%                 dy(3) =  +k2v*z(1) -k2v*z(3)  -d2v(zd1v-zd2v)-d1v(zd1v-hdv )+k1v(z1v-hv);
%                 dy(4) =  k2h*z(2)  - k2h*z(4) -d2h(zd1h-zd1h)-d1h(zd1v-hdh)+k1h(z1h-hh);
%     end
%
%  function n = MM1(t, y)
%                 dy = RM%zeros(4,1); %RM
%                 %n1 = (m2v +m2k(lh/lg)^2)
%
%             end

end