function varargout=dopGraph3lateralityStack(varargin)
% dopOSCCI Mark II: dopGraph3lateralityStack~ 22-9-2010 (last edit)
%
% make laterality graph
%
% accept dop.stack data structure
% OR 'data' input = epochs x epoch length x 2 (left right) matrix
%
% Stack variation receives full stack of multiple input channels but I
% think it'd be useful to have a simpler version that just takes the data -
% all the setup would be done outside in dopGraph. (11-10-2010)
%
fprintf('\n%s\n\t%s',['Running ',mfilename],...
    'Making Graph:');
try
    inputs.varargin=varargin;
    inputs.turnOn={'test'}; % makes the figure visible automatically
    inputs.defaults=struct('data',[],...); send in a stack of data
        'stack',[],... dop.stack structure with event marker labels of data
        'sampleRate',25,...); % down sampleRate in hertz to figure out timing
        'indName','Bruce',... % individual data file identification
        'graphName','Laterality',...); % graph name
        'outPos',[0 .3 1 .7],...); % outer position of figure in normalized units
        'actRej',[70 130],...); % activation rejection range
        'mk',[],...); % dop.mk structure - event names + raw and matched markers
        'epoch',[],... % lower and upper epoch values for events
        'base',[],...); % lower and upper epoch values for events
        'poi',[],...); % period of interest
        'zeroPt',0,...% event marker placement
        'yExt',.6,...); % y-axis extention beyond minMax of data
        'li',[],...); % laterality index structure
        'actWin',2,...); % activation window duration
        'actWinExt',1.2); % .2 portion extention on activation window area
    gl=setGetInputsStruct(inputs);
    
    gl.translation='graph laterality';
    gl.structAbbreviations='gl';
    gl.sampleSec=gl.sampleRate^-1;
    gl.subNames={'Left & Right','Difference (Left - Right)'};
    %% check for events:
    gl.ev.names='event1';
    if ~isempty(gl.stack);
        gl.ev.names=fieldnames(gl.stack);
    end
    gl.ev.n=length(gl.ev.names);
    for j=1:gl.ev.n
        %% setup figure
        gl.ev.name=gl.ev.names{j};
        gl.h(j)=figure('NumberTitle','off','Name',[gl.graphName,' ',gl.ev.name,': ',gl.indName],...
            'Units','normalized','OuterPosition',gl.outPos,'Visible','Off');
        gl.hNames{j}=[gl.indName,'_',gl.ev.name];
        if ~isempty(gl.stack);
            gl.data=gl.stack.(gl.ev.name);
            if isfield(gl.zeroPt,gl.ev.name)
                gl.zeroPt=gl.zeroPt.(gl.ev.name);
            end
            if ~isempty(gl.li)
                gl.ev.li=gl.li.(gl.ev.name).scrnEpochs.all;
                %               note: [1x444 char]
                %               LI: 2.0652
                %         LIepochs: 10
                %     LIpeakSample: 206
                %        LIpeakSec: 8.2400
                %            LIsem: 0.1877
                %            LIsig: 0
                %              LIp: 0.1523
                %          LIcohen: 0.4945
            end
        end
        for jj=1:2
            gl.subplot=jj;
            subplot(2,1,jj);
            gl.areaHandles=[];
            switch jj
                case 1
                    %% left right data = subplot 1
                    gl.lr.translation='left right';
                    gl.lr.labels={'L (left)','R (right)'}; % TX data
                    gl.lr.colours={'b','r'};
                    gl.xdata=((1:length(gl.data))-gl.zeroPt)*gl.sampleSec;
                    
                    for i=1:size(gl.data,3)
                        gl.ydata(i,:)=mean(gl.data(:,:,i));
                        plot(gl.xdata,gl.ydata(i,:),'color',gl.lr.colours{i},'DisplayName',gl.lr.labels{i})
                        gl.minMax(i,:)=[min(gl.ydata(i,:)) max(gl.ydata(i,:))];
                        if i==1; hold; end
                    end
                case 2
                    gl.ydata=diff(mean(gl.data),1,3);
                    plot(gl.xdata,gl.ydata,'color','k','DisplayName','L - R'); hold;
                    gl.minMax=[min(gl.ydata)...
                        max(gl.ydata)];
            end
            gl.yLim=[min(min(gl.minMax))-max(max(abs(gl.minMax)))*gl.yExt ...
                max(max(gl.minMax))+max(max(abs(gl.minMax)))*gl.yExt];
            %% baseline
            if ~isempty(gl.base)
                gl.areaHandles(end+1)=area(gl.base(j,:),[gl.yLim(1) gl.yLim(1)],...
                    gl.yLim(2),'FaceColor',[.9 .9 .9],...
                    'DisplayName','base');
            end
            %% period of interest
            if ~isempty(gl.poi)
                gl.areaHandles(end+1)=area(gl.poi(j,:),[gl.yLim(1) gl.yLim(1)],...
                    gl.yLim(2),'FaceColor',[.1 .8 .3],...
                    'DisplayName','POI');
            end
            %% peak difference
            if isfield(gl,'ev') &&isfield(gl.ev,'li')
                gl.liSel=(gl.ev.li.LIpeakSample-(1/gl.sampleSec)...
                    :gl.ev.li.LIpeakSample+(1/gl.sampleSec))+gl.zeroPt;
                gl.liXdata=gl.xdata(gl.liSel);
                switch jj
                    case 1
                        gl.liYdata(1,:)=ones(1,length(gl.liXdata))...
                            *(max(max(abs(gl.ydata(:,gl.liSel)))*gl.actWinExt)); %gl.ydata(:,gl.liSel);
                        gl.areaHandles(end+1)=area(...
                            gl.liXdata,... % x values
                            gl.liYdata(1,:),... % y values
                            -(max(max(abs(gl.ydata(:,gl.liSel))))*gl.actWinExt));%min(min(abs(gl.ydata(:,gl.liSel))))...
                    case 2
                        gl.liYdata=gl.ydata(gl.liSel);
                        gl.areaHandles(end+1)=area(...
                            gl.liXdata,gl.liYdata,0);
                end
                set(gl.areaHandles(end),'FaceColor','y',...
                    'DisplayName','act. win.','EdgeColor','r',...
                    'LineStyle',':');
                %% > peak point
                plot(ones(1,2)*gl.ev.li.LIpeakSec,gl.yLim,...
                    'color','c','LineStyle','-','LineWidth',2,...
                    'DisplayName','peak dif.');
            end
            %% legend
            gl.lgd=legend('location','northWest','orientation','horizontal');
            if jj==1
                % just left and right for top section
                gl.lgdString=get(gl.lgd,'String');
                set(gl.lgd,'String',gl.lgdString([1 2]));
            end
            %             set(gl.lgd,'Box','off');
            %% set axes
            gl.axes=get(gl.h,'CurrentAxes');
            set(get(gl.axes,'Title'),'string',[gl.subNames{jj},': ',gl.indName]);
            
            set(gl.axes,'XLim',[min(gl.xdata) max(gl.xdata)]);
            set(gl.axes,'YLim',gl.yLim);
            %             set(gl.axes,'XTick',(1:size(gl.data,1)*(gl.sampleSec/60))*60)
            set(get(gl.axes,'YLabel'),'String','Doppler Velocity (cm/s)');
            if jj==2
                set(get(gl.axes,'XLabel'),'String',...
                    ['Epoch Time in seconds ( 0 = ',gl.ev.name,')']);
            end
            %% zero points
            plot([0 0],gl.yLim,'color','g',...
                'LineWidth',2,'DisplayName',gl.ev.name);
            plot(gl.epoch,[0 0],'color','k','LineStyle',':',...
                'LineWidth',1,'DisplayName','y = 0');
            %% place areas under
            if ~isempty(gl.areaHandles) %exist(gl.baseH,'var')
                for i=1:length(gl.areaHandles);
                    uistack(gl.areaHandles(length(gl.areaHandles)+1-i),'bottom')
                end
            end
        end
    end
    %% test
    if gl.test; set(gl.h,'Visible','on'); end
    %         fprintf('%s\n','Test procedure on - throwing error for debugging');
    %         plang; % should generate an error
    %     end
    
    %% set ouputs
    %    gl.out.note='li = laterality index output';
    varargout{1}=gl.h;
    varargout{2}=gl.hNames;
    fprintf('\t%s\n','Done');
catch err
    fprintf('\t%s\n','Error...');
    % who = set of available variables as a cell array
    dopOSCCIcrash(err,gl,mfilename);
end

clear gr inputs