function RAPTOR_plot_state(xk,model,varargin)
% varargin contains v
% v required in case of imposed pedestal in any kinetic profile
if ~isempty(varargin)
  v = varargin{1};
else
  v = [];
end
chanlist = {'psi','te','ti','ni','ne','n1','n2','n3'};

nstates = number_of_states(model,chanlist);

[nr,nc] = plot_config(nstates);

jchan=0;
for ichan=1:numel(chanlist)
  mychan = chanlist{ichan};
  
  if strcmp(model.(mychan).method,'state')
    jchan = jchan+1;
    evalfun = sprintf('eval_%s(xk,[],v,model,false)',mychan);
    evalfungauss = sprintf('eval_%s(xk,[],v,model,true)',mychan);
    
    mychanvalue = eval(evalfun);
    mychanvaluegauss = eval(evalfungauss);
    
    rho = model.rgrid.rho;
    rhogauss = model.rgrid.rhogauss;

    ax = subplot(nr,nc,jchan);
    plot(ax,rho,mychanvalue,'o',rhogauss,mychanvaluegauss,'.');
    ylabel(ax,mychan); xlabel(ax,'rho');
  end
end
end

function nstates = number_of_states(model,chanlist)

nstates = 0; % init
for ichan=1:numel(chanlist)
  mychan = chanlist{ichan};
  if strcmp(model.(mychan).method,'state');
    nstates = nstates+1;
  end
end
end

function [nr,nc] = plot_config(nstates)
if nstates<=3
  nr=1; nc=nstates;
else
  nr = ceil(nstates/3);
  nc = ceil(nstates/nr);
end
end
