function [t,y,breakIndices] = runSSystem2( Lss,simData)
%RUNSSYSTEM2 integrates a linear switched system with output-dependent switching
% Given an output-dependent switching function, integrate a linear switched
% system with a given control function with the implicit euler method
%
% Input arguments:
%   Lss:        Struct containing the switched system
%       Lss.dim         (int) state dimension (n) of the system matrices 
%                       (must be the same for all subsystems)
%       Lss.dimInput    (int) dimension (m) of the input function 
%                       (must be the same for all subsystems)
%       Lss.dimOutput   (int) dimension (p) of the output function 
%                       (must be the same for all subsystems)
%       Lss.nSystems    (int) number of subsystems
%       Lss.As          cell array of length FOM.nSystems containing the A
%                       matrices, i.e. FOM.As = {A_1,...,A_{nSystems}}
%       Lss.Bs          cell array of length FOM.nSystems containing the B
%                       matrices, i.e. FOM.Bs = {B_1,...,B_{nSystems}}
%       Lss.Cs          cell array of length FOM.nSystems containing the C
%                       matrices, i.e. FOM.Cs = {C_1,...,C_{nSystems}}
%   simData:    Struct containing all information for the time simulation
%       simData.switchingConditions:     1d array containing the switching order
%                                   Ex: [1 2 1] start with mode 1, switch
%                                   to mode 2 and then back to mode 1
%       simData.u:                  the control function for the LSS
%       simData.initialCondition:   the initial condition
%       simData.startSystem:        the number of the mode to start with
%       simData.timeStep:           the step size for the time
%                                   discretization
%       simData.maxTime:            integrate the system from [0 maxTime]
%
% Output arguments:
%   t:              the time points used in the integration
%   y:              the approximation of the ouput at the time points
%   breakIndices:   the indices at which the switching happens
%   xx:             the approximation of the state variable at the time
%                   points
%
% Authors:
%   Philipp Schulze, pschulze@math.tu-berlin.de
%   Benjamin Unger, unger@math.tu-berlin.de
%
% Version:
%   1.0 (January 08, 2018)

    %% create shortcuts for the required Data
    As                  = Lss.As ;
    Bs                  = Lss.Bs ;
    Cs                  = Lss.Cs ;
    Ks                  = Lss.Ks ;
    n                   = Lss.dim ;
    p                   = Lss.dimOutput ;

    switchingConditions = simData.switchingConditions ;
    u                   = simData.u ;
    x0                  = simData.initialCondition ;
    startSystem         = simData.startSystem ;
    dt                  = simData.timeStep ;
    T                   = simData.maxTime ;
    
    %% initialization of Variables
    breakIndices        = [] ;
    t                   = 0:dt:T ;
    nt                  = length(t) ;
    x                   = zeros(n, nt) ;
    x(:,1)              = x0 ;
    mode                = startSystem ;
    y                   = zeros(p, nt) ;
    y(:,1)              = Cs{mode}*x0 ;
    
    %% time integration with backward Euler
    for i=2:nt
        [switchAction, switchTo]    = switchingConditions(y(:,i-1), mode) ;
        if(switchAction)
            x(:,i-1)                = Ks{mode}{switchTo}*x(:,i-1) ;
            mode                    = switchTo ;
            breakIndices            = [breakIndices i] ;
        end
        x(:,i)      = (eye(n)-As{mode}*dt)\(x(:,i-1)+dt*Bs{mode}*u(t(i)));
        y(:,i)      = Cs{mode}*x(:,i);
    end
    breakIndices    = [breakIndices, nt] ;
