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;