function [ M_brow, M_dir, M_conf, M_imm, Stat, Ext, Dcoeff, S_mss, pourc, Resume, EXT, beta] = motion_classification_simu(TRACKS,  T, t_msd, tau, tau_space, name, model)

% It computes the classification of every sub-tracks
% 
% INPUT
% - TRACKS : structure array containing the coordinates of tracks (TRACKS(i).coord)
% - T      : lenght of the sub-tracks
% - t_msd  : parameter to clalculate the linear regression of the
%            msd-fucntion on [1, t_msd]
%
% OUTPUT
% - M_brow (dir, conf, imm) : vector containing the number (wrt the order
%                             of TRACKS) of brownian (directed, confined, immobile) tracks
% - Stat                    : value of the statistical test for each track of TRACKS
% - Ext                     : value of the extent for each track of TRACKS
% - Stat_mss                : slope of the linear regression of  MSS-function
% - D_msd                   : slope of the linear regression of msd-function (for brownian motion this corresponds to the diffusion coeff
% - pourc                   : percentage for every class of mouvement (see line 105)
% - Resume                  : # track | # track orig | mvmnt | Ext | D_msd | S_mss


n_tracks= size(TRACKS,2);
T=size(TRACKS(1).coord,2);

M_brow=[];
M_dir=[];
M_conf=[];
M_imm=[];
M_NC=[];
Stat=[];
Dcoeff=[];
beta=[];
msd=[];
Dcoeff=[];
SqDispl=[];
pourc=[];
test_imm=[];
EXT=[];
Ext=[];
lenght= num2str(T-1);
Resume= zeros(n_tracks,6);

% quantiles of order alpha/2 and (1-alpha/2) with alpha =0.05
% quantiles depend on the lenght of the trajectories !

switch  lenght  
    case '3'
     q1 = .743;
     q2 = 2.247;
    case '5'
     q1=  .724; %.731;
     q2=  2.464; % 2.379;
    case '10'
     q1=.725;
     q2=2.626; 
    case '30'
     q1=.754;
     q2=2.794;
end
q1
q2


n_analysed_tracks=size(TRACKS,2);

for i=1:n_tracks
    
    X=tau_space.*(TRACKS(i).coord);
          
    Y = X(:,1:end) - repmat(X(:,1),1,T);
    S = max(sqrt(sum(Y.*Y,1)));
                   
    Z = X(:,2:end)-X(:,1:end-1);
    sigma = sqrt(sum(sum(Z.*Z))/2);
    Stat(i) = S/sigma; 

%% Extent

    
      test_immo=[];
      count=1;
      for l=1:size(X,2)
          for k=l:size(X,2)
              test_immo(count)=norm(X(:,l)-X(:,k),2);
              count=count+1;
          end
      end
      
      Ext(i) = max(test_immo);
      extent = max(test_immo);
      
      EXT=[EXT, [extent;i;i]];




%% MSS
        kmax=floor(T/3)+1;
        MSS=zeros(kmax,1);
        
        for k=1:kmax
        moment(k).coord=zeros(T-1,1);
        for j=1:T-1   
        A=X(:,j+1:T)-X(:,1:T-j);
        moment(k).coord(j) = (sum(sqrt(diag(A'*A)).^(k))/(T-j));        
        end
        moment(k).coord=((1/tau)^(k)).*moment(k).coord; % change of variables t->tau*t
        ymss= log(moment(k).coord);
        tmss=[ones(T-1,1) ,log(tau.*[1:T-1]')];
        nu=tmss\ymss; 
        MSS(k) = real(nu(2));
        end

        nu_bis=[ones(kmax,1),[1:kmax]']\MSS; 
        S_mss(i)=real(nu_bis(2));

        
%% Classification

switch model
    
    case 'Stat'
        %%Diffusion coefficient
        
        SqDispl = sum((X(:,2:end)-X(:,1:end-1)).^2,1);
        Dcoeff(i) = mean(SqDispl)/(4*tau);
        taille=2;
        % Classif
        if Stat(i)< q1
            type=1;
            M_conf = [M_conf, i];
        elseif Stat(i) > q2
            type=2;
            M_dir = [M_dir, i];
        elseif Stat(i) >= q1 & Stat(i) <= q2
            type=3;
            M_brow=[M_brow, i];
        end

        
    case 'MSD'
        
        % Linear regression of msd-function on [1,t_msd]
        msd(i).coord = zeros(1, T-1);
        
        for j=1:T-1
            msd(i).coord(j) = sum(sum((X(:,j+1:T)-X(:,1:T-j)).^2))/(T-j);
        end
        t=[ones(t_msd,1) ,log(tau.*[1:t_msd]')];
        y=log(msd(i).coord(1:t_msd)');
        param=t\y;
        beta(i)=param(2);
        
        t=[ones(t_msd,1) ,tau.*[1:t_msd]'];
        y=msd(i).coord(1:t_msd)';
        param=t\y;
        Dcoeff(i)=param(2)/4;
        
        % Classif
        if  0<=Dcoeff(i) & Dcoeff(i)< .1
            type=4;
            M_imm = [M_imm, i];
        elseif .1<=Dcoeff(i) & Dcoeff(i)< .9
            type=1;
            M_conf = [M_conf, i];
        elseif Dcoeff(i)>=1.1
            type=2;
            M_dir = [M_dir, i];
        elseif .9<=Dcoeff(i) & Dcoeff(i)< 1.1
            type=3;
            M_brow=[M_brow, i];
        else % nonclassified
            type=5;
        end
end



%% Resume (# track | # track orig | mvmnt | Ext | D_msd | S_mss)

   Resume(i,1)=i;
   Resume(i,2)=i;
   Resume(i,3)=type;
   Resume(i,4)=Ext(i);
   Resume(i,5)=Dcoeff(i);
   Resume(i,6)= S_mss(i);
end
    pourc(1)=numel(M_conf)*100/n_analysed_tracks;
    pourc(2)=numel(M_dir)*100/n_analysed_tracks;
    pourc(3)=numel(M_brow)*100/n_analysed_tracks;
    pourc(4)=numel(M_imm)*100/n_analysed_tracks;
    pourc(5) = n_analysed_tracks;
    
%     disp('(% conf, % dir, % brow, % imm)')
%     disp(pourc)



end

