function [J, e] = x0estTimeErrorJacForCTModel(sys, yu, Ts, x0, ...
   InterSample, FreeIndex, maxsize, sqrlam)
% Error/Jacobian calculation for initial state estimation in time domain.
% For continuous-time systems.
% Returns J, e for solving J*x0 = e.
%
% Inputs:
%  sys: state-space model (ltipack.ssdata)
%  yu: [y, u]: matrix of output/input data for one experiment.
%  Ts: data sample time
%  x0: initial guess of initial states ([] can be used)
%  InterSample: cell array of 'zoh', 'foh', 'bl' strings
%  FreeIndex: index of free elements of x0.
%  maxsize: options.Advanced.MaxSize: size of largest matrix allowed to be
%           formed (default: 250e3)
%  sqrlam: sqrtm(OutputWeight)

% Copyright 1986-2022 The MathWorks, Inc.

sys = ss(sys);
nx = order(sys); [ny, nu] = size(sys.d);
nxf = sum(FreeIndex);
Ncap = size(yu,1);
M = floor(maxsize/(nxf+nu+ny));
if isempty(x0), x0 = zeros(nx,1); end
id = sys.Delay.Input;
Shift  = 0;
if any(id<0)
   % cut out the trailing max(-id) samples because those were obtained
   % assuming the input was zero.
   Shift = floor(max(-id)/Ts);
end

% Simulate using initial guess of x0; that way J\e gives delta_x0.
Ncap2 = Ncap-Shift;
if Ncap2<=0
   ctrlMsgUtils.error('Ident:estimation:x0estNegDelStructInit')
end

if ny>1 || M<Ncap
   R1 = [];
   syspsi = ltipack.ssdata(sys.a,zeros(nx,1),sys.c,zeros(ny,1),[],0);
   for kc = 1:M:Ncap
      jj = (kc:min(Ncap,kc-1+M));
      if jj(length(jj))<Ncap
         jjz = [jj,jj(length(jj))+1];
      else
         jjz = jj;
      end
      
      psi = zeros(ny*length(jj),nxf);
      t = (0:Ts:(length(jjz)-1)*Ts)';
      %[yh,xkc] = lsim(sys,yu(jjz,ny+1:end),t,x0,InterSample);
      [yh,xkc,ug] = idpack.ltidata.hybridlsim(sys,yu(jjz,ny+1:end),t,x0,InterSample);
      x0v = xkc(end,:).';
      if numel(ug)>1
         x0 = struct('ZOH', [], 'FOH', [], 'BL', []);
         L_ = 0;
         if ~isempty(ug{1}), x0.ZOH = x0v(L_+(1:nx)); L_ = L_+nx; end 
         if ~isempty(ug{2}), x0.FOH = x0v(L_+(1:nx)); L_ = L_+nx; end
         if ~isempty(ug{3}), x0.BL  = x0v(L_+(1:nx)); end
      else
         x0 = x0v;
      end
      e = (yu(jj,1:ny)-yh(1:length(jj),:))*sqrlam;
      evec = e(:);
      kl = 1;
      
      for kx = 1:nx
         if ~FreeIndex(kx), continue; end
         if kc==1
            x0psi = zeros(nx,1); x0psi(kx) = 1;
         else
            x0psi = xpsi(end,:).';
         end
         % Compute initial condition response. ISB is not important here.
         [psitemp, xpsi] = idpack.ltidata.hybridlsim(...
            syspsi,zeros(length(jjz),1),t,x0psi,{'zoh'});
         %[~,~,~,psitemp] = timeresp(syspsi,'impulse',t,RespConfig(InitialState=x0psi));
       
         psitemp = psitemp(1:length(jj),:)*sqrlam;
         psi(:,kl) = psitemp(:);
         kl = kl+1;
      end
      if ~isempty(R1)
         if size(R1,1)<nx+1
            ctrlMsgUtils.error('Ident:estimation:x0estSmallMaxSize')
         end
         R1 = R1(1:nx+1,:);
      end
      H1 = [R1;[psi,evec] ];
      if any(id<0)
         % cut out the last max(-id) samples because those were obtained
         % assuming the input was zero.
         H1 = H1(1:min(length(jj),Ncap2),:);
      end
      R1 = triu(qr(H1));
   end
   
   nr = min([size(R1),nxf]);
   J = R1(1:nr,1:nxf); e = R1(1:nr,nxf+1);
else
   % Direct estimation without segmentation
   % Jacobian system: model with ny=1 inputs and 0 outputs
   syspsi = ltipack.ssdata(sys.a.',sys.c.',zeros(0,nx),zeros(0,ny),[],0);
   t = (0:Ts:(size(yu,1)-1)*Ts)';
   y0 = idpack.ltidata.hybridlsim(sys,yu(:,2:end),t,x0,InterSample);
   %[~,psix0] = lsim(syspsi,[Ts;zeros(Ncap,1)],t,x0*0,'zoh');
   [~,~,~,psix0] = timeresp(syspsi,'impulse',t,RespConfig);
   J = psix0(1:end-Shift,FreeIndex)*sqrlam; 
   e = (yu(1:end-Shift,1)-y0(1:end-Shift))*sqrlam;
end
