function varargout = bst_wavedata_display(D,OPTIONS);
%BST_WAVEDATA_DISPLAY - Various ways of loooking at time series
% function varargout = bst_wavedata_display(D,OPTIONS);
% Inputs:
%
% Outputs:

%<autobegin> ---------------------- 09-Jul-2004 22:16:52 -----------------------
% --------- Automatically Generated Comments Block Using AUTO_COMMENTS ---------
%
% CATEGORY: Visualization
%
% Alphabetical list of external functions (non-Matlab):
%   toolbox\bestfitsph.m
%   toolbox\bst_message_window.m
%   toolbox\bst_static_taskbar.m
%   toolbox\carto_sph.m
%   toolbox\findclosest.m
%   toolbox\get_user_directory.m
%   toolbox\interp_mail.m
%   toolbox\makeuswait.m
%   toolbox\norlig.m
%
% At Check-in: $Author: Mosher $  $Revision: 51 $  $Date: 7/09/04 8:42p $
%
% This software is part of BrainStorm Toolbox Version 2.0 (Alpha) 09-Jul-2004
% 
% Principal Investigators and Developers:
% ** Richard M. Leahy, PhD, Signal & Image Processing Institute,
%    University of Southern California, Los Angeles, CA
% ** John C. Mosher, PhD, Biophysics Group,
%    Los Alamos National Laboratory, Los Alamos, NM
% ** Sylvain Baillet, PhD, Cognitive Neuroscience & Brain Imaging Laboratory,
%    CNRS, Hopital de la Salpetriere, Paris, France
% 
% See BrainStorm website at http://neuroimage.usc.edu for further information.
% 
% Copyright (c) 2004 BrainStorm by the University of Southern California
% This software distributed  under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPL
% license can be found at http://www.gnu.org/copyleft/gpl.html .
% 
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%<autoend> ------------------------ 09-Jul-2004 22:16:52 -----------------------



% /---Script Author--------------------------------------\
% |                                                      |
% | *** Sylvain Baillet, Ph.D.                           |
% | Cognitive Neuroscience & Brain Imaging Laboratory    |
% | CNRS UPR640 - LENA                                   | 
% | Hopital de la Salpetriere, Paris, France             |
% | sylvain.baillet@chups.jussieu.fr                     |
% |                                                      |
% | John C. Mosher, Ph.D.                                |
% \------------------------------------------------------/


% Date of creation: July 2003
%
% Script History----------------------------------------------------
% SB     Jul-2003  Creation
% JCM 08-Sep-2003  Debugging, commenting
% JCM 11-Dec-2003  change "prod" to "all" for "ishandle" compatability in
%                  Matlab 6.5
% SB  15-Dec-2003  Fixed issues when plotting EEG only
% ------------------------------------------------------------------

% Default display options
Def_OPTIONS = struct(...
	'AbsoluteColormap', 0,...
	'Axes',[],...
	'Channel', [],...
	'Colorbar', 1,...
    'Colormap', [],...
    'DisplayTypes', [1 0 0], ...
    'DataType','surface',...
	'Fmax',[],...
	'hSurface', NaN*ones(1,length(D)), ...
	'hSensors', NaN, ...
	'WaveLabel',' ',...
	'IndepTopo',[],...
	'LineColor',[],...    
	'ModalityLabel',' ',...
	'Modality', 1,...
    'MoveCursor','on',...
	'nColumnLayout',1,...
	'SelectedChannels', {[]},...
	'SensorMarkers', 1,...
	'SensorLabels', 0,...
	'ShowContours', 1, ...
	'CurrentTime',[],...
	'Time', []...
	);

fontcolor = [1 1 1]; % Default font color for axes ticks and labels.

if nargin == 0 % Return default disply options
	if nargout > 1
		PlotHandles = Def_OPTIONS;
	else
		PlotHandles = [];
	end
	return
elseif nargin < 1 | nargin > 2
	errordlg(sprintf('Wrong number of arguments when calling %s', mfilename))
	return
elseif nargin == 1 % Use default option values
	
	OPTIONS = Def_OPTIONS;
	
elseif nargin == 2
	
	% Check field names of passed OPTIONS and fill missing ones with default values
	DefFieldNames = fieldnames(Def_OPTIONS);
	for k = 1:length(DefFieldNames)
		if ~isfield(OPTIONS,DefFieldNames{k}) | strcmp(DefFieldNames{k},'BEM')
			if ~isfield(OPTIONS,DefFieldNames{k})
				
				OPTIONS = setfield(OPTIONS,DefFieldNames{k},getfield(Def_OPTIONS,DefFieldNames{k}));
				
			end
		end
	end
end

switch(OPTIONS.DataType)
case 'surface'
	OPTIONS.DataType = 1;
case 'source'
	OPTIONS.DataType = 2;
end

if isempty(OPTIONS.Time)
	OPTIONS.Time = 1:size(D{1},2);
end

if isempty(OPTIONS.CurrentTime)
	OPTIONS.CurrentTime = mean(OPTIONS.Time);
end


nSets = length(D); % Number of data sets
nDispTypes = length(find(OPTIONS.DisplayTypes~=0)); % Number of display types

OverlayAllChannels = 0;
ColumnLayout = 0;

PlotHandles.TimeCursor = [];
PlotHandles.TimeZeroLine = [];

