function plotRaster_pub(myTimes, timeLims)
%function [sdfs, sdfTimes, allAxH] = plotRaster_pub(myTimes, timeLims)
% PLOTRASTER_PUB makes raster plot
%
%   PLOTRASTER(MYTIMES) uses times MYTIMES (in sec), can be a vector or a
%       cell
%
%   PLOTRASTER(MYTIMES, TIMELIMS) plot from TIMELIMS(1) to TIMELIMS(2) (in
%       sec)

rasterColor = [0 0 0]; % plot spikes in black

if isempty(myTimes)
    warning('nothing to plot')
    return
end

% initalize
timeEdges = timeLims(1):0.001:timeLims(2);
nTimes = length(timeEdges);

myTimes = myTimes(:);
nTrials = length(myTimes);

% create raster matrix
rasterMat = zeros(nTrials, nTimes);
for irow = 1:size(rasterMat,1)
    if ~isempty(myTimes{irow})
        rasterMat(irow,:) = histc(myTimes{irow}, timeEdges);
    end
end

% compute burst indices
% check if rasters are plotted in red
burstColor = repmat([1 0 0], size(rasterColor, 1), 1);
if any(cellfun(@(x) isequal(x, [1 0 0]), num2cell(rasterColor, 2)))
    warning('plotRaster:labelBursts', 'Burst spikes will be labelled red, and rasterColor is red as well!');
end
dtsilent = 0.100; % 100 ms silence before spike
dtburst  = 0.004; % next spike within 4 ms
iburst = cell(1, nTrials);
itonic = cell(1, nTrials);
for itrial = 1 : nTrials
    [iburst{itrial},itonic{itrial}] = FindBurstSpikes(myTimes{itrial},dtburst,dtsilent);
end
    
% compute burst matrix
burstMat = zeros(nTrials, nTimes);
for irow = 1:size(burstMat,1)
    if isempty(iburst{irow})
        continue
    end
    burstMat(irow,:) = histc(myTimes{irow}(iburst{irow}), timeEdges);
end

% plot the rasters
[trials,timebins] = find(rasterMat); % Find the trial (yPoints) and timebin (xPoints) of each spike
trials = trials(:)';
timebins = timebins(:)';

spikeHeight = 0.5;

% the first two dimensions represent the coordinates for the vertical lines
% the third dimension containing NaNs works as a separator between the
% vertical lines
xPoints = [timebins; timebins; NaN(size(timebins))]*0.001 + timeEdges(1);
yPoints = [trials - spikeHeight; trials + spikeHeight; NaN(size(trials))]*(1/nTrials);
plot(xPoints, yPoints, 'Color',rasterColor);

[trials,timebins] = find(burstMat);
trials = trials';
timebins = timebins';

xPoints = [timebins; timebins; NaN(size(timebins))]*0.001 + timeEdges(1);
yPoints = [trials - spikeHeight; trials + spikeHeight; NaN(size(trials))]*(1/nTrials);
plot(xPoints, yPoints, 'Color', burstColor, 'linewidth', 1);
hold off

allAxH = gca;

spaceHolder = spikeHeight*1/nTrials; % adds space for the first row of spikes
set(allAxH, 'xlim', timeLims, 'ylim', [-inf 1+spaceHolder], 'TickDir', 'out', 'Box', 'off');

xlabel(allAxH(end), 'Time (s)');