sources_mvt(X_ini , F , G , Q , dt , L);

function X = sources_mvt(X_ini , F  , G , Q , dt , L);

%  Generate Target trajectories
%
%  X = sources_mvt(X_ini , F  , G , Q , dt);
%
% Inputs
% -----
%
% X_ini          (nx x 1 x M) initial state vector
% F              inline transition matrix (nx x nxM). F depends from delta_t.
% G              inline covariance noise process (nx x nxM).
% Q              constant covariance noise process (nx x nxM).
% dt             (1 x T) or (T x 1) time vector (t_k - t_{k - 1}).
%
% Optionnal input 
%
%  X = sources_mvt(X_ini , F  , G , Q , dt , L);
%
% L              number of Mth sources to generate (default L = 1).
%
% Output
% -----
%
% X              (nx x 1 x M x T x L)
%
%***************************************************************************%%
%  Auteur Sbastien PARIS ([e-mail]), Septembre 2002                   %
%***************************************************************************%%


if (nargin < 6)
   
   L                  = 1;
   
end


T                         = length(dt);

[nx , q , M]              = size(X_ini);

%------------------------ Allocation -------------------------%

X                         = zeros(nx , 1 , M , T , L);                

W                         = randn(nx , 1 , M , T - 1 , L);

%-------------------- Initialisation -----------------%

X(: , : , : , 1 , :)      = X_ini(: , : , : , ones(1 , L));  

%-------------------- Itrations ---------------------%


for t  = 2 : T
        
  X(: , : , : , t , :)    = ndtimes(F(dt(t - 1)), X(: , : , : , t - 1 , :)) + ndtimes(permute(ndchol(G( dt(t - 1) ).*Q ) , [2 1 3]) , W(: , : , : , t - 1 , :));

end