
function poincare_map(phases,PRC,T,N)
% Plots the cobweb of the Poincare map
% given the phases and phase resetting curve (PRC)
% T is the forced stimuli period
% N number of iterations of the map
figure;

PTC=mod(phases+PRC+T,phases(end));
hold on;
PTC1=PTC;
PTC1(abs(diff(PTC))>2)=NaN; % skip the discontinuities 
% Set up the plot
axis('square'); 
axis([0 phases(end) 0 phases(end)]);
grid off;
x0=rand(1,1)*phases(end);   % Random initial condition 
x(1) = x0;    
for ic=1:N
 x(ic+1) = interp1(phases,PTC,x(ic)); % PTC is interpolated 
end
% Draw the map
line([x(1) x(1)],[0 x(2)],'Color','k') 
for ic=1:N-1
  line([x(ic) x(ic+1)],[x(ic+1) x(ic+1)],'Color','k')  
  line([x(ic+1) x(ic+1)],[x(ic+1) x(ic+2)],'Color','k')
end
line([x(N) x(N+1)],[x(N+1) x(N+1)],'Color','k') 
plot(phases,PTC1,'b');
plot(phases,phases,'k');
% Draw the spiketrains
figure;
hold on;
s=x0;
k=0;
for ic=1:N-1
   k=k+T;
   s=k+(T-x(ic));
   line([s s],[0.2 0.8]);
   line([k k],[0 0.2],'Color','k','linewidth',2);
end;
line([0 s],[0.3, 0.3]);
line([0 s],[0, 0],'Color','k','linewidth',2);
axis([(s-60) s -0.1 1.5]);
 
 
 