function [L,LX,LY] = liuqekalman(shot,t,varargin)
%
%LIUQEKALMAN  Development for estimating vessel currents with Kalman filter
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

 % Open shot
 assert(nargin >= 1,'LIUQE:Call:Specify shot')
 assert(~isempty(mdsopen(shot)),'LIUQE:MDS:Cannot open shot #%d',shot)

 % LIUQE parameters
 L.P = liuptcv(varargin{:});
 % LIUQE geometrical data
 L.G = liugtcv;
 Msa = mdsdata('STATIC("MUT_S_A")[$1,$2]',L.G.dims,L.G.dima);
 Mss = mdsdata('STATIC("MUT_S_S")[$1,$2]',L.G.dims,L.G.dims);
 Msx = mdsdata('STATIC("MUT_S_X")[$1,* ]',L.G.dims         );

 % LIUQE coefficients
 L = liuc(L);
 Msy = Msx(:,L.lxy);

 % LIUQE TCV measurements
 LX = liuxtcv(t,L);
 T = LX.t(2)-LX.t(1);

 % LIUQE reconstruction
 clf
 LX.Is = LX.Is*0;
 sIa = filter([1 -1],T,LX.Ia,[],2); sIa(:,1) = sIa(:,2);
 LX.Is = diag(-1./L.G.Rs)*(Msa*sIa);

 for kiter = 1:10
  LY = liut(L,LX);
%   sIa = filter([1 -1],T,LY.Ia,[],2); sIa(:,1) = sIa(:,2);
%   sIs = filter([1 -1],T,LY.Is,[],2); sIs(:,1) = sIs(:,2);
%   sJy = filter([1 -1],T,LY.Jy,[],3); sJy(:,:,1) = sJy(:,:,2);
%   Us = Msa*sIa + Mss*sIs + Msy*reshape(sJy,[],numel(LX.t));
%   LX.Is   = bsxfun(@rdivide,Us,-L.G.Rs(:));
  k = iround(LX.t,LY.t);
  norm(LX.Is(:,k)-LY.Is)
  LX.Is(:,k) = LY.Is;
  plot(LX.t,sum(LX.Is)), hold on, drawnow
 end
 
end