function [fig,storeData,timesSec]= PlotTC(dat,plotl,settings)
% Plots either all regions with all designs and all regressors, or only the
% selected ones (depending on plotlist.ploteverything)

% check whether we are specifying how many rows and columns we want
if isfield(plotl,'DesRegionsRegsFigure') && size(plotl.DesRegionsRegsFigure,2)>4
    specifyRowsCols =1;
else
    specifyRowsCols=0;
end

switch plotl.plotWhat
    case 'SortByRegion' % if we want to plot all regressors and regions from selected designs - each new figure is a new region
        c = 0;
        nregions = length(dat.design(1).out.regions);
        ndes     = length(plotl.designList);
        
        % now create plotl.DesRegionsRegsFigure
        plotl.DesRegionsRegsFigure=[];
        cfig=0;
        for ides = 1:ndes
            cdes = plotl.designList(ides);
            % Make a combination of all regressors and regions
            nregs = length(dat.design(cdes).out.labels);
            t=allcomb(1:nregions,1:nregs); % combination of regions and designs
            t = [t nan(length(t),1)];
            for it=1:length(t)
                if t(it,2) == 1 % if new regressor, make new figure
                    cfig = cfig+1;
                end
                t(it,3) = cfig;
            end
            plotl.DesRegionsRegsFigure = [plotl.DesRegionsRegsFigure; [cdes*ones(length(t),1) t]]; 
        end
    
    case 'SortByRegressor' % again all regions and regressors from selected design, each new figure is a new regressor
        c = 0;
        nregions = length(dat.design(1).out.regions);
        ndes     = length(plotl.designList);
        
        % now create plotl.DesRegionsRegsFigure
        plotl.DesRegionsRegsFigure=[];
        cfig=0;
        for ides = 1:ndes
            cdes = plotl.designList(ides);
            % Make a combination of all regressors and regions
            nregs = length(dat.design(cdes).out.labels);
            t=allcomb(1:nregions,1:nregs); % combination of regions and designs
            t= sortrows(t,2); % sort by regressors
            t = [t nan(length(t),1)];
            for it=1:length(t)
                if t(it,1) == 1 % if new regressor, make new figure
                    cfig = cfig+1;
                end
                t(it,3) = cfig;
            end
            plotl.DesRegionsRegsFigure = [plotl.DesRegionsRegsFigure; [cdes*ones(length(t),1) t]]; 
            
        end
        
    case {'TimeCoursePaperFigures','PPI'}
end

% rename for easier indexing
plots = plotl.DesRegionsRegsFigure;

% plot settings
if ~isfield(plotl,'Colors') || isempty(plotl.Colors)
    colours= repmat([0 255 255],size(plotl.DesRegionsRegsFigure,1),1); % colors in case we have more lines on the same plot
else
    colours = plotl.Colors;
end
colours= colours./255;
Liwi=2;
lineform={'-';'-.';'--';':'};


% Now plot the actual data
figcount=1;
for cp = 1:size(plots,1)
    % check whether we need to make a new figure
    
    if cp==1||plots(cp,4)~=plots(cp-1,4)
        
        thisplot=0;
        switch plotl.plotWhat
            case 'SortByRegion' 
                fig(figcount)=figure('color','w','name',[dat.design(plots(cp,1)).out.regions{plots(cp,2)} ' Design' num2str(plots(cp,1))]);
                figcount=figcount+1;
            case 'SortByRegressor'
                fig(figcount)=figure('color','w','name',[dat.design(plots(cp,1)).out.labels{plots(cp,3)} ' Design' num2str(plots(cp,1))]);
                figcount=figcount+1;
            case 'TimeCoursePaperFigures'
                fig(figcount)=figure('color','w');
                figcount=figcount+1;
            case 'PPI'
                fig(figcount)=figure('color','w','name',['Target' dat.design(plots(cp,1)).out.ppiInfo.targetRegion])
                figcount=figcount+1;
        end
        
        % check how many subplots we will need per figure
        nplots = numel(find(plots(:,4)==plots(cp,4)));
        
        if specifyRowsCols==1
            xnum = plotl.DesRegionsRegsFigure(cp,5);
            ynum = plotl.DesRegionsRegsFigure(cp,6);
        else
            if nplots<=5
                xnum=1;ynum=nplots;
            elseif nplots<=10
                xnum=2;ynum=5; 
            elseif nplots<=15               
                xnum=3;ynum=5;
            else
                xnum=5;ynum=7;
            end
        end
            
    end
    % select the subplot
    if isfield(plotl,'DesRegionsRegsFigure') && size(plotl.DesRegionsRegsFigure,2)>6
        subplot(xnum,ynum,plotl.DesRegionsRegsFigure(cp,7));hold on;
    else
        thisplot=thisplot+1;
    subplot(xnum,ynum,thisplot)
    end
    % get the correct data + labels
    pdata = squeeze(dat.design(plots(cp,1)).out.betas(:,plots(cp,2),:,plots(cp,3)));
    % Store the data so that we can print .txt files
    storeData{cp}=pdata;
    
    % Plot the data
    plot(mean(pdata),'Color',colours(cp,:),'LineStyle', lineform{1}, 'Linewidth',Liwi);
    
    hold on
    
    % Add shaded errorbars
    lst=mean(pdata(:,1))-std(pdata(:,1))/sqrt(size(pdata,1));
    ubou=mean(pdata(:,:))+std(pdata(:,:))/sqrt(size(pdata,1));
    lbou=mean(pdata(:,:))-std(pdata(:,:))/sqrt(size(pdata,1));
    patch([1 1:size(pdata,2) size(pdata,2):-1:1],[lst ubou lbou(end:-1:1)],colours(cp,:),'facealpha',0.25,'edgecolor','none')
    
    % do the 0 line
    line([0 size(pdata,2)],[0 0],'Color','k')
    % Add a straight line for where the mean median RT is
