function [LX2] = meqxconvert(L,LX,L2,minimal)
% [LX2] = meqxconvert(L,LX,L2[,minimal])
%
% convert LX data structure based on L
% into LX2 structure based on L2
% L,L2 may have different dima, dimu, bfct, x grid
% If minimal=true, also remove all fields not essential for running L2.code (default=false)
%
% Only the fields Ia, Iu, ag, Fx, Iy, PpQ, TTpQ, PQ, TQ, iTQ, PpQg, TTpQg, are updated.
% In case of doublet target: IpD, bpD, IniD are populated.
% Some specific fields used only in some codes are added if they are absent.
%
% NB: This does *not* yield a consistent equilibruim structure, and should
% only be used to convert LX inputs to be ingested by the various codes.
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

if nargin<4
  minimal=false; % default
end
LX2 = LX; % copy

%% Sanity checks
assert(L2.nD>=L.nD || strcmp(L2.code,'liu'), 'output nD < input nD not valid request')

%% Number of time slices
nt = numel(LX2.t);

%% Convert dima
if ~isequal(L.G.dima,L2.G.dima)
  [lia,ia2] = ismember(L.G.dima,L2.G.dima); % check all coil currents have been picked
  LX2.Ia = zeros(L2.G.na,nt);
  LX2.Ia(ia2(lia),:) = LX.Ia(lia,:);
  coilok = (lia | ~any(LX.Ia,2)); % coil selected or has no current
  assert(all(coilok),'Current in %s was not selected but is nonzero',strjoin(L.G.dima(~coilok),','))
end

%% Convert dimu
if ~isequal(L.G.dimu,L2.G.dimu)
  LX2.Iu = L2.G.Tvu \ LX.Iv;
end

%% Convert bfct parameters
% fit profiles from source to new basis for FGS

if strcmp(L2.code, 'liu')
  LX2.ag = zeros(L2.ng, nt); % LIUQE will find them with LSQ fit. No ag constraint considered yet.
else
  refitag = ~isfield(L,'bfct') || ~isequal(L.bfct,L2.bfct) || ~isequal(L.bfp,L2.bfp) || ~isfield(LX,'ag');
  if refitag
    assert(L.nD==1 && L2.nD==1 , 're-fit of P'' nad TT'' for multiple domains not available')
    LX2.ag = meqfitprof(L2,L.pQ.^2,LX.FA,LX.FB,LX.PpQ,LX.TTpQ);
  else
    if L2.ng > L.ng % ag=0 for additional functions
      LX2.ag = zeros(L2.ng, nt);
      LX2.ag(1:L.ng,:) =  LX.ag;
    end
  end
  
  if refitag || ~isequal(L.pQ,L2.pQ)
    % Re-compute basis functions on new pQ grid or following new BFs
    FN = L2.pQ.^2;
    [PpQ,TTpQ,PQ,TQ,iTQ] = deal(zeros(L2.npq+1,L2.nD,nt)); % init
    [PpQg,TTpQg] = deal(zeros(L2.npq+1,L2.nD,L2.ng,nt));
    
    for kt = 1:nt
      [F0,F1] = meqpdomlim(LX2.FA(:,kt),LX2.FB(:,kt),L2.nD);
      [PpQ(:,:,kt),TTpQ(:,:,kt),PQ(:,:,kt),TQ(:,:,kt),iTQ(:,:,kt),PpQg(:,:,:,kt),TTpQg(:,:,:,kt)] = ...
        meqprof(L2.fPg,L2.fTg,LX2.ag(:,kt),FN,F0,F1,LX2.rBt(:,kt),L2.bfct,L2.bfp,L2.idsx,L2.smalldia);
    end
    if L2.nD == 1, PpQ  = reshape( PpQ,[L2.npq+1,nt]);
      TTpQ = reshape(TTpQ,[L2.npq+1,nt]);
      PQ   = reshape(  PQ,[L2.npq+1,nt]);
      TQ   = reshape(  TQ,[L2.npq+1,nt]);
      iTQ  = reshape( iTQ,[L2.npq+1,nt]);
    end
    PpQg  = reshape( PpQg(:,1,:,:),[L2.npq+1,L2.ng,nt]);
    TTpQg = reshape(TTpQg(:,1,:,:),[L2.npq+1,L2.ng,nt]);
    LX2 = meqlarg(LX2,PpQ,TTpQ,PQ,TQ,iTQ,PpQg,TTpQg);
  end
