function meqplotliu(L,LX,LY, varargin)
% Debugging plot for equilibrium reconstruction
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

ntdata = numel(LX.t);
nttrec = numel(LY.t);

colorlabel = {};

% Prepare subfigure location definition
if nargin>3
  hf = varargin{1};
else
  hf = figure('Name', 'LIUQE');
end

ax(1) = subplot(3,3, [1,4,7]);
ax(2) = subplot(3,3, 2);
ax(3) = subplot(3,3, 3);
ax(4) = subplot(3,3, 5);
ax(5) = subplot(3,3, 8);

for ii=1:numel(ax)
  hold(ax(ii),'on')
end

nB = numel(LY.FB);

%%  Plotting
contour(ax(1), L.rx, L.zx, LY.Fx, 'k--', 'DisplayName', 'Fx')
plot(ax(1),L.G.rl, L.G.zl,'k', 'LineWidth', 5, 'DisplayName', 'limiter')
plot(ax(1),L.G.rv, L.G.zv,'y', 'LineWidth', 5, 'DisplayName', 'vessel')
if nB ==1
  ii = 1;
  contour(ax(1),L.rx, L.zx, LY.Fx, [LY.FB(ii) LY.FB(ii)], 'color', 'b', 'linewidth', 2, 'DisplayName', sprintf('FB domain %d', ii))
elseif nB>1
  if LY.FB(1) == LY.FB(2)
    % Doublet case
    ii = 1;
    contour(ax(1),L.rx, L.zx, LY.Fx, [LY.FB(ii) LY.FB(ii)], 'color', 'b', 'linewidth', 2, 'DisplayName', sprintf('FB(1) = FB(2) doublet'))
    plot(ax(1), LY.rB(1), LY.zB(1), 'bp', 'DisplayName', 'X point', 'MarkerSize', 10, 'MarkerFaceColor', 'b')
  else
    % Droplet case
    ii = 1;
    contour(ax(1),L.rx, L.zx, LY.Fx, [LY.FB(ii) LY.FB(ii)], 'color', 'b', 'linewidth', 2, 'DisplayName', sprintf('droplet FB domain %d', ii))
    plot(ax(1), LY.rB(ii), LY.zB(ii), 'bp', 'DisplayName', sprintf('[rB,zB] domain %d', ii), 'MarkerSize', 10,  'MarkerFaceColor', 'b')
    ii = 2;
    contour(ax(1),L.rx, L.zx, LY.Fx, [LY.FB(ii) LY.FB(ii)], 'color', 'r', 'linewidth', 2, 'DisplayName', sprintf('droplet FB domain %d', ii))
    plot(ax(1), LY.rB(ii), LY.zB(ii), 'rp', 'DisplayName', sprintf('[rB,zB] domain %d', ii),  'MarkerSize', 10, 'MarkerFaceColor', 'r') 
  end
end
if nB == 3
  ii = 3;
  contour(ax(1),L.rx, L.zx, LY.Fx, [LY.FB(ii) LY.FB(ii)], 'color', 'g', 'linewidth', 2, 'DisplayName', sprintf('FB domain %d', ii))
  plot(ax(1), LY.rB(3), LY.zB(3), 'gp', 'DisplayName',  sprintf('[rB,zB] domain %d', ii), 'MarkerSize', 10, 'MarkerFaceColor', 'g')
end
plot(ax(1),LY.rA,LY.zA, 'o', 'DisplayName', 'zA')
ylabel(ax(1),'Z [m]')
xlabel(ax(1),'R [m]')
axis(ax(1), 'equal');
legend(ax(1), 'show')


title(ax(2), 'I_a')
plot(ax(2),LY.Ia, 'r*')
errorbar(ax(2),1:L.G.na,LX.Ia, L.P.Iaerr*L.P.wIa, '-')
xticks(ax(2),1:L.G.na);
xticklabels(ax(2),  L.G.dima)
legend(ax(2),'fit', 'data')


title(ax(3),'I_u')
plot(ax(3),LY.Iu, 'r*')
errorbar(ax(3),1:L.G.nu, LX.Iu, L.P.Iuerr*L.P.wIu(1:L.G.nu),'-')
legend(ax(3),'fit', 'data')

title(ax(4), 'B_m')
plot(ax(4),LY.Bm, 'r*')
errorbar(ax(4),1:L.G.nm,LX.Bm, L.P.Bmerr*L.P.wBm(1:L.G.nm), '-')
legend(ax(4),'fit', 'data')


LX.Ff(L.P.wFf==0) = NaN;
title(ax(5),'F_f')
plot(ax(5),LY.Ff, 'r*')
errorbar(ax(5),LX.Ff, L.P.Fferr*L.P.wFf(1:L.G.nf),'-')
legend(ax(5),'fit', 'data')


%% Add messages
str = sprintf('shot = %d \n t = %.4f \n Iptot data = %.3E \n Iptot fit = %.3E \n idml = %d \n ', LY.shot, LY.t, LX.Ip, LY.Ip, L.P.idml);
if nB >1 
  for ii=1:numel(LY.IpD)
    str = strcat(str, sprintf('\n Ip domain %d = %.3E', ii, LY.IpD(ii)) );
  end
end
annotation(hf,'textbox',[0.7 0.2 0.3 0.3],'String',str,'FitBoxToText','on');


end




