function varargout=dopEventStack(eventData,eventMarkers,varargin)
% dopOSCCI Mark II: dopEventStack ~ 18-Jan-2012 (last edit)
%
% create an event stack of data based upon epoch settings
%
% stack will be an 3D matrix with:
% - x = event
% - y = sample
% - z = left-right data
%
% if data isn't available, mirror* available data to fill space and
% recommend that the trial not be used.
%
% * mirror
% when data missing from sample -50 to 1, the reverse data from sample 1 to
% 50 to fill the gap.
%
% * consider additional methods for filling:
% - leave (use NaN) and make later routines ignore this
% - repeat = repeat closest sample value to fill gap
% > not yet programmed 7-9-2010
%
% also - fill routines have not been checked 7-9-2010
%

fprintf('\n%s\n\t%s',['Running ',mfilename],...
    'Stacking Event Data:');
try
    inputs.varargin=varargin;
    inputs.defaults=struct('epoch',[-10 10],... % epoch around event marker in seconds
        'sample',25,... % sample rate of input (to figure out timing)
        'eventLabels','mark1',...); % names of event channels
        'evChs',0,... % doppler external channel numbers
        'fillMethod','mirror',...% method to fill empty data with
        'record',[]); % record sample rate - used when raw data fed in
    
    st=setGetInputsStruct(inputs);
    
    st.translation='stack - event stack';
    st.structAbbreviations='st';
    st.data=eventData; clear eventData
    st.evMks=eventMarkers; clear eventMarkers
    
    %% downsample raw data
    if ~isempty(st.record)
        if exist('downsample','file') == 2
            st.data = downsample(st.data,st.record/st.sample);
        else
            st.data = dopDownsample(st.data,st.record/st.sample);
        end
    end
    % make sure the data is layed out the way we want
    st.ds=size(st.data); % data size
    if min(st.ds)~=st.ds(2); % want this to be a column
        st.data=st.data'; end % flip it around
    
    st.epochSamples=st.epoch/(st.sample^-1);
    st.epochLength=diff(st.epochSamples,[],2);
    
    % for each set of markers, create extra stack - stacks may be different
    % shape so might as well be separate variables in the structure
    st.out.evNames=st.evMks.evNames;% fieldnames(st.evMks); % event names (output 29-10-10)
    st.nEvs=length(st.out.evNames); % number of event channels
    
    for i=1:st.nEvs
        st.tmp.ev=i; % event channel number
        st.tmp.evn=st.out.evNames{st.tmp.ev}; % event name
        
        st.tmp.data=st.data; % create temp data incase multiple channel markers require different adjustments
        
        st.tmp.mks=st.evMks.(st.tmp.evn).marks; % pull the markers out of the cell array
        st.tmp.nmks=length(st.tmp.mks);
        st.tmp.limits=[1 length(st.tmp.data)]; % upper and lower limits of available data
        
        st.tmp.evEpochMks=[st.epochSamples(st.tmp.ev,1)+st.tmp.mks' st.epochSamples(st.tmp.ev,2)+st.tmp.mks'];
        
%         st.tmp.evName=st.out.evNames{st.tmp.ev};
%         if iscell(st.eventLabels)
%             st.tmp.evName=st.eventLabels{st.tmp.ev};
%         end
                    
        %% check data limits:
        % check if require data outside of available limits
        st.tmp.startData=[]; % data concatenated to above available data (before)
        st.tmp.finishData=[]; % data concatenated to below available data (after)
        if sum(sum(st.tmp.evEpochMks<st.tmp.limits(1)))
            switch st.fillMethod
                case 'mirror'
                    st.tmp.startData=flipud(st.tmp.data(1:-st.tmp.evEpochMks(1,1)+1,:)); % +1 to include 0
            end
        end
        if sum(sum(st.tmp.evEpochMks>st.tmp.limits(2)))
            st.tmp.finishAmt=st.tmp.evEpochMks(end,end)+1-st.tmp.limits(2); % samples required +1 to include 0 (sort of 0)
            switch st.fillMethod
                case 'mirror'
                    st.tmp.finishData=flipud(st.tmp.data(st.tmp.limits(2)-st.tmp.finishAmt:st.tmp.limits(2),:));
            end
        end
        
        %% adjust if necessary
        if ~isempty(st.tmp.startData) || ~isempty(st.tmp.finishData)
            st.tmp.adjLabels={'start','finish'};
            for j=1:length(st.tmp.adjLabels)
                st.tmp.al=st.tmp.adjLabels{j}; % adjust label
                if ~isempty(st.tmp.([st.tmp.al,'Data']))
                    st.tmp.([st.tmp.al,'Pct'])=...
                        length(st.tmp.([st.tmp.al,'Data']))/st.epochLength(st.tmp.ev); % percentage of epoch affected
                    switch st.tmp.al
                        case 'start'
                            st.tmp.data=vertcat(st.tmp.startData,st.tmp.data);
                            % & adjust marker values
                            st.tmp.evEpochMks=st.tmp.evEpochMks+length(st.tmp.startData);
                            % everything shifted down, minimum negative value should become
                            % 1 - with any luck...
                        case 'finish'
                            st.tmp.data=vertcat(st.tmp.data,st.tmp.finishData);
                    end
                    %% report
                    % no indication here whether trials with be excluded or
                    % not - save this for a screening routine using the
                    % immediately below variable (if it exists)
                    % = dop.stack.short (start and/or finish)
                    st.out.short.(st.tmp.evn).(st.tmp.al)=st.tmp.([st.tmp.al,'Pct']);
                    st.out.report.(st.tmp.evn).(st.tmp.al)=...
                        sprintf('%s\n\t%s\n\n',...
                        ['Event channel ',num2str(st.tmp.ev),' (''',st.tmp.evn,...
                        ''' external channel ',num2str(st.evChs(st.tmp.ev)),'):'],...
                        [st.tmp.al,' data short by ',...
                        num2str(st.tmp.([st.tmp.al,'Pct'])),' epochs (',...
                        num2str(length(st.tmp.([st.tmp.al,'Data']))),' samples).']);
                    disp(st.out.report.(st.tmp.evn).(st.tmp.al));
                end
            end
        end % end of adjustment check
        %% create stacks
        for j=1:st.tmp.nmks % loop this
            st.out.(st.tmp.evn)(j,:,1:2)=...
                st.tmp.data(st.tmp.evEpochMks(j,1):st.tmp.evEpochMks(j,2),:);
%             disp(st.tmp.data(st.tmp.evEpochMks(j,1):st.tmp.evEpochMks(j,1)+5,:))
        end
    end
    %% set ouputs
    varargout{1}=st.out;
%     varargout{2}=st.evMks; % may have to modify these
    fprintf('\t%s\n','Done');
catch err
    %% catch dopOSCCI error
    % 
save(dopOSCCIdebug(mfilename));rethrow(err); 
%     dopOSCCIcrash(err,st,mfilename);
end

clear st inputs

