function doall=shilnikov
clear all; close all;
% A return map for a poincare section around a homoclinic orbit of the
% Shilnikov type.  Hopefully showing the horseshoe formed from the
% homoclinic dynamics.  For the course "nonlinear dynamics and chaos",
% fall 2003.  see Guckenheimer & Holmes section 6.5 p 318-323.
% although notations and derivation are a bit different and are from
% Vered Rom Kedar's notes.  Eli Tziperman.

% Two methods are implemented for calculating the map: super imposing
% the maps from Pi0->Pi1 and Pi1->Pi0, or directly.  They don't give
% the same results, so there is still a bug somewhere ...

% Parameters:
global epsilon rho lambda omega a b c d e mu xstar
epsilon=0.01;
rho=-0.8;
lambda=0.9;
omega=4.0;
a=2.0;
b=0.0;
c=0.0;
d=2.0;
e=1.0;
mu=0.0;
xstar=(epsilon/2.0)*(1+exp(rho*2*pi/omega));

% Initialize plotting parameters: 
Nz=50;                                 % number of initial conditions
Nx=50;                                 % number of initial conditions
k=5;                                   % location of i.c. strip
zk  =exp(-2*pi* k   *lambda/omega);    % upper bndry of i.c. strip
zkp1=exp(-2*pi*(k+0.98)*lambda/omega); % lower bndry of i.c. strip
deltaz=(zk-zkp1)/Nz;                    % width of i.c. strip 
MarkerSize=1;

% Initialize plotting arrays:
x_in=zeros(1,(Nx+1)*(Nz+1));
z_in=zeros(1,(Nx+1)*(Nz+1));
x_out=zeros(1,(Nx+1)*(Nz+1));
z_out=zeros(1,(Nx+1)*(Nz+1));

% calculate map of a strip of i.c.:

% (Changing the x-begining of the strip to [zero, larger] makes the
% image [full circle, empty]).
i=0;
for x=epsilon*0.3:epsilon*1.0/Nx:epsilon*2
  for z=zkp1:deltaz:zk
    i=i+1;
    x_in(i)=x;z_in(i)=z;
    [x_out(i),z_out(i)]=P(x,z);
  end
end

% plot i.c. and their mapped image:
axes('position',[0.3 0.3 0.4 0.4]); % [left, bottom, width, height]
h=plot(x_in,z_in,'r.');   set(h,'MarkerSize',MarkerSize);
hold on;
h=plot(x_out,z_out,'b.'); set(h,'MarkerSize',MarkerSize);
h=xlabel('x'); 
  set(h,'FontName','Helvetica','FontWeight','bold','FontSize',12);
h=ylabel('z');
  set(h,'FontName','Helvetica','FontWeight','bold','FontSize',12);
h=title('red: i.c.; blue: M(i.c.)'); 
  set(h,'FontName','Helvetica','FontWeight','bold','FontSize',12);
set(gca,'FontName','Helvetica','FontWeight','bold','FontSize',12);

print -dpsc2 shilnikov.ps


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [r,theta]=P0(x,z)
global epsilon rho lambda omega a b c d e mu xstar
% map from Pi0 to Pi1:
r=x*(epsilon/z)^(rho/lambda);
theta=(omega/lambda)*log(epsilon/z);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x,z]=P1(r,theta)
global epsilon rho lambda omega a b c d e mu xstar
% map from Pi1 to Pi0:
x_in=r*cos(theta);
y_in=r*sin(theta);

x=a*x_in+b*y_in+e*mu+xstar;
z=c*x_in+d*y_in+mu;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x,z]=P(x,z)
global epsilon rho lambda omega a b c d e mu xstar
% from Pi0 to itself, superimposing P0 and P1:

[r,theta]=P0(x,z);
[x,z]=P1(r,theta);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x,z]=P_direct(x,z)
global epsilon rho lambda omega a b c d e mu xstar
% from Pi0 to itself, using vered's solution:

theta=(omega/lambda)*log(epsilon/z);
x=x*(epsilon/z)^(rho/lambda)*(a*cos(theta)+b*sin(theta))+xstar;
z=x*(epsilon/z)^(rho/lambda)*(c*cos(theta)+d*sin(theta));