function hf = RAPTOR_plot_snapshots(outs,plotconfig,twin,times)
% RAPTOR_plot_snapshots(outs,plotconfig,twin,times)

def = RAPTOR_plot_GUI; % get defaults

pp = plotconfig;
pp.plots = assign_defaults(plotconfig.plots,def,twin); 

ppt.style = pp.style; 
ppp.style = pp.style; % copy
ppt.plots = {}; % init plotconfig for time
ppp.plots = {}; % init plotconfig for profiles

axt = []; axp = [];
for ii=1:numel(pp.plots)
    p = pp.plots{ii};
    if ~isempty(strfind(p.xdata,'time'))
        ppt.plots = [ppt.plots,p];
        axt = unique([axt,p.axesnr]);

    else
        ppp.plots = [ppp.plots,p];
        axp = unique([axp,p.axesnr]);
    end
end

%%
hf = figure(1);
%
% number of required axes
ncols = numel(times);
nrows = numel(axp)+1;

clf(hf);
set(hf,'units','pixels','position',[100 100 800 500])
set(hf,'paperpositionmode','auto')

haxx = multiaxes(hf,ncols,nrows,[0.08 0.05],[0.05 0.08],[0.02 0.05]);
haxp = haxx(ncols+1:end); 
haxxx = mergeaxes(haxx(1:ncols));
set(haxxx,'units','normalized','position',get(haxxx,'position') + [0 0.02, 0 -0.02])
haxt = splitaxes(haxxx,0.05,numel(axt),'h');

%% plot time plots
createplots_time(haxt,ppt,outs);
% vertical bars for times
for it = 1:numel(times)
    for iax = 1:numel(haxt);
        ax=haxt(iax);
        yl = get(ax,'ylim');
        plot(ax,times(it)*[1 1],yl,'linestyle','--','color',0.6*[1 1 1]);
    end
end
        
%% plot profile plots

% reassign axes numbering 1:(n different axes)
for ip = 1:numel(ppp.plots);
    ia = find(ppp.plots{ip}.axesnr == axp); % index of axes in unique list
    ppp.plots{ip}.axesnr = ia; % new index for axes w.r.t. haxp
end

nt = numel(times);
for ii = 1:nt
    it = find(outs{1}.time>=times(ii),1,'first'); % find time index
    haxip = haxp((ii-1)+(1:nt:numel(haxp))); % axes for this time index
    createplots_profiles(haxip,ppp,outs,it);
    
    title(haxip(1),sprintf('t=%1.1fs',times(ii)));
end

%%

return

function hp = createplots_time(hax,pp,outs)
%% plots
timeplots = []; profileplots = []; % init

np = numel(pp.plots);
for iout = 1:numel(outs)
    out = outs{iout};
    ps = pp.style(iout);
    for ip=1:np
        px = pp.plots{ip};   % this plot description structure
        ax = hax(px.axesnr); % this axis handle
        
        % x data strings
        
        xdatastr = px.xdata;
        ydatastr = px.ydata;
        
        % optional color override
        if isfield(px,'col') && ~isempty(px.col)
            col = px.col;
        else
            col = ps.col; % default style
        end
        
        % plot
        try
            hpx = plot(ax,0,NaN);
            set(hpx,'XData',eval(xdatastr),...
                'YData',px.yscale*eval(ydatastr),...
                'color',col,...
                'linestyle',px.sty);
            
            % decide what kind of plot this is
        catch err
            warning('could not plot \n Xdata:%s\n Ydata:%s\n',xdatastr,ydatastr)
            hpx = NaN; % set an invalid handle
        end
        
        hold(ax,'on');
        
        if iout==1 % do labels only on first output
            % x label
            labelff(ax,px.xlabel,'x');
            
            % y label
            hly = get(ax,'ylabel');
            ystr = get(hly,'string');
            if isempty(ystr) % generate label
                ystr = sprintf('%s(%s)',px.ylabel,px.sty);
            else % append new plot
                ystr = sprintf('%s,%s(%s)',ystr,px.ylabel,px.sty);
            end
            labelff(ax,ystr,'y');
        end
        
        % axis limits
        if isfield(px,'ylim') && ~isempty(px.ylim)
            set(ax,'ylim',px.ylim);
        end
        if isfield(px,'xlim') && ~isempty(px.xlim)
            set(ax,'xlim',px.xlim);
        end
        
        % store handles
        hp(ip,iout) = hpx;
    end
end

linkaxes(hax,'x')

return


function hp = createplots_profiles(hax,pp,outs,it)
%% plots

np = numel(pp.plots);
for iout = 1:numel(outs)
    out = outs{iout};
    ps = pp.style(iout);
    for ip=1:np
        px = pp.plots{ip};   % this plot description structure
        ax = hax(px.axesnr); % this axis handle
        
        % x data strings
        if ~isempty(strfind(px.xdata,'%i'))
            xdatastr = sprintf(px.xdata,it);
        else
            xdatastr = px.xdata;
        end
        ydatastr = sprintf(px.ydata,it); % plot first profile
        
        % optional color override
        if isfield(px,'col') && ~isempty(px.col)
            col = px.col;
        else
            col = ps.col; % default style
        end
        
        % plot
        try
            hpx = plot(ax,0,NaN);
            set(hpx,'XData',eval(xdatastr),...
                'YData',px.yscale*eval(ydatastr),...
                'color',col,...
                'linestyle',px.sty);
            % decide what kind of plot this is
        catch err
            warning('could not plot \n Xdata:%s\n Ydata:%s\n',xdatastr,ydatastr)
            hpx = NaN; % set an invalid handle
        end
        
        hold(ax,'on');
        
        if iout==1 % do labels only on first output
            % x label
            labelff(ax,px.xlabel,'x');
            
            % y label
            hly = get(ax,'ylabel');
            ystr = get(hly,'string');
            if isempty(ystr) % generate label
                ystr = sprintf('%s(%s)',px.ylabel,px.sty);
            else % append new plot
                ystr = sprintf('%s,%s(%s)',ystr,px.ylabel,px.sty);
            end
            labelff(ax,ystr,'y');
        end
        
        % axis limits
        if isfield(px,'ylim') && ~isempty(px.ylim)
            set(ax,'ylim',px.ylim);
        end
        if isfield(px,'xlim') && ~isempty(px.xlim)
            set(ax,'xlim',px.xlim);
        end
        
        % store handles
        hp(ip,iout) = hpx;
    end
end
linkaxes(hax,'x')

return


function pp_mod = assign_defaults(pp,def,twin)

pp_mod = pp;
fields = fieldnames(def);

for ip=1:numel(pp)
    
    for ii=1:numel(fields)
        field = fields{ii};
        % add default to missing field
        if ~isfield(pp{ip},field) || isempty(pp{ip}.(field))
            pp_mod{ip}.(field) = def.(field);
        end
        % add time window data to time plots
        if ~isempty(twin) && ~isempty(strfind(pp{ip}.xdata,'time'))
            pp_mod{ip}.xlim = twin;
        end
    end
end
return