function [haxa,haxb] = RAPTOR_plot_multi(outs,titstr)
% [haxa,haxb] = RAPTOR_plot_multi(outs,titlestrings)
% multiple plots comparing various RAPTOR runs,
% with contours of q, Te, etc
% Mostly useful for optimization results
%
% IN:
% outs: cell of structures with RAPTOR results from RAPTOR_out.m
% titlestrings: cell of strings with corresponding name
%
% OUT:
% haxa,haxb: handles of axes, time plots resp. profile plots

ndata = numel(outs);

%%
% set_defaults_ff('pub')
colors = get(0,'defaultaxescolororder');

figure(18);
set(gcf,'name','H-mode comparison','position',[100 100 650 650]); clf;

hax = multiaxes(ndata,5,[0.08,0.05],[0.1 0.05],[0.03 0.01])';
haxb = splitaxes(hax(end,:),0.04,3,'h');
haxa = hax(1:end-1,:);

for iax=1:numel(haxa)
    set(haxa(iax),'position',get(haxa(iax),'position')+[0 0.05 0 0],'Box','on')
    hold(haxa(iax),'on');
end

% Plot reference profiles & Actuators
for id = 1:ndata
    out = outs{id};
    col = colors(id,:);
    
    ax=haxa(1,id);
    haxa2(1,id) = axes('Position',get(ax,'Position'));
    ax2 = haxa2(1,id);
    
    hp1=plot(ax,out.time,1e-6*out.U(1,:),...
        'linestyle','-','color',col); % Ip
    hp2(1)=plot(ax2,out.time,1e-6*out.Pnb(end,:),...
        'linestyle','--','color',col); % NBI
    hold(ax2,'on');
    hp2(2) = plot(out.time,1e-6*out.Pec(end,:),...
        'linestyle','-.','color',col); % EC
    
    set(ax2,'color','none','YaxisLocation','right');
    set(ax,'ylim',[0 18]);
    set(ax2,'ylim',2*get(ax,'ylim'));
    if id==1;
        ylabel(ax,'I_p [MA]');
    elseif id==3
        ylabel(ax2,'P_{NBI},P_{EC} [MW]');
    end
    hleg = legend(ax,[hp1,hp2],'I_p','P_{NB}','P_{EC}','location','southeast');
    set(hleg,'box','off');
    
    set(ax,'ylim',[0 17],'YTick',[0:5:17],'YColor',[0 0 0],'Xticklabel','')
    set(ax2,'ylim',[0 34],'YTick',[0:10:34],'YColor',[0 0 0],'Xticklabel','')
    
    if isfield(out,'Pec')
        Peccd = bsxfun(@times,1./diff(out.rho),diff(out.Pec));
        %Peccd = diff(out.Pec)
        Peccd = [zeros(size(Peccd(1,:)));Peccd];
        
        
        ax = haxa(2,id);
        cmap=colormap('hot'); colormap(flipud(cmap(:,:)))
        contourf(ax,out.time,out.rho,Peccd,51);
        shading(ax,'flat'); hold(ax,'on');
        set(ax,'clim',[0 100e6]);
    end
    
    if id==1; ylabel(ax,'\rho'); end
    set(ax,'ytick',[0.2:0.2:1]);
    title(ax,'V''P_{ECCD}(\rho,t) contours','units','normalized','position',[0.5 0.8]);
    %%
    
    %% q contour plots
    ax = haxa(3,id); cla(ax);
    cmap=colormap('hot'); colormap(flipud(cmap(:,:)))
    contourf(ax,out.time,out.rho,out.iota,51);
    shading(ax,'flat'); hold(ax,'on');
    labvec = [1,5/4 3/2, 2, 3 4]; % q surfaces to plot
    clab = fliplr(1./labvec);
    [ccf,hcf] = contour(ax,out.time,out.rho,out.iota,clab,'k');
    aa=clabel(ccf,hcf,fliplr(1./labvec),'labelspacing',100);
    for ia=1:numel(aa); set(aa(ia),'string',['q=',num2str(1./eval(get(aa(ia),'string')))]); end
    
    set(ax,'clim',[0.0 1.5])
    if id==1; ylabel(ax,'\rho'); end
    set(ax,'ytick',[0.2:0.2:1]);
    %%
    % Te contour plots
    ax = haxa(4,id); cla(ax);
    cmap=colormap('hot'); colormap(flipud(cmap(:,:)))
    contourf(ax,out.time,out.rho,out.te/1e3,101);
    shading(ax,'flat'); hold(ax,'on');
    labvec = [2,8,10,15,20,25];
    clab = labvec;
    [ccf,hcf] = contour(ax,out.time,out.rho,out.te/1e3,clab,'k');
    aa=clabel(ccf,hcf,labvec);
    for ia=1:numel(aa); set(aa(ia),'string',['Te=',num2str(eval(get(aa(ia),'string'))),'keV']); end
    
    set(ax,'clim',[0 40])
    if id==1; ylabel(ax,'\rho'); end
    set(ax,'ytick',[0:0.2:1]);
    
    %% Profile plots
    
    % q
    ax = haxb(1);
    plot(ax,out.rho,out.q(:,end),'-','color',col)
    set(ax,'ylim',[0.5 5.5]);
    hold(ax,'on')
    plot(ax,[0 1],[1 1],'linestyle',':','color',0.6*[1 1 1]);
    labelff(ax,'q','y');
    
    
    % Upl
    ax = haxb(2);
    plot(ax,out.rho,out.Upl(:,end),'-','color',col);
    set(ax,'ylim',[-0.2 0.2]);
    labelff(ax,'U_{pl}[V]','y');
    
    % Te/ne
    ax = haxb(3);
    
    plot(ax,out.rho,1e-3*out.te(:,end)','-','color',col);
    for iax=1:numel(haxb); hold(haxb(iax),'on'); end
    plot(ax,out.rho,1e-19*out.ne(:,end)','--','color',col);
    set(ax,'ylim',[0 35]);
    labelff(ax,'T_e [keV] (-), n_e 10^{-19}[m^{-3}]','y');
    
    
    title(haxa(1,id),titstr{id});
end


%% additional plot of jcd, P_EC

% % J_aux
% iax = 1;
% axes(haxb(iax));
% hpa = patch(out.rhogauss,out.jauxgauss(:,end)/(max(out.jauxgauss(:,end))),[0.8 0.8 0.6]);
% set(hpa,'edgeColor','none');
% plot(haxb(iax),out.rho,ones(model.rgrid.nrho,1),'k:','linewidth',1)
% text(0.5,0.5,'j_{aux}');
%
% % P_aux
% iax = 2;
% axes(haxb(iax));
% hpa=patch(model.rgrid,rhogauss,out.pauxgauss(:,end)/1e5,[0.6 0.6 0.6]); set(hpa,'edgecolor','none');
% text(0.4,0.7,'P_{aux}');

% Final cosmetics

% axes locations

for iha=1:numel(haxa)
    ax=haxa(iha);
    %figannotate(ax,['(',char(96+iha),')'],'rightabove',0.9,0.87)
    set(ax,'xlim',[20 100]);
    set(ax,'box','on')
    set(ax,'layer','top')
end

for iha=1:size(haxa,2);
    xlabel(haxa(end,iha),'time [s]');
end

for ihb=1:size(haxb,2);
    xlabel(haxb(end,ihb),'\rho');
end

for ihb=1:numel(haxb)
    ax = haxb(ihb);
    %figannotate(ax,['(',char(96+numel(haxa)+ihb),')'],'rightabove',0.9,0.8)
    set(ax,'position',get(haxb(ihb),'position') + [0 0 0 0],'XTick',[0:.5:1])
    set(ax,'layer','top')
    
    set(ax,'box','on')
end

set(haxa(1:end-1,:),'XTicklabel','');
set(haxa(:,2:end),'YTicklabel','');
%set(haxb(:,2:end),'YTicklabel','');
set(haxa2(1:end-1),'Yticklabel','');

for ibx=1:numel(haxb)
    set(haxb(ibx),'Box','on');
end

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