function est = kest_asyn(N,M,A,C,Q,R,x0_,x0,P0,T,u,sys,t)

% ASYMPTOTIC
fprintf('\nAsymptotic Kalman filtering\n')
vai1 = 0; vai2 = 0;
r=rank(obsv(A,C));
if r == length(A) %Coppia osservabile
    fprintf('(A,C) oss => (A,C) riv => P(t) converge ad una soluzione s.d.p. di ARE\n')
    vai1 = 1;
else
    [Abar,Bbar,Cbar,Tt,k] = obsvf(A,0,C);
    A_non_obs=Abar(1:(length(A)-r),1:(length(A)-r));    %Estraggo la matrice relativa allo spazio non osservabile
    d = eig(A_non_obs);         %Determino i suoi autovalori
    if length(find(abs(d)>=1)>0)%controllo che stiano tutti strettamente dentro al cerchio unitario
        vai1 = 0;               
        fprintf('(A,C) NON riv\n');
    else
        fprintf('(A,C) riv =>  P(t) converge ad una soluzione s.d.p. di ARE\n');
        vai1 = 1;
    end
end

r=rank(ctrb(A,sqrtm(Q)));
if  r == length(A) %Coppia raggiungibile
    fprintf('(A,sqrtm(Q)) ragg => (A,sqrtm(Q)) stab => esiste unica P_inf soluzione di ARE d.p. stabilizzante\n')
    vai2 = 1;
else
    [Abar,Bbar,Cbar,Tt,k] = ctrbf(A,sqrtm(Q),0);
    A_non_ragg=Abar(1:(length(A)-r),1:(length(A)-r));%Estraggo la matrice relativa allo spazio non raggiungibile
    d = eig(A_non_ragg);         %Determino i suoi autovalori
    if length(find(abs(d)>=1)>0)%controllo che stiano tutti strettamente dentro al cerchio unitario
        vai2 = 0;                
        fprintf('(A,sqrtm(Q)) NON stab\n');
    else
        fprintf('(A,sqrtm(Q)) stab => esiste unica P_inf soluzione di ARE s.d.p. stabilizzante\n');
        vai2 = 1;
    end
end

if vai1 & vai2
    e_post=zeros(M,4,N);
    fprintf('Compute parameters...')
    [P,eigenvalues,L] = dare(A',C',Q,R); L = L';
    K=A*P*C'*pinv(C*P*C'+R);
    fprintf('done\n')
    fprintf('Compute estimate...')
    for m=1:M
        [y,t,x] = lsim(sys,squeeze(u(m,:,:)),t,x0); y = y'; x = x';
        x0_post = x0_ + P0*C'*pinv(C*P0*C' + R)*( y(:,1) - C*x0_ ); % da xo a priori a x0 a posteriori
        x_post=x;
        x_post(:,1)=x0_post;
        for j=2:N
           x_post(:,j)=A*x_post(:,j-1)+K*(y(:,j)-C*A*x_post(:,j-1));
        end
        e_post(m,:,:)=x-x_post(:,1:end); %Errore di stima a posteriori
    end
    fprintf('done\n')
    est.P_prior = 0;
    est.P_post = P;
    est.L = L;
    est.x_prior = 0;
    est.x_post = x_post;
    
    e=squeeze(e_post(M,:,:));
    % Calcolo della varianza campionaria dell'errore di filtraggio tramite algoritmo NON RICORSIVO 
    est.P_emp_post=zeros(4,4,N);
    for j = 1:N
         X=e(:,1:j)'-ones(j,1)*mean(e(:,1:j)',1);
         est.P_emp_post(:,:,j)=(X'*X)./j;
    end
    % Calcolo della varianza campionaria dell'errore di filtraggio SU M SIMULAZIONI tramite algoritmo NON RICORSIVO 
    est.PM_post=zeros(4,4,N);
    for j = 1:N
         X=squeeze(e_post(:,:,j))-ones(M,1)*mean(squeeze(e_post(:,:,j)),1);
         est.PM_post(:,:,j)=(X'*X)./M;
    end
    est.NPM_post=zeros(N,1);
    est.NP_emp_post=zeros(N,1);
    for j = 1:N
        est.NP_emp_post(j)=trace(est.P_emp_post(:,:,j));
        est.NPM_post(j)=trace(est.PM_post(:,:,j));
    end
    est.OK=1;
else
    fprintf('Can''t proceed...\n');
    est.OK=0;
end