end

%% Reinterpolate Fx, Ia if necessary
if ~all([L2.nzx,L2.nrx]==[L.nzx,L.nrx]) && ~strcmp(L2.code,'liu')
  FFx = griddedInterpolant({L.zx,L.rx},LX.Fx(:,:,1));
  FIy = griddedInterpolant({L.zy,L.ry},LX.Iy(:,:,1));
  
  ntIy = size(LX.Iy,3); ntFx = size(LX.Fx,3);
  Fx = zeros(L2.nzx,L2.nrx,ntFx);
  Iy = zeros(L2.nzy,L2.nry,ntIy);
  
  for it=1:ntFx
    if it>1, FFx.Values = LX.Fx(:,:,it); end
    Fx(:,:,it) = FFx({L2.zx,L2.rx});
  end
  for it=1:ntIy
    if it>1, FIy.Values = LX.Iy(:,:,it); end
    Iy(:,:,it) = FIy({L2.zy,L2.ry});
  end
  Iy(isnan(Iy))=0; % grid may be too wide
  Iy = Iy.*reshape(LX.Ip./sum(reshape(Iy,L2.ny,[])),1,1,[]); % ensure same total Ip
  LX2.Fx = Fx; LX2.Iy = Iy;
end


%% Multidomain specials
% Fields that have a 'D' counterpart - integral quantities
for ifield = {'IpD','bpD','IniD'}
  if isfield(LX,ifield{:}) % Copy from input structure if available
    LX2.(ifield{:}) = zeros(L2.nD,1);
    LX2.(ifield{:})(1:L.nD) = LX.(ifield{:})(1:L.nD);
  else
    if ~isfield(LX2,ifield{:}) && isfield(LX2,ifield{:}(1:end-1)) % if field without 'D' exists
      if any(L2.P.idoublet)
        % scaling to translate desired global parameters into scalar ones for
        % a symmetric doublet
        s = 1/2;
        LX2.(ifield{:}) = zeros(L2.nD,numel(LX.t)); % init
        LX2.(ifield{:})(1:2,:) = s*repmat(LX2.(ifield{:}(1:end-1)),2,1);
      else
        LX2.(ifield{:}) = LX2.(ifield{:}(1:end-1)); % copy from global field
      end
    end
  end
end

% qA special, this is always 'per A(xis)' so no D counterpart
if isfield(LX,'qA') && size(LX.qA,1)==1 && any(L2.P.idoublet)
  qA = zeros(L.nD,numel(LX.t)); % init
  qA(1:2,:) = repmat(LX.qA,2,1); % assume same qA for all lobes
  LX2.qA = qA;
end

%% contour post-processing initialization
initW = L2.P.iterq && ...
  ( L.P.iterq==0 || ... % no previous postprocessing
  (L2.G.nW~=L.G.nW) || (L2.P.nFW ~= L.P.nFW) || ... % sizes compatibility
  (L.G.nW==0 || ~(isequal(L2.G.rW,L.G.rW) && isequal(L2.G.zW,L.G.zW) && isequal(L2.G.oW,L.G.oW)))... % gaps gap settings incompatibility
  );

initq =  L2.P.iterq && ...
  ((L.P.iterq==0) || ... % no previous postprocessing
  L2.noq ~= L.noq || L2.npq ~= L.npq); % sizes mismatch

if initW, LX2.aW = []; end
if initq, LX2.aq = []; end

