function s_arr=f_StateSequence(Q,dt,p0,Nsteps,Nr)

% Given Markovian model "Q" (single channel or composite), computes for
% "Nsteps" successive steps of duration "dt" and for "Nr" sweeps
% the state of the model at each time step and for each sweep,
% given initial probability distribution p0

% Input
% Matrix Q of the model
% dt
% Initial probability p0 as column vector (if empty [], will be set to steady state)
% Number Nsteps of time steps (output will include the first and the last value)
% Number of independent realizations Nr (sweeps)

% Output:
% s_arr: Nr-row by Nsteps+1 columns matrix of state indices

% Compute TPM (transition probability matrix)
A=expm(Q*dt); % TPM, left stochastic matrix

% Equilibrium (steady state distribution), if p0 not given explicitly
if isempty(p0)
    p0=findSS(Q);
end

nstates=length(p0);

% Calculate cumulated sums of the columns of the TPM A
CA=cumsum(A);
CA=CA*diag(1./CA(nstates,:)); % normalize to be sure that the sum of each column equals 1

% Allocate array for the result
s_arr=zeros(Nr,Nsteps+1); % includes initial and final state

% Loop through the Nr realizations/sweeps
for i=1:Nr % for each realization / sweep
    
    % Compute initial state
    % a number between 1 and Nstates, drawn from the multinomial distribution given by p0
    s=find(rand()<cumsum(p0),1); 
    
    % Record the initial state
    s_arr(i,1)=s;
    
    % Loop for each time step
    for j=2:Nsteps+1
        s=find(rand()<CA(:,s),1); % transition to the state at the next time step
        s_arr(i,j)=s; % Record s in the array
    end
    
end
