function eval_debugging_plot(caller_name,radial_coordinate_label, x, model, gauss ,  varargin)
% Routine for plotting values when computed by some eval_ function

% caller_name: string indicating name of calling function (obtain via mfilename)
% radial_coordinate = label of radial coordinate 'psiN', 'rhotor', 'rhopsiN', 'psi'
% varargin  = [data, data_label]

% ex. call
% eval_debugging_plot('figure name','rhotor', hf, hax,x,model,gauss, Te, 'Te', Ti, 'Ti', ne, 'ne', ni, 'ni');

nvarargin = numel(varargin);
if mod(nvarargin, 2) ~= 0
    error('Number of varargin not even. Varargin has to be an array of [data, data_label]');
end

nplots = floor(nvarargin/2); 

% Generate the figure and the subplot handles, if not already present
figname = sprintf('%s debugging',caller_name);
hf = findobj('name',figname);
if isempty(hf) || ~ishandle(hf) || numel(hf)>1
  delete(hf); hf = figure('name', caller_name);
end  
hax = findobj(hf,'type','axes'); % find children
if isempty(hax) || ~all(ishandle(hax))
    nrow = 3; ncolumn = ceil(nplots/nrow);
    hax = zeros(nrow*ncolumn,1);
    for ii=1:nrow*ncolumn
        hax(ii) = subplot(nrow,ncolumn,ii, 'parent', hf);
    end
end

% Compute the radial coordinate
switch radial_coordinate_label
    case 'psiN'
        radial_coordinate = eval_psiN(x,[],[],model,gauss);
    case 'rhotor'
        if gauss
            radial_coordinate = model.geom.rhogauss;
        else
            radial_coordinate = model.geom.rho;
        end
    case 'rhopsiN'
        radial_coordinate = sqrt(eval_psiN(x,[],[],model,gauss));
    case 'psi'
        radial_coordinate = eval_psi(x,[],[],model,gauss);
end

% Plot all the requested data
for ii=1:nplots
    ax = hax(ii);
    y = varargin{2*(ii-1)+1};
    label = varargin{2*(ii)};
    plot(ax,radial_coordinate, y , 'b-', radial_coordinate(1), y(1),'x', 1 , y(end) ,'x');
    ylabel(ax,label)
    xlabel(ax,radial_coordinate_label)
    xlim(ax, [-0.01 1.01]) % Shift a bit the axis to see the central value
end
drawnow

return