function spikes = poissonspikes(trials, rate, tmax)
% Generate spike times of a homogeneous poisson process.
%
% spikes = poissonspikes(trials, rate, tmax)
%
% Arguments:
%   trials: number of trials that should be generated
%   rate: the rate of the Poisson process in Hertz
%   tmax: the duration of each trial in seconds
%
% Returns:
%   spikes: a cell array of vectors of spike times in seconds

    dt = 3.33e-5;
    p = rate*dt; % probability of event per bin of width dt
    % make sure p is small enough:
    if p > 0.1
        p = 0.1
        dt = p/rate;
    end
    spikes = cell(trials, 1);
    for k=1:trials
	% uniform random numbers for each bin:
        x = rand(round(tmax/dt), 1);
        spikes{k} = find(x < p) * dt;
    end
end