function DiscrCyclIkeda(p)
% Determination of the periodic orbits of the Ikeda map
% Hybrid Newton-Raphson / relaxation methods
% Map x'=f(x) --> Map x'=x+gam*C[f(x)-x]
% INPUT p = period
% OUTPUT periodic orbits with period p
% Last modified: 12/03/2002 (C. Chandre)
global a
a=0.9*exp(i*0.4);

C(:,:,1)=eye(2);
C(:,:,2)=[-1 0;0 1];
C(:,:,3)=[1 0;0 -1];

gam=1e-3; % small parameter

% Set of initial conditions
z0(1)=0;
for it=2:2000
    z0(it)=ikeda(z0(it-1));
end
x0(1,:)=real(z0(1000:10:2000));
x0(2,:)=imag(z0(1000:10:2000));
for ii=1:p
    for ic=1:3
        Xnn=1;
        x=x0(:,ii);
        Cc=C(:,:,ic);
        while Xnn>=1e-16 & Xnn<=1e+3
            xp=modmap(x,p,gam,Cc);
            Xnn=norm(xp-x);
            x=xp;
        end
        if Xnn<=1e-16 
            z=x(1)+i*x(2);
            for it=1:p
                zik(it)=ikeda(z);
                z=zik(it);
            end
            Str=[real(zik)' imag(zik)'];
            periodic=sprintf('%+.10f %+.10f\r',Str')
        end
    end
end

%-----------------------------------
function xp=modmap(x,p,gam,C)
global a
xt=x(1)+i*x(2);
for it=1:p
    xf=ikeda(xt);
    xt=xf;
end
xp=x+gam*C*([real(xf);imag(xf)]-x);

%-----------------------------------
function zik=ikeda(z)
global a
zik=1+a*z*exp(-6*i/(1+abs(z)^2));

%
% Copyright (c) 2002. C. Chandre
%

% Built December 3, 2002
% Comments? email chandre@cpt.univ-mrs.fr