function [theUnStableManif,theStableManif]=invtManifs(n,npts)

% invtManifs(n,npts) draws the unstable and the stable manifolds of the
% Henon map by iterating a line segment (npt points) n times.  
% (Default n=10, npts=10000.)
% It returns no results.


% Copyright (c) 1999 by Alistair Mees.  
%
% Copy freely but do not sell.
%
%   $Id$

if nargin>=3 | nargout >=3
  error 'This function requires no outputs and at most two inputs.'
end;
if nargin<1
  n = 10;
end
if nargin<2
  npts = 10000;
end

% The Henon parameters
global a b
a = 1.4;
b = 0.3;

% Define lots of points to represent a straight line segment.
s = [0:npts]/npts;

% The end points of the segment:
% First the fixed pt is (xstar,xstar) where
xstar = (b-1+sqrt((b-1)^2+4*a))/2;
% Let's have a line a distance delta or thereabouts to the right and up:
delta = 0.1;
x0 = [xstar-delta; xstar+delta];
x1 = [xstar+2*delta; xstar-delta];
% And another approx orthog to that for the unstable manifold:
delta2=1e-4;
y0 = [xstar-delta2; xstar-delta2];
y1 = [xstar+delta2; xstar+delta2];

% The segment:
unStableManif = [(1-s)*x0(1)+s*x1(1); (1-s)*x0(2)+s*x1(2)];
stableManif   = [(1-s)*y0(1)+s*y1(1); (1-s)*y0(2)+s*y1(2)];

for i=1:n
  unStableManif = henon(unStableManif);
  stableManif = inverseHenon(stableManif);
end
% Add some more points:
%unStableManif = [unStableManif,henon(unStableManif)];
%stableManif = [stableManif,inverseHenon(stableManif)];

% Plot it:

close all
figure
%whitebg('k');
unStableHandle = plot(unStableManif(1,:),unStableManif(2,:),'r');
axis(axis);  % fix the axes
drawnow
hold on
stableHandle = plot(stableManif(1,:),stableManif(2,:),'b');
title 'Unstable (red) and stable (blue) manifolds'

zoom on
drawnow


if nargout>0
  theUnstableManif = unstableManif;
  if nargout>1
    theStableManif = stableNamif;
  end
end




function xvecNew = henon(xvec)
% Vectorised henon, i.e. return images of lots of points.
% xvec is a 2 by m matrix, representing m points.  Returns a 2 by m matrix.
%
% Form used is f(x,y) = (a-x^2+b*y,x)

global a b

xvecNew = zeros(size(xvec));

xvecNew(1,:) = a - xvec(1,:).^2 + b*xvec(2,:);
xvecNew(2,:) = xvec(1,:);



function xvecThen = inverseHenon(xvec)
% Inverse is finv(x,y) = (y,(x-a+y^2)/b)

global a b

xvecThen = zeros(size(xvec));
xvecThen(1,:) = xvec(2,:);
xvecThen(2,:) = (xvec(1,:) - a + xvec(2,:).^2)/b;