%MEQLPACK  Packs an array of structures as a structure of arrays
% Y = MEQLPACK(X,NOPACK) returns a structure with the fields of the array X
% concatenated along the 3rd dimension and squeezed, except those in NOPACK.
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
function Y = meqlpack(X,nopack)
if nargin<2
  nopack={};
end
if numel(X) <= 1
  Y = X;
else
  for k = fieldnames(X)'
    catval = safecat3(X,k{1});
    Y.(k{1}) = cast(catval,class(X(1).(k{1})));
    if ~ismember(k,nopack)
      Y.(k{1}) = squeeze(Y.(k{1}));
      if isvector(Y.(k{1}))
        Y.(k{1}) = reshape(Y.(k{1}),1,[]);
      end
    end
  end
end
end

function catval = safecat3(mystruct,field)
% Uses cat(3,...) for most fields except if the different time slices
% are not the same size (e.g. limited/diverted transitions)

[m,n,o] = getfieldsizes(mystruct,field);
samesize = all(diff(m)==0) && all(diff(n)==0) && all(diff(o)==0);

if samesize
  n_dims = ndims(mystruct(1).(field));
  catdim = max(3,n_dims+1); % cat in 3rd dimension or higher if necessary
  catval = cat(catdim,mystruct.(field));
else
  % find maximum size of field
  V = {mystruct.(field)}; % Cell array with all values
  t = {mystruct.t};
  % Discard slices with empty time field
  mask = ~cellfun(@isempty,t);
  V = V(mask); m = m(mask); n = n(mask); o = o(mask);
  nt = numel(V);
  if max(o)<=1 % cat in 3rd dimension is ok
    % default value
    if islogical(V{1})
      catval = false(max(m),max(n),nt); % default false for logicals
    else
      catval =   NaN(max(m),max(n),nt); % default NaNs otherwise
    end
    for it=1:nt
      catval(1:m(it),1:n(it),it) = V{it};
    end
  else % rare cases with 3D objects (e.g. aq for doublets)
    if islogical(V{1})
      catval = false(max(m),max(n),max(o),nt); % default false for logicals
    else
      catval =   NaN(max(m),max(n),max(o),nt); % default NaNs otherwise
    end
    % fill each time index with field value
    for it=1:nt
      catval(1:m(it),1:n(it),1:o(it),it) = V{it};
    end
  end
end
end

function [m,n,o] = getfieldsizes(mystruct,field)
nt = numel(mystruct);
[m,n,o] = deal(zeros(nt,1));
for ii=1:nt
  [m(ii),n(ii),o(ii)] = size(mystruct(ii).(field));
end
end
