function phi = spikeEodPhase(spike_times, eod_times)
 phi = zeros(length(spike_times), 1);
 period = mean(diff(eod_times));
 for j = 1:length(spike_times)
     t = spike_times(j);
     phase = (t - eod_times(find(eod_times <= t, 1, 'last'))) / period * 2 * pi;
     phi(j) = phase;
 end