function [haxa,haxb,hlegg] = RAPTOR_plot_overview_timeslices(out,varargin)
% function RAPTOR_plot_overview_timeslices(out,varargin)
% out: RAPTOR out structure
% varargin{1}: times to plot profiles
 

%% FIGURE SETUP
% time indices to plot
if isempty(varargin)
    times = linspace(1.05*out.time(1),0.95*out.time(end),4);
else
    times = varargin{1};
end

nit = numel(times);

for it=1:numel(times)
    tind(it) = find(out.time>times(it),1,'first');
end

assert(numel(tind)==nit,'invalid times')

if exist('set_defaults_ff','file')
set_defaults_ff('pub')
end

figname = 'RAPTOR overview';
hf = findobj('name',figname);

if ishandle(hf)
    figure(hf);
    clf(hf);
else
    hf = figure;
    set(hf,'name',figname,'position',[10 0 600 600]);
end

hax = multiaxes(hf,2,5,[0.05 0.08],[0.1 0.1],[0.05 0.01]);
haxb1 = splitaxes(hax([5,6]),0.02,nit,'h');
haxb2 = splitaxes(hax([7,8]),0.02,nit,'h');
haxb3 = splitaxes(hax([9,10]),0.02,nit,'h');

haxa = hax(1:2,1:2);
haxb = [haxb1;haxb2;haxb3];

for iax=1:numel(haxa);
    ax=haxa(iax);
    set(ax,'position',get(ax,'position') + [0 0.08 0 0])
end

% plot gray vertical lines first
for iax=1:numel(haxa)
    ax=haxa(iax);
    gray = [0.6 0.6 0.6];
    for it = 1:nit
        hp = plot(ax,times(it)*[1 1],[0 100],'color',gray,'linewidth',1,'linestyle','--');
        hold(ax,'on')
    end
end

%% Time traces

ax=haxa(1);
% plot Ip
plot(ax,out.time,1e-6*out.U(1,:))
labelff(ax,'I_{p} [MA]','y')

% plot powers
ax=haxa(2);

% NBI
if size(out.U,1)>1
hpnb = plot(ax,out.time,1e-6*out.Pnb(1,:),'r');
hpec = plot(ax,out.time,1e-6*out.Pec(end,:));
hleg=legend([hpnb(1),hpec(1)],'P_{NBI}','P_{EC} [MW]');
set(hleg,'location','northwest','box','off');
end

% plot Te, ne        
ax=haxa(3);clear hp;
hp(1) = plot(ax,out.time,1e-3*out.te(1,:),'-');
hp(2) = plot(ax,out.time,1e-19*out.ne(1,:),'--','color',[0.6 0.6 0.6]);
if norm(out.Q)>0
    hp(3) = plot(ax,out.time,out.Q,'-.','color','k');
end

hleg=legend(ax,hp,'T_{e0} [keV]','n_{e0} [10^{-19}m^{-3}]','Q');
set(hleg,'location','northwest','box','off');

% Powers
ax=haxa(4);

hold(ax,'on'); clear hp;
hp(1) = plot(ax,out.time,out.Pheat*1e-6,'color',[0 0.5 0]);
if norm(out.Palphatot)>0
hp(2) = plot(ax,out.time,out.Palphatot(end,:)*1e-6,'--');
end

hleg=legend(ax,hp,'P_{heating} [MW]','P_\alpha [MW]');
set(hleg,'location','northwest','box','off');

%% Profile plots of timeslices
for ii = 1:nit 
    tt = tind(ii);

    ax = haxb(1,ii);  
    
    plot(ax,out.rho,out.q(:,tt),'b');
    hold(ax,'on');
    
    plot(ax,out.rho,5*out.upl(:,tt),'k')
    hlegg(ii)=legend(ax,'q','5xU_{pl} [V]');
    set(hlegg(ii),'location','south','box','off','orient','hor')
   
    title(ax,sprintf('t=%ds',round(times(ii))));

    %%
    ax =haxb(2,ii);
    plot(ax,out.rho,1e-3*out.te(:,tt),'b-');
    hold(ax,'on')
    plot(ax,out.rho,1e-19*out.ne(:,tt),'--','color',[0.6 0.6 0.6]);
        
    hleg=legend(ax,'T_e [keV]','n_e [10^{-19}m^{-3}]');
    set(hleg,'location','northeast','box','off');

    ax = haxb(3,ii);
    xpa = [out.rho;out.rho(end);out.rho(1)];
    
    JJ = [out.joh(:,tt),out.jnb(:,tt),out.jec(:,tt),out.jbs(:,tt)];
    cJJ = cumsum(JJ,2)*1e-6;
    
    col = [1 0.4 0.4;
           0.5 0.5 0.5;
           0.5 0.5 1;
           0 1 1;
           ];
    
       for ij = [4 3 2 1]
            ypa = [cJJ(:,ij);0;0];
            hpa(ij) = patch(xpa,ypa,'r',...
                'parent',ax,'facecolor',col(ij,:));
       end
       if ii==1
            [hleg,hobj,hpl,hstr] = legend(hpa,'\Omega','NB','EC','BS');
            set(hleg,'location','east','orientation','vert','box','off');
            labelff(ax,'Cumulative j_{||}[MA/m^2]','y');
            for io = 5:8 % move legend patch to left
                oldpos = get(hobj(io),'Xdata'); % move to the left
                newpos = oldpos; newpos(1:2)=0.4;
                set(hobj(io),'Xdata',newpos);
            end            
 
       end
end

%% Axis cosmetics

for iax=1:numel(haxa)
    ax=haxa(iax);
    set(ax,'xlim',[out.time(1),out.time(end)])
end

set(haxa(1),'ylim',[0 16]);
set(haxa(2),'ylim',[0 34]);
set(haxa(3),'ylim',[0 28]);
set(haxa(4),'ylim',[0 65]);

set(haxa(:,1),'Xticklabel','');
set(haxb(1:end-1,:),'Xticklabel','')
set(haxb(1,:),'ylim',[-1 5],'ytick',[-1:1:5],'ygrid','on');
set(haxb(2,:),'ylim',[0 35]);
set(haxb(3,:),'ylim',[0 1.9]);
set(haxb(:,2:end),'yticklabel','');
set(haxb,'box','on');
%%
hxl = xlabel(haxa(1,2),'time [s]');
set(hxl,'units','normalized'); set(hxl,'pos',[0.5,-0.2],'horizontalalignment','center');

hxl = xlabel(haxa(2,2),'time [s]');
set(hxl,'units','normalized'); set(hxl,'pos',[0.5,-0.2],'horizontalalignment','center'); 

for iax = 1:size(haxb,2)
xlabel(haxb(3,iax),'\rho','pos',[0.5,-0.25],'horizontalalignment','center');
end

set(gcf,'paperpositionmode','auto');

return

