function plot_convergence_study(x,y,varargin)
%% plot_convergence_study(x,y,varargin)
% Function for plotting a convergence study. x is the iteration variable
% (for example dt or the number of grid points), y a 3 dimensional vector
% where the first dimension must be the same as the number of iterations
% (the same dimension as x), the second dimension the number of lines
% plotted on the same subplot and the third dimension the number of column
% subplot in the figure. The plot has always two rows, the first display in
% semilogx the arrays (x,y), the second row displays in loglog scale
% (x(1:end-1),abs(y(1:end-1)-y(end))), that is the error wrt the better
% approximation (x must increase of precision going from index 1 to end).
% varargin contains plotting options.
%
% Input
%   x: one dimensional array, with increasing precision as increasing index
%   y: three dimensional array (but can be given as 1 or 2 dimensional 
%      array). The first dimension must be the same as x, the second are
%      the data that must be plotted in the same subplot, the third the
%      data to be plotted in different subplot.
%   varargin: -titles: cell array of strings with the same length as the 
%               third dimension of y. Gives the titles for each subplot
%               column. Default empty.
%             -legends: cell array of strings with the same length of the 
%               second*third dimension of y. Gives the legends at each line
%               in each subplot. Default empty.
%             -plot_flag: array of logicals with the same length of the
%               second dimension of y. Enable to show or to hide lines in 
%               the subplots. Default show all lines.
%             -isfigure: if true (default), create a new figure.
%             -lineSpecifics: cell array of strings containing plotting
%               features for each line. The length must be superior to the
%               second dimension of y. Default {'xb','or','+g','dk','pm','sc'}
%             -xlabels: array of char displayed below the x-axis of each
%               column. The x-axis label is the same for all the plots. 
%               Default is ''.
%
% As an example, you can plot the output of the function convergence_study
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
              
%% Initialise subplots dimensions
nlines = size(y,2);
lines  = cell(nlines,1);
ncolumn = size(y,3);
assert(ncolumn < 10,'Too many plots. Cannot afford to make more than plot per rows.');
nstart = 201 + ncolumn*10;
nfin = nstart + 2*ncolumn - 1;
nplot = nstart:nfin;

%% Set default varargin
D.titles = {};
D.legends = {};
D.plot_flag = ones(true,nlines);
D.isfigure = true;
D.linesSpecifics = {'xb','or','+g','dk','pm','sc'};
D.xlabels = '';
D.plotMarkerSize = 6;
D.plotMarkerWidth = 0.5;
D.titleFontSize = 11;
D.titleFontWeight = 'normal';
D.titleInterpreter = 'tex';
D.legendLocation = 'best';
D.legendOrientation = 'vertical';
D.legendFontSize = 6;
D.legendNumColumns = 1;
D.legendInterpreter = 'tex';
D.xlabelFontSize = 11;
D.xlabelFontWeight = 'normal';
D.xlabelInterpreter = 'tex';
P = struct;
for k = 1:2:numel(varargin), P.(varargin{k}) = varargin{k+1}; end
for k = fieldnames(D)'
  if ~isfield(P,k{1}), P.(k{1}) = D.(k{1}); end
end

assert(iscell(P.titles),'The titles for the plot must be inserted in a cell array of strings');
assert(iscell(P.legends),'The legends for the plot must be inserted in a cell array of strings');
assert(numel(P.titles) == ncolumn || isempty(P.titles),'The number of titles must be empty or equal to the number of plots');
assert(numel(P.legends) == nlines*ncolumn || isempty(P.legends),'The number of legends must be empty or equal to the number of lines plotted inj each plot');
assert(numel(P.plot_flag) == nlines,'The number of plot_flag selection must be equal to the number of lines plotted in each plot');
assert(ischar(P.xlabels),'The labels for the x axis must be a single character array (same x axis for all convergence study)');

if ~(isequal(size(P.legends),[nlines,ncolumn]) || isempty(P.legends))
  P.legends = reshape(P.legends,[nlines,ncolumn]);
end

if P.isfigure
  figure
end

%% Create the plots
for i = 1:ncolumn
  subplot(nplot(i));
  for j = 1:nlines
    lines{j} = semilogx(x,y(:,j,i),P.linesSpecifics{j},'Visible',plot_onoff_convert(P.plot_flag(j)),'MarkerSize',P.plotMarkerSize,'LineWidth',P.plotMarkerWidth); hold on;
  end
  if ~isempty(P.titles); title(P.titles{i},'FontSize',P.titleFontSize,'FontWeight',P.titleFontWeight,'Interpreter',P.titleInterpreter); end
  if ~isempty(P.legends); legend([lines{P.plot_flag}],P.legends{P.plot_flag,i},'Location',P.legendLocation,'Orientation',P.legendOrientation,'FontSize',P.legendFontSize,'NumColumns',P.legendNumColumns,'Interpreter',P.legendInterpreter); end
  grid on;

  subplot(nplot(i)+ncolumn);
  for j = 1:nlines
    lines{j} = loglog(x(1:end-1),abs((y(1:end-1,j,i)-y(end,j,i))),P.linesSpecifics{j},'Visible',plot_onoff_convert(P.plot_flag(j)),'MarkerSize',P.plotMarkerSize,'LineWidth',P.plotMarkerWidth); hold on;
  end
  if ~isempty(P.legends); legend([lines{P.plot_flag}],P.legends{P.plot_flag,i},'Location',P.legendLocation,'Orientation',P.legendOrientation,'FontSize',P.legendFontSize,'NumColumns',P.legendNumColumns,'Interpreter',P.legendInterpreter); end
  grid on;
  if ~isempty(P.xlabels); xlabel(P.xlabels,'FontSize',P.xlabelFontSize,'FontWeight',P.xlabelFontWeight,'Interpreter',P.xlabelInterpreter); end
end
end

%% Switch to show or hide a plot line
function plot_on = plot_onoff_convert(logical_on)
  if logical_on
    plot_on = 'on';
  else
    plot_on = 'off';
  end
end