if nargin < 2
	PlotHandles.LineColor = [];
	PlotHandles.Label = cellstr([1:length(OPTIONS.Modality)]'); 
else
	PlotHandles.LineColor = OPTIONS.LineColor;
	PlotHandles.Label = OPTIONS.ModalityLabel;
end


% Display preparation
if isempty(OPTIONS.Fmax)
	for iSet = 1:nSets
		for mod = OPTIONS.Modality %1:3
			if ~isempty(OPTIONS.SelectedChannels{mod})
				if OPTIONS.DataType == 1
					OPTIONS.Fmax(iSet,mod) = max(max(abs(D{iSet}(OPTIONS.SelectedChannels{mod},:))));
				elseif OPTIONS.DataType == 2
					OPTIONS.Fmax(iSet,mod) = max(max(abs(D{iSet}(:))));
				end            
			else
				OPTIONS.Fmax(iSet,mod) = 0;
			end
		end
	end
end

OPTIONS.FmaxOrig = OPTIONS.Fmax;
DOrig = D;

fontsize = 8;%min([28/(nColumnLayout*(log(power(length(OPTIONS.Modality),.1)+1))),7]);

% Creation of graphics when they do not exist yet
% time series OVERLAY %--------------------------------------------------------------

makeuswait('start')

if OPTIONS.DisplayTypes(1) == 1
	
	OverlayAllChannels = 1;
	
	for iSet = 1:nSets % For each dataset
		
		F = D{iSet};
		
		iplot = 0;
		imod = 0;
		
		for mod = OPTIONS.Modality % For every selected data type
			imod = imod+1;
			iplot= iplot+1;
			
            if isempty(OPTIONS.SelectedChannels)
                OPTIONS.SelectedChannels{mod} = [1:size(F,1)]; % Assume F contains waveforms of a single type
            end
            
            % Data normalization
			if iSet == 1
				%                 tmp = [OPTIONS.Fmax{:}];
				%                 tmp = reshape(tmp,3,nSets)';
				%                 maxData(mod) = max(tmp(:,mod));
				if OPTIONS.DataType == 1
                    maxData(mod) = max(OPTIONS.Fmax(:,mod));
                else
                    maxData(mod) = max(OPTIONS.Fmax(:));
                end
                
				% Sanity check
				if isnan(maxData(mod)) | maxData == 0
					maxData(mod) = eps;
					OPTIONS.Fmax(iSet,mod) = eps;
					if OPTIONS.DisplayTypes == 1
						DataToDisplay = zeros(size(F(OPTIONS.SelectedChannels{mod},:)));
					elseif OPTIONS.DisplayTypes == 2
						DataToDisplay = zeros(size(F(:)));
					end
				end
            
                    
            end
			
			if OPTIONS.DataType == 1
				DataToDisplay = F(OPTIONS.SelectedChannels{mod},:);
			elseif OPTIONS.DataType == 2 
				DataToDisplay = F;
			end
			
			DataToDisplay = DataToDisplay/maxData(mod);
			
			% Display parameters
			N = length(OPTIONS.Modality); % Maximum number of channels to visualize per column
			delta = .95/(N+1); % Space between plots (in normalized units)  
			
			if iplot == 1 & iSet == 1
				M = 1;%max(max(F(channels,:))); % Useful as a scaling factor for display purposes
				PlotHandles.PlotOverlay(mod).axes = subplot(1,nDispTypes,1);
				hold on
			else
				PlotHandles.PlotOverlay(mod).axes = PlotHandles.PlotOverlay(OPTIONS.Modality(1)).axes;%subplot(1,ncol,col+PlotHandles.OverlayAllChannels);
			end
			
			PlotHandles.PlotOverlay(mod).lines{iSet} = plot(OPTIONS.Time,DataToDisplay'+(2*M*iplot));
		                  
            
			if length(OPTIONS.Modality) > 1
                Line(mod) = PlotHandles.PlotOverlay(mod).lines{iSet}(1); % Keep track of line handle for each modality
            else
                Line(iSet) = PlotHandles.PlotOverlay(mod).lines{iSet}(1); % Keep track of line handle for each set for proper legend (see below)    
            end
            
            if nSets > 1 | length(OPTIONS.Modality) > 1 % Assign a single color to all waveforms of each data set / modality
				
				if iSet == 1 & isempty(OPTIONS.LineColor)
					% Generate a set of colors for time series
					PlotHandles.LineColor = get(0,'DefaultAxesColorOrder'); % Get default line colors
               % JCM: Gotta fix this cheat
               if 0
                  PlotHandles.LineColor = PlotHandles.LineColor(1:nSets,:); % CHEAT , nSets need to be smaller than size(get(0,'DefaultAxesColorOrder'),1)
               else
                  % JCM fix 16-June-2004
                  temp_nSets = rem(0:(nSets-1),size(PlotHandles.LineColor,1))+1; % modulo scheme for line colors
                  PlotHandles.LineColor = PlotHandles.LineColor(temp_nSets,:); 
               end
               %PlotHandles.LineColor(PlotHandles.LineColor>.8) = .8; % Avoid too bright colors
               PlotHandles.LineColor = num2cell(PlotHandles.LineColor,2);
            elseif length(OPTIONS.LineColor)~= nSets
               PlotHandles.LineColor = get(0,'DefaultAxesColorOrder');
               % JCM: Gotta fix this cheat
               if 0
                  PlotHandles.LineColor = PlotHandles.LineColor(1:nSets,:);
               else
                  % JCM fix 16-June-2004
                  temp_nSets = rem(0:(nSets-1),size(PlotHandles.LineColor,1))+1; % modulo scheme for line colors
                  PlotHandles.LineColor = PlotHandles.LineColor(temp_nSets,:); 
               end
               
               %PlotHandles.LineColor = rand(nSets,3); 
               %PlotHandles.LineColor(PlotHandles.LineColor>.8) = .8; % Avoid too bright colors
               PlotHandles.LineColor = num2cell(PlotHandles.LineColor,2);
				end
				
			else
				
				PlotHandles.LineColor{1} = .2*ones(1,3); % Default time series color
				
			end
						
			set(PlotHandles.PlotOverlay(mod).lines{iSet},'Color',PlotHandles.LineColor{iSet},'ButtonDownFcn','dataplot_cb NameSelectedChannel',...
				'Visible','on','Userdata',iSet)
			
			if iSet == 1 % Display data scale for each modality
				ratio = .2;
				% Take ratio X maxData amplitude as a marker of data amplitude
				itime = OPTIONS.Time(ceil(length(OPTIONS.Time)*ratio));
                if OPTIONS.DataType == 1
                    PlotHandles.PlotOverlay(mod).YScale.Line = errorbar(itime,2.5*M*iplot+ratio,ratio/2,'k');
                    if mod == 1 % MEG
						sstring = 'fT'; factor = 1e-15;
					elseif mod == 2 % EEG
						sstring  = 'mV'; factor = 1;%1e-3;
                    else 
						sstring = 'a.u.'; factor = 1;
					end
                    
                    PlotHandles.PlotOverlay(mod).YScale.Label = text(itime,2.5*M*iplot+1.5*ratio,sprintf('%4.1f %s',(ratio*maxData(mod)/factor),sstring));
                    set(PlotHandles.PlotOverlay(mod).YScale.Label,'HorizontalAlignment','left','fontsize', fontsize,'VerticalAlignment','bottom')
                    
                elseif OPTIONS.DataType == 2 & mod == OPTIONS.Modality(1)
                    
                    PlotHandles.PlotOverlay(mod).YScale.Line = errorbar(itime,2.5*M*iplot+ratio,ratio/2,'k');
                    sstring = 'A.m'; factor = 1; % CBB | manage ZScores
                    
                    %                     elseif mod == 2 % Source time series - ZScored
                    %                         sstring = texlabel('sigma'); factor = 1;
                    %end
                    
                    PlotHandles.PlotOverlay(mod).YScale.Label = text(itime,2.5*M*iplot+1.5*ratio,sprintf('%4.1e %s',(ratio*maxData(mod)/factor),sstring));    
                    
                    set(PlotHandles.PlotOverlay(mod).YScale.Label,'HorizontalAlignment','left','fontsize', fontsize,'VerticalAlignment','bottom')
                    
                end
                
                
                YZeroLine = line(get(PlotHandles.PlotOverlay(mod).axes,'Xlim'),2*M*iplot*[1 1],'color',.8*[1 1 1]);
                
            end
        end 
        
    end % For each selected dataset
    
    
    set(PlotHandles.PlotOverlay(mod).axes,'units','normalized','ygrid','on','xgrid','on','Box','on',...
        'Yticklabel',PlotHandles.Label,'Ytick',2*M*[1:iplot],...
		'FontName','default','FontSize',fontsize,'FontUnits','Normal','Fontweight','light',...
		'xColor',fontcolor, 'ycolor', fontcolor,'XMinorGrid','on',...
		'Xlim',[OPTIONS.Time(1),OPTIONS.Time(end)],'Ylim',[0 2*M*(N+1)],...
		'Visible','on')
	
	a =  xlabel('Time (s)');
	set(a,'Visible','on','fontsize',fontsize,'Fontunits','normal') % Apparently has no effect in Matlab 6.3 
	
	axes(PlotHandles.PlotOverlay(mod).axes)
	[LEGH,OBJH,OUTH,OUTM] = legend(Line,OPTIONS.WaveLabel);
	legend boxoff
	
	LegendMenu = uicontextmenu;
	uimenu(LegendMenu,'Label','Hide/Show Legend','Callback','dataplot_cb HideShowMouseSelectedDisplayLegend')
	uimenu(LegendMenu,'checked','off','Label','Zoom ON','Callback','dataplot_cb ToggleWindowFromZoom')
	uimenu(LegendMenu,'checked','off','Label','Rotate','Callback','dataplot_cb ToggleRotate3D')
	uimenu(LegendMenu,'Separator','on','Label','Copy Current Plot','Callback','dataplot_cb CopyCurrentPlot')
	set([gcf;LEGH;OBJH],'UIcontextMenu',LegendMenu)
	set(findobj(OBJH,'type','text'),'Fontsize',7,'fontunits','points')       
	
	PlotHandles.TimeZeroLine(1) = line([0 0],get(gca,'Ylim'));
	set(PlotHandles.TimeZeroLine(1),'Visible','on')
	PlotHandles.TimeCursor(1) = line([OPTIONS.CurrentTime OPTIONS.CurrentTime],get(gca,'Ylim'),'linewidth',1);
	
end


% Column layout of time series % ---------------------------------------------------------------------------
if OPTIONS.DisplayTypes(2) > 0
	
	ColumnLayout = 1; % Flag
	OPTIONS.nColumnLayout = OPTIONS.DisplayTypes(2); % Number of columns
	
	% Initialize column plot structure
	tmp = deal(struct('lines',{cell(1,nSets)},'axes',zeros(OPTIONS.nColumnLayout,1),'ytick',[]));
	for mod = OPTIONS.Modality % For every selected data type
		PlotHandles.PlotColumn(mod) = tmp;
     end
	clear tmp
	
    
	if OPTIONS.DataType == 2 % Source column-display 

        nSets = length(D);
		        
        % Display preparation
		OPTIONS.Fmax = [];
		for iSet = 1:nSets
			imod = 0;
            for mod = OPTIONS.Modality
                imod = imod + 1;
        
                if isempty(OPTIONS.SelectedChannels)
                    OPTIONS.SelectedChannels{mod} = [1:size(D{iSet},1)]; % Assume F contains waveforms of a single type
                end

                if ~isempty(OPTIONS.SelectedChannels{mod})
					OPTIONS.Fmax(iSet,imod) = max(max(abs(D{iSet}(:))));
				else
					OPTIONS.Fmax(iSet,imod) = 0;
				end
			end
		end
        
        D = {cat(1,D{:})};
        nSets = length(D);
	       
    end
	
	for iSet = 1:nSets % For each data set (MEG/EEG or Original/Modeled/Residuals) or each source
		
		F = D{iSet};
		
		% Data normalization
		channels = [];
		imod = 0;
		for mod = OPTIONS.Modality % For every selected data type
			imod = imod+1;
			if iSet == 1
				%                 tmp = [OPTIONS.Fmax{:}];
				%                 tmp = reshape(tmp,3,nSets)';
				%                 maxData(mod) = max(tmp(:,mod));
				maxData(mod) = max(OPTIONS.Fmax(:,mod));
				
				% Sanity check
				if isnan(maxData(mod))
					maxData(mod) = eps;
					OPTIONS.Fmax(iSet,mod) = eps;
					if OPTIONS.DisplayTypes == 1
						DataToDisplay = zeros(size(F(OPTIONS.SelectedChannels{mod},:)));
					elseif OPTIONS.DisplayTypes == 2
						DataToDisplay = zeros(size(F(:)));
					end
				end
			end
			
			if OPTIONS.DataType == 1
				F(OPTIONS.SelectedChannels{mod},:) = D{iSet}(OPTIONS.SelectedChannels{mod},:)/maxData(mod);
				channels = [channels,OPTIONS.SelectedChannels{mod}];  
			elseif OPTIONS.DataType == 2 
				DataToDisplay = D{iSet}/maxData(mod);
				channels = [1:size(F,1)];
			end
		end
		
		
		if iSet == 1
			M = 1;%max(max(abs(F(channels,:)))); % Useful as a scaling factor for display purposes
			N = ceil(length(channels)/OPTIONS.nColumnLayout); % Maximum number of channels to visualize per column
			fontsize = min([28/(OPTIONS.nColumnLayout*(log(power(N,.1)+1))),7]);
			delta = .95/(N+1); % Space between plots (in normalized units)  
		end
		
		for col = 1:OPTIONS.nColumnLayout % For each column of the layout
			
			if col < OPTIONS.nColumnLayout
				try
					IDscol = channels((col-1)*N+1 : N*col); % Labels of channels to plot in each column 
				catch
					errordlg('Please decrease the number of columns.')
					return
				end
				
			else
				IDscol = channels((col-1)*N+1 : end); % Labels of channels to plot in each column 
			end
			
			i = 0;
			
			for k = IDscol % For each channel in current column
				
                if OPTIONS.DataType == 1 % For sensor data - check modality type for each waveform
                  
                    if ~isempty(find(OPTIONS.SelectedChannels{1}==k))
                        mod = 1;
                    elseif ~isempty(find(OPTIONS.SelectedChannels{2}==k)) 
                        mod = 2;
                    elseif ~isempty(find(OPTIONS.SelectedChannels{3}==k)) 
                        mod = 3;
                    end
                    
                else % For source time series - force modality
                
                    mod = OPTIONS.Modality;
                    
                end
                    
				
				if k == IDscol(1) & iSet == 1; % First channel of current subset : create new axes
					PlotHandles.PlotColumn(mod).axes(col) = subplot(1,nDispTypes*OPTIONS.nColumnLayout,OPTIONS.nColumnLayout*OverlayAllChannels + col);%axes;
					hold on
				elseif iSet == 1
					PlotHandles.PlotColumn(mod).axes(col) = gca;
				elseif iSet >1
					axes(PlotHandles.PlotColumn(mod).axes(col));% = gca;
				end
				
				i = i+1; % Index of channel plot in current axes
				
				if iSet > 1
					axes( PlotHandles.PlotColumn(mod).axes(col))
				end
				
				if OPTIONS.DataType == 1
					tmp = plot(OPTIONS.Time,F(k,:)+(2*M*i));
				elseif OPTIONS.DataType == 2
					tmp = plot(OPTIONS.Time,DataToDisplay(k,:)+(2*M*i));
                    
                    YZeroLine = line(get(PlotHandles.PlotColumn(mod).axes(col),'Xlim'),...
                        2*M*i * [1,1],'color',[.8 .8 .8]);
          
                end
		
                
                
				Line(iSet) = tmp(1); % Keep track of a sample line handle from each set (see legend display below)                    
				if nSets > 1
					if col == 1 & k == IDscol(1)
						if isempty(PlotHandles.LineColor) % Line color for current data set not defined 
							PlotHandles.LineColor = rand(nSets,3); 
							PlotHandles.LineColor(PlotHandles.LineColor>.8) = .8; % Avoid too bright colors
							PlotHandles.LineColor = num2cell(PlotHandles.LineColor,2);
						end
					end
				else
                    if isempty(PlotHandles.LineColor) % Line color for current data set not defined 
                        PlotHandles.LineColor{iSet} = .2*ones(1,3); % Default color
                    end
                end
                
				set(tmp,'color',PlotHandles.LineColor{iSet},'ButtonDownFcn','dataplot_cb NameSelectedChannel') % Display channel Name when time series is selectedend
				
				PlotHandles.PlotColumn(mod).lines{iSet} = [PlotHandles.PlotColumn(mod).lines{iSet}; tmp];
				
				PlotHandles.PlotColumn(mod).ytick = [PlotHandles.PlotColumn(mod).ytick,2*M*i];
				
				hold on
				
			end % for each channel
			
			if iSet == 1
				if OPTIONS.DataType == 1
					YTickLabel = {OPTIONS.Channel(IDscol).Name};
				elseif OPTIONS.DataType == 2
					YTickLabel = OPTIONS.WaveLabel(IDscol);
				end
				
				set(gca,'Visible','on',...
					'Yticklabel',YTickLabel,'Ytick',2*M*[1:length(IDscol)],...
					'units','normalized','ygrid','on','xgrid','on','Box','on',...
					'FontName','default','FontSize',fontsize,'FontUnits','Normal','Fontweight','normal',...
					'xColor',fontcolor, 'ycolor', fontcolor,'XMinorGrid','on',...
					'Xlim',[OPTIONS.Time(1),OPTIONS.Time(end)],'Ylim',[0 2*M*(N+1)])
				
				%set(gca,'fontunits','points') % Move back to normal because channel names are not readable enough when window size is increased
				
				PlotHandles.TimeCursor(end+1) = line([OPTIONS.CurrentTime OPTIONS.CurrentTime],get(gca,'Ylim'),'linewidth',1);
				set(PlotHandles.TimeCursor(end),'EraseMode','normal','Visible','on')
				PlotHandles.TimeZeroLine(end+1) = line([0 0],get(gca,'Ylim'),'linewidth',1);
                
            
                
                if OPTIONS.DataType == 2
                    if mod == 1 % Source time series - A.m.
                        sstring = 'A.m'; factor = 1;
                    elseif mod == 2 % Source time series - ZScored
                        sstring = texlabel('sigma'); factor = 1;
                    end

                    ratio = .2;
                    % Take ratio X maxData amplitude as a marker of data amplitude
                    iplot = 1;
                    itime = OPTIONS.Time(ceil(length(OPTIONS.Time)*ratio));
                    PlotHandles.PlotColumn(mod).YScale.Line = errorbar(itime,2.5*M*iplot+ratio,ratio/2,'k');

                    PlotHandles.PlotColumn(mod).YScale.Label = text(itime,2.5*M*iplot+1.5*ratio,sprintf('%4.1e %s',(ratio*maxData(mod)/factor),sstring));
                    set(PlotHandles.PlotColumn(mod).YScale.Label,'HorizontalAlignment','left','fontsize', fontsize,'VerticalAlignment','bottom')
                end
                
                
            end
			
			xlabel('Time (s)')
			
		end % for each column
		
		
		set(PlotHandles.PlotColumn(mod).lines{iSet},'UserData',iSet)
		
		
	end % for each data set
	

    % Display Legend
    axes(PlotHandles.PlotColumn(mod).axes(col))
    if OPTIONS.DataType == 1

		[LEGH,OBJH,OUTH,OUTM] = legend(Line,OPTIONS.WaveLabel);
        
    elseif OPTIONS.DataType == 2

		[LEGH,OBJH,OUTH,OUTM] = legend(PlotHandles.PlotColumn(mod).lines{1},OPTIONS.WaveLabel);
	      
    end
	
    
    try 
        legend boxoff
    end
    LegendMenu = uicontextmenu;
    uimenu(LegendMenu,'Label','Hide/Show Legend','Callback','dataplot_cb HideShowMouseSelectedDisplayLegend')
    uimenu(LegendMenu,'checked','off','Label','Zoom ON','Callback','dataplot_cb ToggleWindowFromZoom')
    uimenu(LegendMenu,'Separator','on','Label','Copy Current Plot','Callback','dataplot_cb CopyCurrentPlot')
    set([gcf;LEGH;OBJH],'UIcontextMenu',LegendMenu)
    
    set(findobj(OBJH,'type','text'),'Fontsize',7,'fontunits','points')       

    
end


% Surface Topography 
if OPTIONS.DisplayTypes(3) > 1 % Spatial data display requested
	
	if OPTIONS.DataType == 2 % Results
		D = DOrig;
		nSets = length(D);
		OPTIONS.Fmax =  OPTIONS.FmaxOrig;
	end
	
	for iSet = 1:nSets % For each dataset /source
		
		if OPTIONS.DataType == 1
			F = D{iSet};
		elseif OPTIONS.DataType == 2;
			F = zeros(length(OPTIONS.Channel),size(D{iSet},2));
			F(OPTIONS.SelectedChannels{mod},:) = OPTIONS.IndepTopo(:,iSet)*D{iSet};
			if iSet == 1 &~OPTIONS.AbsoluteColormap 
				% Compute absolute maximum of surface topography for all sources 
				tmp = [D{:}]'; tmp = reshape(tmp,length(D{1}),length(D));
				Topo_maxx = max(max(abs(OPTIONS.IndepTopo*tmp'))); clear tmp
			end
		end
		
		% Surface data mapping (topography)
		iplot = 0;
		imod = 0;
		for mod = OPTIONS.Modality % For every selected data type
			imod = imod+1;
			iplot= iplot+1;
			% Data normalization
			
			if OPTIONS.DataType == 1
				chanlocs = [OPTIONS.Channel(OPTIONS.SelectedChannels{mod}).Loc]';
			elseif OPTIONS.DataType == 2
				chanlocs = [OPTIONS.Channel.Loc]';
			end
			
            if isempty(chanlocs)
				return
			end
			
            
			if OPTIONS.DataType == 1
				if size([OPTIONS.Channel(OPTIONS.SelectedChannels{mod}(1)).Loc],2) == 2 % Gradiometers
					chanlocs = chanlocs(1:2:end,:);
				end
			elseif OPTIONS.DataType == 2
				if size(OPTIONS.Channel(1).Loc,2) == 2 % Gradiometers
					chanlocs = chanlocs(1:2:end,:);
				end
			end
			
            if size(chanlocs,1) ~= length(OPTIONS.SelectedChannels{mod})
                errordlg('Channel locations do not match the number of channels in Channel structure.','Perverted Channel structure.')
                return
            end
            
          
            
			maxData(mod) = max(max(abs(chanlocs(:,2))));
			chanlocs_orig = chanlocs;
			chanlocs = chanlocs/maxData(mod);
			
			% Display parameters
			N = length(OPTIONS.Modality); % Maximum number of channels to visualize per column
			delta = .95/(N+1); % Space between plots (in normalized units)  
			
			if iplot == 1
				M = (max(abs([chanlocs]))); % Useful as a scaling factor for display purposes
				M = M(2);
				if isempty(OPTIONS.Axes)
					PlotHandles.TopoPlot(mod).axes{iSet} = ...
						subplot(length(OPTIONS.Modality),nDispTypes*nSets,nSets*(ColumnLayout+ OverlayAllChannels) + iSet);
				else
					PlotHandles.TopoPlot(mod).axes{iSet} = OPTIONS.Axes;
				end
				hold on
			else
				if isempty(OPTIONS.Axes)
					PlotHandles.TopoPlot(mod).axes{iSet} = ...
					subplot(length(OPTIONS.Modality),nDispTypes*nSets,nDispTypes*nSets + nSets*(ColumnLayout + OverlayAllChannels) + iSet);
				else
					PlotHandles.TopoPlot(mod).axes{iSet} = OPTIONS.Axes;
				end
				
				hold on
				pos1 = get(PlotHandles.TopoPlot(mod-1).axes{iSet},'Position');
				pos2 = get(PlotHandles.TopoPlot(mod).axes{iSet},'Position');
			end
			PlotHandles.TopoPlot(mod).Colorbar{iSet} = [];
			set(PlotHandles.TopoPlot(mod).axes{iSet},'Visible','off')
			
			tindx = findclosest(OPTIONS.CurrentTime,OPTIONS.Time');
			
			switch(OPTIONS.DisplayTypes(3)) 
			case 1 % Do nothing
			case 2 % Spherical approximation of the scalp
				if ~ishandle(OPTIONS.hSurface(iSet))
					if iSet == 1
						[PlotHandles.PChanLocs{mod},hsens,hsph,CData,Wmat] = carto_sph(chanlocs,F(OPTIONS.SelectedChannels{mod},tindx));
						PlotHandles.TopoPlot(mod).surf{iSet} = hsph;
						PlotHandles.TopoPlot(mod).Wmat{iSet} = Wmat;
						%PlotHandles.TopoPlot(mod).hsens{iSet} = hsens;
						clear hpsh Wmat hsens;
					else % Just copy previous surface object (much faster display !)
						PlotHandles.TopoPlot(mod).surf{iSet} = copyobj(PlotHandles.TopoPlot(mod).surf{1},PlotHandles.TopoPlot(mod).axes{iSet});
						PlotHandles.TopoPlot(mod).Wmat{iSet} = PlotHandles.TopoPlot(mod).Wmat{1};
						%PlotHandles.TopoPlot(mod).hsens{iSet} = PlotHandles.TopoPlot(mod).hsens{1};
						CData = PlotHandles.TopoPlot(mod).Wmat{iSet}*F(OPTIONS.SelectedChannels{mod},tindx);
						CData = reshape(CData,sqrt(length(CData)),sqrt(length(CData)));
						CData = CData * max(abs(F(OPTIONS.SelectedChannels{mod},tindx)))/max(abs(CData(:))); % pr
						set( PlotHandles.TopoPlot(mod).surf{iSet} ,'CData',CData); clear CData
					end
				else
					CData = OPTIONS.TopoPlot(mod).Wmat{iSet}*F(OPTIONS.SelectedChannels{mod},tindx);
					set(OPTIONS.hSurface(iSet),'CData',CData); clear CData
				end
				
				if OPTIONS.SensorMarkers & ~ishandle(OPTIONS.hSensors)
					if iSet == 1
						for kk = 1:size(PlotHandles.PChanLocs{mod}(:,1),1) % SB: I don't use scatter here because ButtonDownFcn is not effective as for Matlab 6.1
							PlotHandles.TopoPlot(mod).hsens{iSet}(kk) = plot3(1.05*PlotHandles.PChanLocs{mod}(kk,1),1.05*PlotHandles.PChanLocs{mod}(kk,2),1.05*PlotHandles.PChanLocs{mod}(kk,3),...
								'.','Markersize',2);
						end
					else % Just copy previous sensor markers (faster display)
						PlotHandles.TopoPlot(mod).hsens{iSet} = copyobj(PlotHandles.TopoPlot(mod).hsens{1},PlotHandles.TopoPlot(mod).axes{iSet});
					end
					
				end
				if OPTIONS.SensorLabels
					if iSet == 1
						for  kk = 1:size(PlotHandles.PChanLocs{mod}(:,1),1) % SB: I don't use scatter here because ButtonDownFcn is not effective as for Matlab 6.1
							PlotHandles.TopoPlot(mod).SensorLabels{iSet}(kk) = text(1.05*PlotHandles.PChanLocs{mod}(kk,1),1.05*PlotHandles.PChanLocs{mod}(kk,2),1.05*PlotHandles.PChanLocs{mod}(kk,3),...
								OPTIONS.Channel(OPTIONS.SelectedChannels{mod}(kk)).Name);
						end
					else
						PlotHandles.TopoPlot(mod).SensorLabels{iSet} = copyobj(PlotHandles.TopoPlot(mod).SensorLabels{1},PlotHandles.TopoPlot(mod).axes{iSet});
					end
					
				end
				
			case 3 % Sensor Cap
				
				if iSet == 1
					chanlocs(:,3)=chanlocs(:,3)-max(chanlocs(:,3));
					[TH,PHI,R]=cart2sph(chanlocs(:,1),chanlocs(:,2),chanlocs(:,3));
					%PHI2=zeros(size(PHI));
					R2=R./cos(PHI).^.2;
					[Y,X]=pol2cart(TH,R2);
					%Y et X sont les coordonnées projetees dans le plan 'tangent superieur' pour chaque capteur
					bord = convhull(Y,X);
					ncapt = size(chanlocs,1);
					[center,R] = bestfitsph(chanlocs);
					coordC = chanlocs-(ones(ncapt,1)*center');
					tri = convhulln(coordC./(norlig(coordC)'*ones(1,3)));
					keep = find(~(ismember(tri(:,1),bord) & ismember(tri(:,2),bord)& ismember(tri(:,3),bord)));
					tri = tri(keep,:);
					
					PlotHandles.TopoPlot(mod).surf{iSet} = trisurf(tri,chanlocs(:,1),chanlocs(:,2),chanlocs(:,3),...
						F(OPTIONS.SelectedChannels{mod},tindx));
					
					hold on
					set(PlotHandles.TopoPlot(mod).surf{iSet},'edgecolor','none','facecolor','interp')
				else
					PlotHandles.TopoPlot(mod).surf{iSet} = copyobj(PlotHandles.TopoPlot(mod).surf{1},PlotHandles.TopoPlot(mod).axes{iSet});
					set( PlotHandles.TopoPlot(mod).surf{iSet} ,'FaceVertexCData',F(OPTIONS.SelectedChannels{mod},tindx)); 
				end
				
				
				if OPTIONS.SensorMarkers
					if iSet == 1
						for kk = 1:size(chanlocs,1) % SB: I don't use scatter here because ButtonDownFcn is not effective as for Matlab 6.1
							PlotHandles.TopoPlot(mod).hsens{iSet}(kk) = plot3(chanlocs(kk,1),chanlocs(kk,2),chanlocs(kk,3),'.');
						end
					else
						PlotHandles.TopoPlot(mod).hsens{iSet} = copyobj(PlotHandles.TopoPlot(mod).hsens{1},PlotHandles.TopoPlot(mod).axes{iSet});
					end
					
				end
				if OPTIONS.SensorLabels
					if iSet == 1
						for  kk = 1:size(PlotHandles.PChanLocs{imod}(:,1),1) % SB: I don't use scatter here because ButtonDownFcn is not effective as for Matlab 6.1
							PlotHandles.TopoPlot(mod).SensorLabels{iSet}(kk) = text(chanlocs(kk,1),chanlocs(kk,2),chanlocs(kk,3),...
								OPTIONS.Channel(OPTIONS.SelectedChannels{mod}(kk)).Name);
						end
					else
						PlotHandles.TopoPlot(mod).SensorLabels{iSet} = copyobj(PlotHandles.TopoPlot(mod).SensorLabels{1},PlotHandles.TopoPlot(mod).axes{iSet});    
					end
					
				end
				
			case 4 % Scalp Surface
				
				% Load scalp envelope information
				[TessFile,h] = bst_static_taskbar('GET','TESS'); clear h
				if isempty(TessFile)
                    errordlg('No tessellated envelope is available for this subject. Please edit this subject''s entry in database')
                    return
                end
                
                bst_message_window('wrap',{...
						' ',...
						sprintf('Loading %s tessellation file',TessFile)});
				Users = get_user_directory;
				load(fullfile(Users.SUBJECTS,TessFile),'Comment')
				
				PlotHandles.TopoPlot(mod).iGrid = 1;
                
                if length(Comment)>1 % Several surfaces in this tessellation file
					[PlotHandles.TopoPlot(mod).iGrid,OK] = listdlg('PromptString','Select an envelope from list:',...
						'SelectionMode','single',...
						'ListString',Comment);
					if OK == 0
						return
					end
				end
				  
                PlotHandles.TopoPlot(mod).SurfaceName{iSet} = Comment{PlotHandles.TopoPlot(mod).iGrid};
                clear Comment
              
                load(fullfile(Users.SUBJECTS,TessFile),'Faces','Vertices')
                Faces = Faces{PlotHandles.TopoPlot(mod).iGrid};
                Vertices = Vertices{PlotHandles.TopoPlot(mod).iGrid}';
                % Warn user if tessellation is huge
                if size(Vertices,1) > 20000
                    ButtonName = questdlg({sprintf('Tessellation has %d vertices; computation time for rendering might be excessive.',size(Vertices,1)),'Do you still want to proceed ?'},'Long computation expected','Yes','No','No');
                    switch(ButtonName)
                    case 'No'
                        return
                    end
                end
                
                % Load scalp envelope information <- DONE
                
                PlotHandles.PChanLocs{mod} = chanlocs_orig;
                % 				% Project sensor locations on scalp envelope
                % 				j = 0; minn = []; imin = [];
                % 				for chan = 1:size(chanlocs_orig,1) 
                % 					celec = ones(size(Vertices,1),1)*chanlocs_orig(chan,:);
                % 					dist = norlig(Vertices-celec);
                % 					[minn(chan) imin(chan)] = min(dist);
                % 				end
                % 				PlotHandles.PChanLocs{mod} = [Vertices(imin,1),Vertices(imin,2),Vertices(imin,3)];
                
                
                axes(PlotHandles.TopoPlot(mod).axes{iSet})
				PlotHandles.TopoPlot(mod).surf{iSet} = patch('faces',Faces,'vertices',Vertices);  clear Faces              
                set(PlotHandles.TopoPlot(mod).surf{iSet},'visible','off')
				
                %vertxcolor = interp_mail(Vertices,PlotHandles.PChanLocs{mod},F(OPTIONS.SelectedChannels{mod},tindx));
                bst_message_window('wrap',{...
                        '-',...
                        'Interpolating data on scalp surface. . .',...
                        '-'})

                allfigs = findobj(0,'type','figure'); % Find all figures in workspace...
                     
                [vertxcolor,PlotHandles.TopoPlot(mod).Wmat{iSet}] = interp_mail(Vertices,PlotHandles.PChanLocs{mod},F(OPTIONS.SelectedChannels{mod},tindx));
                bst_message_window('wrap',...
                        '-> Done')

                
                set(PlotHandles.TopoPlot(mod).surf{iSet},'FaceVertexCData',vertxcolor,'facecolor','interp','edgecolor','none','backfacelighting','unlit','visible','on');
				material dull, lighting none
				%camlight 
				            
				if OPTIONS.SensorMarkers
					if iSet == 1
						for kk = 1:size(chanlocs,1) % SB: I don't use scatter here because ButtonDownFcn is not effective as for Matlab 6.1
							PlotHandles.TopoPlot(mod).hsens{iSet}(kk) = plot3(chanlocs_orig(kk,1),chanlocs_orig(kk,2),chanlocs_orig(kk,3),'.');
						end
					else
						PlotHandles.TopoPlot(mod).hsens{iSet} = copyobj(PlotHandles.TopoPlot(mod).hsens{1},PlotHandles.TopoPlot(mod).axes{iSet});
					end
					
				end
				
				if OPTIONS.SensorLabels
					if iSet == 1
						for  kk = 1:size(PlotHandles.PChanLocs{imod}(:,1),1) % SB: I don't use scatter here because ButtonDownFcn is not effective as for Matlab 6.1
							PlotHandles.TopoPlot(mod).SensorLabels{iSet}(kk) = text(chanlocs_orig(kk,1),chanlocs_orig(kk,2),chanlocs_orig(kk,3),...
								OPTIONS.Channel(OPTIONS.SelectedChannels{mod}(kk)).Name);
						end
					else
						PlotHandles.TopoPlot(mod).SensorLabels{iSet} = copyobj(PlotHandles.TopoPlot(mod).SensorLabels{1},PlotHandles.TopoPlot(mod).axes{iSet});    
					end
				end
				
				camlight 
				lighting gouraud
				
			case 5 % Disc
				
				if iSet == 1
					[PlotHandles.PChanLocs{mod},hsens,c,CData,e] = carto_sph(chanlocs,F(OPTIONS.SelectedChannels{mod},tindx));
					PlotHandles.TopoPlot(mod).surf{iSet} = c;
					PlotHandles.TopoPlot(mod).Wmat{iSet} = e; clear c e
					
					xs = get(PlotHandles.TopoPlot(mod).surf{iSet},'XData');
					ys = get(PlotHandles.TopoPlot(mod).surf{iSet},'YData');
					zs = get(PlotHandles.TopoPlot(mod).surf{iSet},'ZData');
					
					chanlocs = [xs(:),ys(:),zs(:)];
					chanlocs(:,3) =  chanlocs(:,3) - max(chanlocs(:,3));
					
					[TH,PHI,R] = cart2sph(chanlocs(:,1),chanlocs(:,2),chanlocs(:,3));
					R2 = R./cos(PHI).^.2;
					[Y,X] = pol2cart(TH,R2);
					
					tri = delaunay(Y,X);
					PlotHandles.TopoPlot(mod).surf{iSet} = patch('faces',tri,'vertices',[Y,X,0*X],'CData',get(PlotHandles.TopoPlot(mod).surf{iSet},'CData'),'edgecolor','none',...
						'FaceColor','interp');
					PlotHandles.TopoPlot(mod).Yo{iSet} = Y; 
					PlotHandles.TopoPlot(mod).Xo{iSet} = X;
					
				else
					
					PlotHandles.TopoPlot(mod).surf{iSet} = copyobj(PlotHandles.TopoPlot(mod).surf{1},PlotHandles.TopoPlot(mod).axes{iSet});
					PlotHandles.TopoPlot(mod).Wmat{iSet} = PlotHandles.TopoPlot(mod).Wmat{1};
					PlotHandles.TopoPlot(mod).Yo{iSet} = PlotHandles.TopoPlot(mod).Yo{1}(:); 
					PlotHandles.TopoPlot(mod).Xo{iSet} = PlotHandles.TopoPlot(mod).Xo{1}(:);
					
					CData = PlotHandles.TopoPlot(mod).Wmat{iSet}*F(OPTIONS.SelectedChannels{mod},tindx);
					CData = reshape(CData,sqrt(length(CData)),sqrt(length(CData)));
					CData = CData * max(abs(F(OPTIONS.SelectedChannels{mod},tindx)))/max(abs(CData(:))); 
					set( PlotHandles.TopoPlot(mod).surf{iSet},'CData',CData); 
					
				end
				
				
				PlotHandles.PChanLocs{mod}(:,3) =  PlotHandles.PChanLocs{mod}(:,3) - max(PlotHandles.PChanLocs{mod}(:,3));
				[TH,PHI,R] = cart2sph(PlotHandles.PChanLocs{mod}(:,1),PlotHandles.PChanLocs{mod}(:,2),PlotHandles.PChanLocs{mod}(:,3));
				R2 = R./cos(PHI).^.2;
				[Y,X] = pol2cart(TH,R2);
				
				if OPTIONS.SensorMarkers
					if iSet == 1
						for kk = 1:length(Y) % SB: I don't use scatter here because ButtonDownFcn is not effective as for Matlab 6.1
							PlotHandles.TopoPlot(mod).hsens{iSet}(kk) = plot3(Y(kk),X(kk),0.05,'.');
						end
					else
						PlotHandles.TopoPlot(mod).hsens{iSet} = copyobj(PlotHandles.TopoPlot(mod).hsens{1},PlotHandles.TopoPlot(mod).axes{iSet});                        
					end
				end
				if OPTIONS.SensorLabels
					if iSet == 1
						for  kk = 1:length(Y) % SB: I don't use scatter here because ButtonDownFcn is not effective as for Matlab 6.1
							PlotHandles.TopoPlot(mod).SensorLabels{iSet}(kk) = text(Y(kk),X(kk),0.05,...
								OPTIONS.Channel(OPTIONS.SelectedChannels{mod}(kk)).Name);
						end
					else
						PlotHandles.TopoPlot(mod).SensorLabels{iSet} = copyobj(PlotHandles.TopoPlot(mod).SensorLabels{1},PlotHandles.TopoPlot(mod).axes{iSet});    
					end
				end
				
				if OPTIONS.ShowContours % Add contour plot
					n = sqrt(size(PlotHandles.TopoPlot(mod).Yo{iSet},1));
					PlotHandles.TopoPlot(mod).Yo{iSet} = reshape(PlotHandles.TopoPlot(mod).Yo{iSet},n,n);
					PlotHandles.TopoPlot(mod).Xo{iSet} = reshape(PlotHandles.TopoPlot(mod).Xo{iSet},n,n);
					[PlotHandles.TopoPlot(mod).Contour{iSet}.C,PlotHandles.TopoPlot(mod).Contour{iSet}.H] = contour(PlotHandles.TopoPlot(mod).Yo{iSet},...
						PlotHandles.TopoPlot(mod).Xo{iSet},...
						CData,'w-');
				end 
				
				
			case 6 % 2D Sensor Cap
				
				if iSet == 1
					chanlocs(:,3)=chanlocs(:,3)-max(chanlocs(:,3));
					[TH,PHI,R]=cart2sph(chanlocs(:,1),chanlocs(:,2),chanlocs(:,3));
					R2=R./cos(PHI).^.2;
					[Y,X]=pol2cart(TH,R2);
                    %Y et X sont les coordonnées projetees dans le plan
                    %'tangent superieur' pour chaque capteur
                    tri = delaunay(Y,X);
                    PlotHandles.TopoPlot(mod).surf{iSet} = patch('faces',tri,'vertices',[Y,X,0.0*ones(size(Y))],'FaceVertexCData',F(OPTIONS.SelectedChannels{mod},tindx),...
                        'edgecolor','none',...
                        'FaceColor','interp');
                    
                    set(PlotHandles.TopoPlot(mod).surf{iSet},'edgecolor','none','facecolor','interp')
                else
                    PlotHandles.TopoPlot(mod).surf{iSet} = copyobj(PlotHandles.TopoPlot(mod).surf{1},PlotHandles.TopoPlot(mod).axes{iSet});
                    set( PlotHandles.TopoPlot(mod).surf{iSet} ,'FaceVertexCData',F(OPTIONS.SelectedChannels{mod},tindx)); 
                end
                
                % Use plot3 here as a trick to let the sensor markers be visible
                if OPTIONS.SensorMarkers
                    if iSet == 1
                        for kk = 1:length(Y) % SB: I don't use scatter here because ButtonDownFcn is not effective as for Matlab 6.1
                            PlotHandles.TopoPlot(mod).hsens{iSet}(kk) = plot3(Y(kk),X(kk),0.05,'.');
                        end
                    else
                        PlotHandles.TopoPlot(mod).hsens{iSet} = copyobj(PlotHandles.TopoPlot(mod).hsens{1},PlotHandles.TopoPlot(mod).axes{iSet});                        
                    end
                end
                if OPTIONS.SensorLabels
                    if iSet == 1
                        for  kk = 1:length(Y) % SB: I don't use scatter here because ButtonDownFcn is not effective as for Matlab 6.1
                            PlotHandles.TopoPlot(mod).SensorLabels{iSet}(kk) = text(Y(kk),X(kk),0.05,...
                                OPTIONS.Channel(OPTIONS.SelectedChannels{mod}(kk)).Name);
                        end
                    else
                        PlotHandles.TopoPlot(mod).SensorLabels{iSet} = copyobj(PlotHandles.TopoPlot(mod).SensorLabels{1},PlotHandles.TopoPlot(mod).axes{iSet});    
                    end
                end
                
                
            case 7 % 2D layout of time series
                
                timescale = 10;
                
                chanlocs = chanlocs'; % Pasted here some older code with different array arrangement for channel locations (CBB) 
                
                minn = min(chanlocs')';
                maxx = max(chanlocs')';
                chanlocs(3,:) =  chanlocs(3,:) - maxx(3);
                
                [TH,PHI,R] = cart2sph(chanlocs(1,:),chanlocs(2,:),chanlocs(3,:));
                PHI2 = zeros(size(PHI));
                R2 = R./cos(PHI).^.2;
                
                [Y,X] = pol2cart(TH,R2);
                dat = F(OPTIONS.SelectedChannels{mod},:);
                M = max(abs(dat(:)));
                
                My = max(abs(Y));
                Mx = max(abs(X));
                Mt = OPTIONS.Time(end)-OPTIONS.Time(1);                
                
                
                for i = 1:length(OPTIONS.SelectedChannels{mod})
                    dat = F(OPTIONS.SelectedChannels{mod}(i),:);
                    PlotHandles.TopoPlot(mod).Lines{iSet}(i) = plot(Mx*[OPTIONS.Time]/(timescale*Mt)-X(i),(My*dat/(10*M)+Y(i)),'color','w');
                    %set(h,'ButtondownFcn','dataplot_cb line_select')
                    hold on
                    ttl = text(Mx*OPTIONS.Time(1)/(timescale*Mt)-X(i),.5*My*M/(10*M)+Y(i),OPTIONS.Channel(OPTIONS.SelectedChannels{mod}(i)).Name);
                    set(ttl,'fontsize',8,'color',.8*[1 1 1])
                end
                axis equal
                axis off
                
                OPTIONS.SensorMarkers = 0;   % Unrelevant for 2D layout
                OPTIONS.SensorLabels = 0;   % Unrelevant for 2D layout
                OPTIONS.Colorbar = 0;    
            end
			
			
			switch(OPTIONS.DisplayTypes(3)) 
			case{2}
				CData = get(PlotHandles.TopoPlot(mod).surf{iSet},'CData');
				CData = CData * max(abs(F(OPTIONS.SelectedChannels{mod},tindx)))/max(abs(CData(:))); % properly rescale to data amplitude
				set(PlotHandles.TopoPlot(mod).surf{iSet},'CData',CData);
			end
			
			if OPTIONS.AbsoluteColormap 
				tmp = max(abs(F(OPTIONS.SelectedChannels{mod},tindx)));
				caxis([-tmp tmp]), clear tmp
			else
				if OPTIONS.DataType == 1
					% 					tmp = [OPTIONS.Fmax{:}]; % Regular surface data topography
					% 					tmp = reshape(tmp,3,nSets)';
					% 					maxData(mod) = max(tmp(:,mod));
					maxData(mod) = max(OPTIONS.Fmax(:,mod));
					
				else OPTIONS.DataType = 2; % Source independent topography
					maxData(mod) = Topo_maxx;
				end
				%caxis([-OPTIONS.Fmax(mod) OPTIONS.Fmax(mod)]),
				caxis([-maxData(mod) maxData(mod)]);
			end
			
			axis equal, axis tight, axis vis3d
			
			xx = get(gca,'Xlim'); 
			yy = get(gca,'Ylim'); 
			if nSets == 1
				sstring = OPTIONS.ModalityLabel{imod};
			elseif nSets > 1
				sstring = {OPTIONS.ModalityLabel{imod},OPTIONS.WaveLabel{iSet}};
			end
			
			text(1.1*yy(2),xx(2),sstring,'FontUnits','Points','Fontsize',fontsize,'HorizontalAlignment','left','Color','w','VerticalAlignment','bottom')
			
			if OPTIONS.SensorMarkers
				if all(ishandle(PlotHandles.TopoPlot(mod).hsens{iSet})) % Sanity check | valid handles
					set(PlotHandles.TopoPlot(mod).hsens{iSet},'Marker','o','markeredgecolor','w','markerfacecolor','none','markersize',3,...
						'ButtonDownFcn','dataplot_cb NameSelectedChannel','selected','off','Tag','Sensor') % Display channel Name when time series is selected
				end
			end
			
			if OPTIONS.SensorLabels
				if all(ishandle(PlotHandles.TopoPlot(mod).SensorLabels{iSet})) % Sanity check | valid handles
					set(PlotHandles.TopoPlot(mod).SensorLabels{iSet},'HorizontalAlignment','center','Fontsize',fontsize,'fontunits','normal')
				end
			end
			
			if OPTIONS.Colorbar % Display Colorbar
				PlotHandles.TopoPlot(mod).Colorbar{iSet} = colorbar('horiz');
				set(PlotHandles.TopoPlot(mod).Colorbar{iSet},'fontsize',fontsize,'fontunits','normal')
				
				tmp = get(PlotHandles.TopoPlot(mod).axes{iSet},'CLim');
				
				if 1%OPTIONS.DataType == 1
					if mod == 1
						tmp = tmp/1e-15; % fT
					else
						tmp = tmp/1e-6; % microV
					end
					%                 elseif OPTIONS.DataType == 2
					%                     tmp = tmp/1e-9; % nAm
				end
				
				
				m = num2str(tmp(2),2);
				[m,e] = strtok(m,'e');
				e = str2num(['1',e]);
				m = str2num(m);
				
				tmp2 = linspace(-m*e,m*e,3)';
				tmp2(abs(tmp2)<0.001) = 0;
				YtickLabel = num2str(tmp2,'%3.0f');
				
				Ytick = linspace(-tmp(2),tmp(2),3)';
				a = get(PlotHandles.TopoPlot(mod).Colorbar{iSet},'Children');
				set(a,'Xdata',tmp)
				
				posSurfAxes = get( PlotHandles.TopoPlot(mod).axes{iSet},'Position');
				posColBarAxes = get( PlotHandles.TopoPlot(mod).Colorbar{iSet},'Position');
				
				set(PlotHandles.TopoPlot(mod).Colorbar{iSet},...
					'Position',[posColBarAxes(1),posColBarAxes(2),posSurfAxes(3),posColBarAxes(4)],...
					'Xaxislocation','bottom','Ticklength', [0 0],'Xlim',tmp,'Xtick',Ytick,'XTickLabel',YtickLabel,...
					'Yaxislocation','right','Ytick',2,'fontunits','points') 
				
				if 1%OPTIONS.DataType == 1
					if mod == 1
						YTickLabel = 'fT';
					else
						YTickLabel = 'mV';
					end
					
				elseif OPTIONS.DataType == 2
					YTickLabel = 'nAm';
				end
				
				set(PlotHandles.TopoPlot(mod).Colorbar{iSet},'YtickLabel',YTickLabel)
				
			end
			
			
			if isfield(PlotHandles.TopoPlot(mod),'surf')
                if ishandle(PlotHandles.TopoPlot(mod).surf{iSet})
                    set(PlotHandles.TopoPlot(mod).surf{iSet},'Visible','on')
                    view(-90,90) % depends on data type and coordinate system
                end
            end
            
			
			PlotHandles.TopoPlot(mod).timeText{iSet} = text(0.05,0.05,0,sprintf('%3.1f ms',1000*OPTIONS.Time(tindx)),'units','normalized');
			set(PlotHandles.TopoPlot(mod).timeText{iSet},'fontweight','normal','color','w','Fontname','helvetica','Fontsize',8,'FontUnits','Point');
			
		end
		
		
	end
	
    if isempty(OPTIONS.Colormap)
    
        % Create Default BrainStorm Colormap 
        Wwidth = 0;  % Width of the white zero level
        Length = 64; % Number of color levels
        Cmin = .2;
        Cmax = .8;
        
        Half = floor(Length/2);    HLength = Half * 2;
        
        R1  = linspace(Cmin,Cmax,Half-Wwidth);
        R2  = ones(1,Wwidth*2);
        R3  = ones(1,Half-Wwidth)*Cmax;
        R   = [R1 R2 R3];
        B   = fliplr(R);
        G3  = fliplr(R1);
        G   = [R1 R2 G3];
        RBWmap = [R' G' B'];
        PlotHandles.Colormap = RBWmap;
        colormap(PlotHandles.Colormap)
        
        OPTIONS.Colormap = PlotHandles.Colormap;
        
    else % User has specified a colormap

        colormap(OPTIONS.Colormap)
        
    end
    
end


switch OPTIONS.MoveCursor
case 'on'
	set(PlotHandles.TimeCursor,'color','r','Tag','TimeCursor','Visible','on','ButtonDownFcn','dataplot_cb MoveTimeCursor start')
	set(PlotHandles.TimeZeroLine,'color',.8*[1 1 1],'Tag','TimeZeroLine','linewidth',1,...
		'linestyle','--','Visible','on')
	set([PlotHandles.TimeCursor],'EraseMode','xor')
end


if nargout == 1
	varargout{1} = PlotHandles;
elseif nargout == 2
	varargout{1} = PlotHandles;
	varargout{2} = OPTIONS;
end

makeuswait('stop')
