function hax=RAPTOR_plotNF(varargin)
% RAPTOR_plot({out},twindow,plottimes,{colorlist},{symbollist},{legends})

if nargin>0
        out = varargin{1};
else
    error
end

% check out is cell
nout = numel(out);
if ~iscell(out) && nout==1
        out = {out};
end

if nargin>1 && ~isempty(varargin{2})
    twindow = varargin{2};
end

if nargin>2 && ~isempty(varargin{3})
    plottimes = varargin{3};
else
    plottimes = linspace(out{1}.time(1),out{1}.time(end),3);
end

if nargin>3 && ~isempty(varargin{4})
    colorlist = varargin{4};
else
    colorlist = {'k','r','b','g'};
end

if nargin>4 && ~isempty(varargin{5})
    symbollist = varargin{1};
else
    symbollist = {'-','--','-.',':'};
end

if nargin>5 && ~isempty(varargin{5})
    legends = varargin{5};
else
    legends= {'Results 1','Results 2','Results 3','Results 4'};
end

assert(numel(colorlist) >= nout,'legends has wrong size');
assert(numel(symbollist) >= nout,'legends has wrong size');
assert(numel(legends) >= nout,'legends has wrong size');


%% start plot
figure(111)
hax = multiaxes(4,4,[0.05,0.05],[0.05,0.08],[0.05,0.02]);
set(gcf,'position',[50 50 700 500]);

%% loop for plotting each simulation result in one figure
for ii=1:nout
    
    % plot inputs
    for iu=1:size(out{ii}.U,1)-1
        plot(hax(1,iu), out{ii}.time,out{ii}.U(iu,:)/1e6,symbollist{ii},'color',colorlist{ii},'linewidth',2)
        hold(hax(1,iu),'on')
        if iu==1
            labelff(hax(1,iu),'I_p [MA]','y')
        elseif iu==2            
            labelff(hax(1,iu),['P_{EC}(\rho=0.2) [MW]'],'y')
            elseif iu==3            
            labelff(hax(1,iu),['P_{EC}(\rho=0.4) [MW]'],'y')
            elseif iu==4            
            labelff(hax(1,iu),['P_{EC}(\rho=0.55) [MW]'],'y')
        end
        set(hax(1,iu),'box','on')
        xlim(hax(1,iu),twindow)
        if iu<4
        set(hax(1,iu),'XTickLabel',[])
        end
    end
    title(hax(1,1),'Input evolutions')
    legend(hax(1,1),legends,'Location','SouthEast')
    xlabel(hax(1,4),'t [s]')
    
    
    % iota evolutions
    plot(hax(2,1),out{ii}.time,out{ii}.iota(1,:),symbollist{ii},'color',colorlist{ii},'linewidth',2)
    hold(hax(2,1),'on')
    set(hax(2,1),'box','on')
    xlim(hax(2,1),twindow)
    labelff(hax(2,1),'\iota (\rho=0)','y')
    title(hax(2,1),'\iota and q evolutions')
    set(hax(2,1),'XTickLabel',[])
    
    
    plot(hax(2,2),out{ii}.time,out{ii}.iota(7,:),symbollist{ii},'color',colorlist{ii},'linewidth',2)
    hold(hax(2,2),'on')
    set(hax(2,2),'box','on')
    xlim(hax(2,2),twindow)
    labelff(hax(2,2),'\iota (\rho=0.2)','y')
    set(hax(2,2),'XTickLabel',[])
    
    plot(hax(2,3),out{ii}.time,out{ii}.iota(13,:),symbollist{ii},'color',colorlist{ii},'linewidth',2)
    hold(hax(2,3),'on')
    set(hax(2,3),'box','on')
    xlim(hax(2,3),twindow)
    labelff(hax(2,3),'\iota (\rho=0.4)','y')
    set(hax(2,3),'XTickLabel',[])
    
    % q95
    plot(hax(2,4),out{ii}.time,out{ii}.q95,symbollist{ii},'color',colorlist{ii},'linewidth',2)
    hold(hax(2,4),'on')
    set(hax(2,4),'box','on')
    xlim(hax(2,4),twindow)
    labelff(hax(2,4),'q_{95}','y')
    
    xlabel(hax(2,4),'t [s]')
    
    
    % Te0,Ti0
    plot(hax(3,1),out{ii}.time,out{ii}.te(1,:)/1e3,'-','color',colorlist{ii},'linewidth',2)
    hold(hax(3,1),'on')
    plot(hax(3,1),out{ii}.time,out{ii}.ti(1,:)/1e3,':','color',colorlist{ii},'linewidth',2)
    set(hax(3,1),'box','on')
    xlim(hax(3,1),twindow)
    labelff(hax(3,1),'T_{e0} (-), T_{i0} (:) [keV]','y')
    title(hax(3,1),'Several quantity evolutions')
    set(hax(3,1),'XTickLabel',[])
    
    % ne0,ni0
    plot(hax(3,2),out{ii}.time,out{ii}.ne(1,:)/1e19,'-','color',colorlist{ii},'linewidth',2)
    hold(hax(3,2),'on')
    plot(hax(3,2),out{ii}.time,out{ii}.ni(1,:)/1e19,':','color',colorlist{ii},'linewidth',2)
    set(hax(3,2),'box','on')
    xlim(hax(3,2),twindow)
    labelff(hax(3,2),'n_{e0} (-), n_{i0} (:) [10^{19} m^{-3}]','y')
    set(hax(3,2),'XTickLabel',[])
    
    % Wtot
    plot(hax(3,3),out{ii}.time,out{ii}.Ptotin/1e6,symbollist{ii},'color',colorlist{ii},'linewidth',2) 
    hold(hax(3,3),'on')
     plot(hax(3,3),out{ii}.time,out{ii}.U(size(out{ii}.U,1),:)/1e6,':','color',colorlist{ii},'linewidth',2)
    set(hax(3,3),'box','on')
    xlim(hax(3,3),twindow)
    labelff(hax(3,3),'P_{totin} [MW]','y')
    set(hax(3,3),'XTickLabel',[])
    
    % Ptot
    plot(hax(3,4),out{ii}.time,out{ii}.H98,symbollist{ii},'color',colorlist{ii},'linewidth',2)
    hold(hax(3,4),'on')
    set(hax(3,4),'box','on')
    xlim(hax(3,4),twindow)
    labelff(hax(3,4),'H_{98}','y')
    xlabel(hax(3,4),'t [s]')
    % profiles at certain timestamp
