function varargout = RAPTOR_plot_step(out)
% RAPTOR_plot_step(out)
% params: RAPTOR parameter structure

% set some defaults (fix this later)
module_params =  struct(...
    'plotconstraints', 0,  ... % option to plot constraints or not
    'manualaxislimits',0, ... % use manually set axis limits in RAPTOR_plot_step.m
    'color',[0 0 1],... % default profile color
    'clearaxes',1,... % default to clear axes before drawing
    'plotmulti',1 ... % plot multiple lines on one axis (e.g. j|| components).
    );

%%
persistent hax

nrow = 4; ncol=3;
figname = 'RAPTOR step plot';
hax = get_persistent_hax(hax, figname, nrow, ncol, module_params.clearaxes);

%%
col = module_params.color;
plotmult = module_params.plotmulti;
rho = out.rho;

% do plots
ax=hax(1,1);
plot(ax,rho,out.psi,'color',col);
labelff(ax,'\psi','y');

ax=hax(1,2);
plot(ax,out.rhogauss,out.jauxgauss/1e6,'color',col); hold(ax,'on');
plot(ax,out.rhogauss,out.jbsgauss/1e6,'color',col,'linestyle','--')
labelff(ax,'j_{aux,bs} [MA/m^2]','y');

ax=hax(1,3);
plot(ax,out.rho,1e-6*out.jpar,'color',col);
if plotmult
    hold(ax,'on');
    plot(ax,out.rho,1e-6*out.jaux,'r');
    plot(ax,out.rho,1e-6*out.jbs,'c');
    plot(ax,out.rho,1e-6*out.joh,'m');
end
labelff(ax,'j_{par,aux,bs,oh} [MA/m^2]','y');

ax=hax(2,1);
plot(ax,rho,out.iota,'color',col)
labelff(ax,'iota','y');

ax=hax(3,1);
plot(ax,rho,out.q,'color',col)
hold(ax,'on');
plot(ax,rho,out.shear,'--k')
labelff(ax,'q, shear','y');

ax=hax(4,1);
plot(ax,out.rho,out.upl,'color',col);
hold(ax,'on');
labelff(ax,'U_{pl}','y');
axis(ax,'tight');

ax=hax(2,2);
plot(ax,out.rho,1e-3*out.te,'color',col); hold(ax,'on');
plot(ax,out.rho,1e-3*out.ti,'color',col,'linestyle','--');
labelff(ax,'T_e [keV]','y');

ax=hax(3,2);
    plot(ax,out.rhogauss,out.chiegauss,'color',col); hold(ax,'on');
    plot(ax,out.rhogauss,out.chiigauss,'color',col,'linestyle','--'); hold(ax,'on');
    labelff(ax,'\chi_e,\chi_i(--)','y');

ax=hax(4,2);
plot(ax,out.rho,1e-6*out.penet,'color',col);
if plotmult
    hold(ax,'on');
    plot(ax,out.rho,1e-6*out.pein,'r');
    plot(ax,out.rho,-1e-6*out.peout,'--');
end
labelff(ax,'P_{dens net, in, out(--)} [MW m^{-3}]','y');

ax=hax(4,3);
plot(ax,out.rho,1e-6*out.Ip,'color',col);
if plotmult
    plot(ax,out.rho,1e-6*out.Iaux,'k');
    plot(ax,out.rho,1e-6*out.Ibs,'c');
    plot(ax,out.rho,1e-6*out.Ioh,'r');
end
labelff(ax,'I_p(\rho) [MA]','y');
    
ax = hax(2,3);
plot(ax,out.rho,1e-19*out.ne,'color',col); hold(ax,'on');
plot(ax,out.rho,1e-19*out.ni,'color',col,'linestyle','--');
labelff(ax,'n_e,n_i [10^{19}m^{-3}]','y');

ax = hax(3,3);
plot(ax,out.rhogauss,out.dnegauss,'color',col); hold(ax,'on');
plot(ax,out.rhogauss,out.vnegauss,'color',col,'linestyle','--'); hold(ax,'on');
labelff(ax,'D_e,V_e(--)','y');

% grid off
for ii=1:numel(hax)
    grid(hax(ii),'off');
end

% time
title(hax(1,1),sprintf('t = %g',out.time));
drawnow
return