function  [obsres,simres] = RAPTOR_EKF_sim(G,V,U,H,model,params,KFmodel,KFparams)

% initial plasma state derivative (for solver)

if isa(KFparams,'Simulink.Parameter')
    KFp = KFparams.Value;
else
    KFp = KFparams;
end

zk_km1 = [KFp.xk_0;KFp.dk_0];

% initialize xdot
xdotk = initialize_xdot(model,params);

nt = numel(params.tgrid);

ink = zeros(KFmodel.ny,nt);
res = zeros(1,nt);
exitflag = zeros(1,nt);

Xpred  = zeros(KFmodel.nx,nt);
Hhat   = zeros(KFmodel.ny,nt);
Xhat   = zeros(KFmodel.nx,nt); % state estimate
Xk_km1 = zeros(KFmodel.nx,nt);
Xk1_k  = zeros(KFmodel.nx,nt);
Xdotk1 = zeros(KFmodel.nx,nt);
Dhat   = zeros(KFmodel.nd,nt);
Shat   = zeros(KFmodel.nz,KFmodel.nz,nt);
Fkhat  = zeros(KFmodel.nz,KFmodel.nz,nt);
Gkhat  = zeros(KFmodel.nz,KFmodel.nw,nt);
Hkhat  = zeros(KFmodel.ny,KFmodel.nz,nt);
srcs   = cell(nt,1);

traps = cell(1,nt);
staps = cell(1,nt);
geops = cell(1,nt);

for it = 1:nt;
    % geometry and kinetic profiles
    gk = G(:,it);
    vk = V(:,it);
    uk = U(:,it);
    hk = H(:,it);
    
    % state covariances
    if it==1
        [Qx,Qd,Sx0,Sd0] = KFcovariances(KFmodel,KFp,model);
        Qk = blkdiag(Qx,Qd);
        Sk_km1 = blkdiag(Sx0,Sd0);
    end

    % state at previous time step
    xk_km1 = zk_km1(KFmodel.ix);
  
    %%
    diagfun = model.diagnostics.diagfun;
    [hk_km1,dhk_dx,~,Rk,~,~] = diagfun(xk_km1,gk,vk,uk,it,model);
    Hk = [dhk_dx,zeros(numel(hk_km1),KFmodel.nd)];
    
    %%
    % cla; plot(hk_km1,'ob'); hold on; plot(hk,'xr'); drawnow
    % determine dt
    if it~=nt
        dt = params.tgrid(it+1)-params.tgrid(it);
    else
        dt = params.tgrid(it)-params.tgrid(it-1); % extrap last one
    end
    
    % time derivative of kinetic profiles, enters in some equation
    if it==1
        vdot = zeros(size(vk));
    else
        vdot = (V(:,it)-V(:,it-1))/dt;
    end
    
    % time derivative of geometry profiles, enters in some equations
    if it==1
        gdot = zeros(size(gk));
    else
        gdot = (G(:,it)-G(:,it-1))/dt;
    end
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Observer step - kalman filter
    % This part is common for RT and off-line part
     [zk1_k,zk_k,xdotk1,Sk1_k,Sk_k,~,~,Fk,Gk,stap,geop,trap,res(it),icrash,exitflag(it)] = ...
            RAPTOR_EKF_step(zk_km1,xdotk,gk,gdot,vk,vdot,uk,hk,hk_km1,Hk,Sk_km1,Qk,Rk,dt,it,model,params,KFmodel,KFp);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    xk_k = zk_k(1:KFmodel.nx);
    dk_k = zk_k((KFmodel.nx+1):end);
    
    %update the equilibrium from chease
    if model.equi.updateequi > 0
        gk_new = update_g_from_CHEASE(xk_k,gk,vk,xdotk1,uk,stap,geop,trap,it,dt,model,params);
        G(:,it:end) = repmat(gk_new,1,size(G,2)-it+1); %use this
        % update the equilibrium structure in diagnostics that
        % require it.
        model.diagnostics.diagfun = update_diagnostic_equilbria(xk_k,G,V,U,it,model, params.tgrid);
    end
    
    % storage
    ink(:,it) = sqrt(Rk)\(hk-hk_km1); % innovation sequence

    Hhat(:,it)   = hk_km1;
    Xhat(:,it)   = xk_k; % state estimate
    Xk_km1(:,it) = xk_km1;
    Xk1_k(:,it)  = zk1_k(KFmodel.ix);
    Xdotk1(:,it) = xdotk1;
    
    staps{it} = stap;
    geops{it} = geop;
    traps{it} = trap;

    
    Dhat(:,it)   = dk_k; % distubrance estimate (physical units)
    Shat(:,:,it) = Sk_k; % state error covariance
    Fkhat(:,:,it)   = Fk;   % state transition matrix;
    Gkhat(:,:,it)   = Gk;   % input matrix;
    Hkhat(:,:,it)   = Hk;   % output matrix;
    if it~=nt
        Xpred(:,it+1) = zk1_k(1:KFmodel.nx); % one-step-ahead predicted state
    end
    
    % next step
    Sk_km1 = Sk1_k;
    zk_km1 = zk1_k;
    xdotk = xdotk1;    