%     tstamp=out{ii}.time(end);
    timestamp=find(out{ii}.time>=plottimes,1,'first');
    
    % Upl
    plot(hax(4,1),out{ii}.rho,out{ii}.Upl(:,timestamp),symbollist{ii},'color',colorlist{ii},'linewidth',2)
    hold(hax(4,1),'on')
    set(hax(4,1),'box','on')
    xlim(hax(4,1),[0 1])
    labelff(hax(4,1),'U_{pl}','y')
    title(hax(4,1),['Profile quantities @t=' num2str(plottimes) 's'])
    set(hax(4,1),'XTickLabel',[])
    % iota
    plot(hax(4,2),out{ii}.rho,out{ii}.iota(:,timestamp),symbollist{ii},'color',colorlist{ii},'linewidth',2)
    hold(hax(4,2),'on')
    set(hax(4,2),'box','on')
    xlim(hax(4,2),[0 1])
    labelff(hax(4,2),'\iota','y')
    set(hax(4,2),'XTickLabel',[])
    
    % q
    plot(hax(4,3),out{ii}.rho,out{ii}.q(:,timestamp),symbollist{ii},'color',colorlist{ii},'linewidth',2)
    hold(hax(4,3),'on')
    set(hax(4,3),'box','on')
    xlim(hax(4,3),[0 1])
    labelff(hax(4,3),'q','y')
    set(hax(4,3),'XTickLabel',[])
    
    % Te,Ti
    plot(hax(4,4),out{ii}.rho,out{ii}.te(:,timestamp)/1e3,'-','color',colorlist{ii},'linewidth',2)
    hold(hax(4,4),'on')
    plot(hax(4,4),out{ii}.rho,out{ii}.ti(:,timestamp)/1e3,':','color',colorlist{ii},'linewidth',2)
    set(hax(4,4),'box','on')
    xlim(hax(4,4),[0 1])
    labelff(hax(4,4),'T_{e} (-), T_{i} (:) [keV]','y')
    xlabel(hax(4,4),'\rho')    

end

annotation('textbox',[0.62 0.255 0.1 0.1],'String',{'P_{NBI} [MW]'},'FitBoxToText','on','EdgeColor','none');

return