%% Non-defined quantities for specific code targets
if any(strcmp(L2.code,{'fge','rzp'}))
  if ~isfield(LX2,'Ini'   ), LX2.Ini    = zeros(1,numel(LX2.t)); end % Ini total
  if ~isfield(LX2,'IniD'  ), LX2.IniD   = zeros(L2.nD,numel(LX2.t)); end % Ini total
  if ~isfield(LX2,'signeo'), LX2.signeo = repmat((1-0.9*L2.pQ')*L2.P.signeodflt,1,L2.nD,numel(LX2.t)); end
  if ~isfield(LX2,'Lp'    ), Fyp        = reshape(LX2.Fx(2:end-1,2:end-1,:), L2.ny ,numel(LX2.t))- [L2.G.Mxa(L2.lxy(:),:),L2.G.Mxu(L2.lxy(:),:)]*[LX2.Ia;LX2.Iu];...
      LX2.Lp = dot(Fyp,reshape(LX2.Iy,L2.ny ,numel(LX2.t)))./LX2.Ip.^2; end % Global plasma self inductance
  if ~isfield(LX2,'SQ'    ), LX2.SQ     = evalSQ(L2,LX2); end
  if ~isfield(LX2,'Rp'    ), LX2.Rp     = repmat( 1/LX2.SQ(end)  ,1, numel(LX2.t)); end % Global resistance of the plasma.
end

%% Keep only relevant fields - toss the rest
if minimal
  switch L2.code
    case 'liu' % Implemented only for liu at this time
      required_fields = {'t','Bm','Ff','Ia','Ft','rBt','Iu','Ip'};
      optional_fields = {'Uf','Is','shot'};
      
      % check that all required fields are there
      ireqok = contains(required_fields,fieldnames(LX2));
      assert(all(ireqok),'missing required field(s): %s',strjoin(required_fields(~ireqok)));
      
      % remove unneeded fields
      for field = fieldnames(LX2)'
        myfield = field{:};
        if ~contains([required_fields,optional_fields],myfield)
          LX2 = rmfield(LX2,field{:});
        end
      end
    otherwise
      warning('''minimal=true'' option not implemented yet for code=%s',L2.code)
  end
end

end

function SQ = evalSQ(L,LX)
% Compute SQ = int_0^Phi(pQ) (sigma/T^2) dPhi
% This is a bit cumbersome but not all loaded equilibria have Opy and all
% the FA, FB etc, so need to recompute everything to finally get SQ

SQ = zeros(L.npq+1,L.nD,numel(LX.t));
if ~any(L.icde) || all(LX.Ip)==0,return; end

for kt=1:numel(LX.t)
  LXt = meqxk(LX,kt);
  LXt.Rp = 0; % Used to allow computint SQ also for non cde0Dss cases
  
  % take from LXt
  Fx = LXt.Fx; sIp = sign(LXt.Ip); ag = LXt.ag; IpD = LXt.IpD; Iy = LXt.Iy; Fy0 = zeros(L.ny,1);
  
  % meqpdom
  [rA,zA,FA,dr2FA,dz2FA,drzFA,rX,zX,FX,~,~,~,~,~,...
    rB,zB,FB,lB,lX,Opy] = meqpdom(Fx,sIp,L.P.isaddl,L); %#ok<ASGLU>
  [F0,F1] = meqpdomlim(FA,FB,L.nD); % flux limits of various domains
  
  [IpQ,~,~,~,FtQ,Ft0Q] = meqintQ(L,F0,F1,LXt.rBt,ag,Fx,Opy,sIp);
  FtPQ = FtQ+Ft0Q; % total Plasma toroidal flux on pQ grid
  [~,~,~,TQ] = meqprof(L.fPg,L.fTg,ag,L.pQ(:).^2,F0,F1,LXt.rBt,L.bfct,L.bfp,L.idsx);
  
  for iD=1:numel(FB)
    [~,SQ(:,iD,kt)] = meqcde(L,LXt,LXt,F0,F1,Opy,Iy,Fx,Fy0,IpD(iD),iD,IpQ(:,iD),FtPQ(:,iD),TQ(:,iD));
  end
end
SQ = squeeze(SQ);
end
