function LX = meqinterp(LY,t,method)
%MEQINTERP Interpolate LY fields on to time basis t.
% LX = meqinterp(LY,t,method);
% method can be 'linear' or 'nearest'.
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

if nargin<3
  method = 'linear';
end

switch method
  case 'nearest'
    LX = meqxk(LY,iround(LY.t,t));
    LX.t = t;
  otherwise
    LX = meqinterp1(LY,t,method);
end

end

function LX = meqinterp1(LY,t,method)
  assert( max(t)<=max(LY.t) && min(t)>=min(LY.t),...
     'Range t [%2.2f,%2.2f] exceeds range of LY.t [%2.2f,%2.2f]',min(t),max(t),min(LY.t),max(LY.t));
  for field=fieldnames(LY)'
    fieldLY = LY.(field{:});
    timedim = find(size(fieldLY)==numel(LY.t),1,'last');
    if isempty(timedim)
      LX.(field{:}) = fieldLY; % pass direct
    else % interpolate
      % permute since time dimension must be first
      N=ndims(fieldLY);
      p=[timedim,1:timedim-1,timedim+1:N];
      LX.(field{:}) = ipermute(interp1(LY.t',permute(fieldLY,p),t',method),p);
    end
  end
end
