function [phi,beats]=poincare(phizero,b,tau,niter)

%   [PHI,BEATS] = POINCARE(PHIZERO,B,TAU,NITER) iterates the 1D Poincare 
%   oscillator. PHIZERO is the initial phase, B is the stimulation strength,
%   TAU is the period of the stimulation, and NITER is the number of 
%   iterations.
%   The output consists of two arrays: 
%   PHI is a listing of the succesive phases during periodic stimulation, 
%   and BEATS is a listing of the number of beats that occur between
%   successive stimuli.

%   Copyright 2003
%   Centre for Nonlinear Dynamics in Physiology and Medicine

phi=zeros(1,niter);
phi(1)=phizero;

for i=2:niter
 angle= 2*pi*phi(i-1);
 rprime=sqrt(1+b^2+2*b*cos(angle));
 argument=(cos(angle)+b)/rprime;
 phi(i)=acos(argument)/(2*pi);
 if phi(i-1) > 0.5
  phi(i)=1-phi(i);
 end;
phi(i)=phi(i)+tau;
beats(i)=phi(i)-rem(phi(i),1);
phi(i)=rem(phi(i),1);
end;
end;