%     meanRT=dat.design(plots(cp,1)).inp(1).meanMedianDec1RT;
%     meanRT = settings.pseudosamples/settings.window*meanRT;
%     
%     %dat.design(plots(cp,1)).inp(1).SeparateMedianDec1RT
%     meanRTSTD=settings.pseudosamples/settings.window*mean(dat.design(plots(cp,1)).inp(1).SeparateMedianDec1RTSTD);
% 
%     % Here an alternative illustration for reviewer
%     %keyboard;
%     
%     
%     if ~isfield(plotl,'ylims') || isempty(plotl.ylims)
%         line([meanRT meanRT],[min(lbou) max(ubou)],'Color','k')
%        %  errorbar(meanRT,max(ubou),meanRTSTD,'k','horizontal','Linewidth',3)
%        % errorbar(meanRT,min(ubou),meanRTSTD,'k','horizontal','Linewidth',3)
%     else
%         line([meanRT meanRT],[plotl.ylims(cp,1)+0.1*abs(plotl.ylims(cp,1)) plotl.ylims(cp,2)-0.1*abs(plotl.ylims(cp,2))],'Color','k')
%          %errorbar(meanRT,plotl.ylims(cp,1)+0.1*abs(plotl.ylims(cp,1)),meanRTSTD,'k','horizontal','Linewidth',3)
%         %errorbar(meanRT,plotl.ylims(cp,2)-0.1*abs(plotl.ylims(cp,2)),meanRTSTD,'k','horizontal','Linewidth',3)
%     end
%    
%             h=errorbar(meanRT,0,meanRTSTD,'k','horizontal','Linewidth',3,'Markersize',10);
            h.CapSize = 24;
    % add xtick label
    xlabel(['Time in seconds' ],'Fontweight','Bold','FontSize',14);
    
    switch plotl.plotWhat
        case 'SortByRegion'
            title(dat.design(plots(cp,1)).out.labels{plots(cp,3)})
            %keyboard  % maybe a mistake here with the addressing??
        case 'SortByRegressor'
            title(dat.design(1).out.regions{plots(cp,2)})
        case 'PPI'
            title(dat.design(plots(cp,1)).out.labels(plots(cp,3)));
        case 'TimeCoursePaperFigures'
            title([dat.design(1).out.regions{plots(cp,2)} ' - ' dat.design(plots(cp,1)).out.labels(plots(cp,3))])
    end
    
        set(gca,'XTick',0:(settings.pseudosamples/settings.window*2):(settings.pseudosamples))
    set(gca,'XTickLabel',0+settings.TCShift:2:settings.window+settings.TCShift,'FontSize',16)
%     set(gca,'XTick',0:(settings.pseudosamples/settings.window*2):(settings.pseudosamples))
%    set(gca,'XTickLabel',0:2:settings.window,'FontSize',16)
     
    % Store the time (in sec) related to each pseudosample
    timesSecT=0:(settings.window/settings.pseudosamples):settings.window;
    timesSec=[];
    for it = 1:size(pdata,2)
        timesSec{it}= [num2str(timesSecT(it),'%.2f'), ' sec'];
    end
    
    % Determine how long to make the epoch
    if ~isfield(plotl,'TimeLimit')
        xlim([0 settings.pseudosamples])
    else
        
        xlim([0 floor(settings.pseudosamples/settings.window*plotl.TimeLimit(cp))])
    end
    if ~isfield(plotl,'ylims')
        ylim([min(lbou)-0.1*abs(min(lbou)) max(ubou)+0.1*abs(max(ubou))])
    else
        ylim([plotl.ylims(cp,1) plotl.ylims(cp,2)])
    end
end

        
      