end

simres.X = Xhat;
simres.G = G;
simres.V = V;
simres.Xdot = Xdotk1;
simres.U = U;
simres.staps = staps;
simres.traps = traps;
simres.geops = geops;

obsres.Xhat = Xhat;
obsres.Xk_km1 = Xk_km1;
obsres.Xk1_k = Xk1_k;
obsres.Xdotk1 = Xdotk1;
obsres.Dhat = Dhat;
obsres.Shat = Shat;
obsres.Hhat = Hhat;
obsres.Fkhat = Fkhat;
obsres.Gkhat = Gkhat;
obsres.Hkhat = Hkhat;
obsres.res = res;
obsres.exitflag = exitflag;
obsres.ink = ink;
obsres.Qk = Qk;
obsres.Rk = Rk;
% obsres.Xpred = Xpred;
% obsres.Ypred = Ypred;
% obsres.Innov = Innov;
% obsres.Omega = Omega;

return
% 
% function KF_plot_matrices(Fk,Hk,Sk_km1,Sk_k,Sk1_k,Qk,Rk,Omk,Lk,KFmodel)
% 
% LH = Lk*Hk';
% 
% nz = KFmodel.nz;
% 
% figure(1); clf; % set(gcf,'position',[100 600 400 1200]);
% hax = multiaxes(3,4,[],[],[]);
% % plot covariance matrices
% 
% ax = hax(1);
% imagesc(Hk'*Sk_km1*Hk,'parent',ax); title(ax,'H_k^TS_{k|k-1}H_k'); colorbar('peer',ax);
% 
% ax = hax(2);
% imagesc(Rk,'parent',ax); title(ax,'R_k'); colorbar('peer',ax);
% 
% ax = hax(3);
% imagesc(Omk,'parent',ax); title(ax,'\Omega_k'); colorbar('peer',ax);
% 
% ax = hax(4);
% imagesc(Sk_km1,'parent',ax); title(ax,'S_{k|k-1}'); colorbar('peer',ax);
% 
% ax = hax(5);
% imagesc(eye(nz) - LH(1:nz,1:nz),'parent',ax); title(ax,'-L_kH_k^TS_{k|k-1}'); colorbar('peer',ax);
% 
% ax = hax(6);
% imagesc(Sk_k,'parent',ax); title(ax,'S_{k|k}'); colorbar('peer',ax);
% 
% ax = hax(7);
% imagesc(Fk(1:nz,1:nz)*Sk_k(1:nz,1:nz)*Fk(1:nz,1:nz)','parent',ax); title(ax,'F_kS_{k|k}F_k^T'); colorbar('peer',ax);
% 
% ax = hax(8);
% imagesc(KFmodel.Gk*Qk*KFmodel.Gk','parent',ax); title(ax,'G_kQ_kG_k^T'); colorbar('peer',ax);
% 
% ax = hax(9);
% imagesc(Sk1_k,'parent',ax); title(ax,'S_{k+1|k}'); colorbar('peer',ax);
% 
% ax = hax(9);
% imagesc(Sk1_k,'parent',ax); title(ax,'S_{k+1|k}'); colorbar('peer',ax);
% 
% ax = hax(10);
% imagesc(Lk,'parent',ax); title(ax,'L_k'); colorbar('peer',ax);
% 
% set(hax(1:9),'xticklabel','','yticklabel','');
% 
% return
