% Antti Hannukainen 2.3.2006
%
% Compute periodic orbits of the  Lorenz - system
%
%
%
%
% 


% Define function and it's derivatives
s = 10;
b = 8/3;


f   = @(x,a)( [ -s*x(1)+s*x(2);  
	         a*x(1)-x(1)*x(3)-x(2);
	         x(1)*x(2) - b*x(3)]);

Dxf = @(x,a)( [ -s    s     0;
		a-x(3)  -1    -x(1);
		x(2)    x(1)  -b]);

Daf =  @(x,a)( [ 0 ; x(1) ; 0] );


% Define the system needed in standard path following

fu  = @(u)( f(u(1:3),u(4)));
Dfu = @(u)( [Dxf(u(1:3),u(4)) Daf(u(1:3),u(4))] );
Dxfu = @(u)( Dxf(u(1:3),u(4)) );


x0 = [ 0 ; 0 ; 0];
a  = 0;

u0 = [x0;a];


% First find the "pitchfork" (Look's like Fold)


% Take the "pitchfork" and follow other curves 
% Build new direction for the curve following !

% Follow both branches
% Branch 1

% And Branch 2




% Plot the situation
figure;

plot3(B1_U(1,:),B1_U(2,:),B1_U(3,:));

hold on;
plot3(B2_U(1,:),B2_U(2,:),B2_U(3,:));

% Now Both Hopf's are found => compute periodic solutions !
% 
% Periodic solutions are computed with the same continuation as 
% other solutions.

% build mesh in time

t = 0:0.025:1;

% Build  initial guess ( compute eigenvalues of hopf + use ideas of reduction )
%
%

%-------------------------------------
%
% Branch 1
%
%-------------------------------------
[X,D] = eig(Dxfu(B1_Hopf_point));

w = abs(D(2,2));

u = X(:,2);
v = X(:,3);

z  = kron(sin(w*t)',real(u) ) + kron(cos(w*t)',real(v)) ;

Del = 0.01;

old_u = zeros(3*length(t)+2,1);


per_u = [ Del*z + repmat(B1_Hopf_point(1:3),length(t),1) ; 2*pi/w ; B1_Hopf_point(4)];


% Take the first step with gauss-Newton
% 
% 

u = gnewton(per_u, @(u)(abu_f(u,old_u,t,f)), @(u)(abu_Df(u,old_u,t,f)),@pifun);

xxx = bvp_u2x(u(1:(length(u)-2)),length(t) );
plot3(xxx(1,:),xxx(2,:),xxx(3,:),'r*');


% Do the continuation with parc !

[u,v] = parc(u ,@(u)(abu_f(u,per_u,t,f)), @(u)(abu_Df(u,per_u,t,f)),0.0001);

for i=1:1:20

  % Take very long steps
  [u,v] = parc(u ,@(u)(abu_f(u,per_u,t,f)), @(u)(abu_Df(u,per_u,t,f)),10,v);
  
  

  hold on;
  xxx = bvp_u2x(u(1:(length(u)-2)),length(t) );
  plot3(xxx(1,:),xxx(2,:),xxx(3,:),'r*-');
 
end

%-------------------------------------
%
% Branch 2
%
%-------------------------------------

% Build  initial guess ( compute eigenvalues of hopf + use ideas of reduction )
%
%

[X,D] = eig(Dxfu(B2_Hopf_point));

w = abs(D(2,2));

u = X(:,2);
v = X(:,3);

z  = kron(sin(w*t)',real(u) ) + kron(cos(w*t)',real(v)) ;

Del = 0.01;

old_u = zeros(3*length(t)+2,1);


per_u = [ Del*z + repmat(B2_Hopf_point(1:3),length(t),1) ; 2*pi/w ; B2_Hopf_point(4)];


% Take the first step with gauss-Newton
% 
% 

u = gnewton(per_u, @(u)(abu_f(u,old_u,t,f)), @(u)(abu_Df(u,old_u,t,f)),@pifun);

xxx = bvp_u2x(u(1:(length(u)-2)),length(t) );
plot3(xxx(1,:),xxx(2,:),xxx(3,:),'r*');


% Do the continuation with parc !

[u,v] = parc(u ,@(u)(abu_f(u,per_u,t,f)), @(u)(abu_Df(u,per_u,t,f)),0.0001);

for i=1:1:20

  % Take very long steps
  [u,v] = parc(u ,@(u)(abu_f(u,per_u,t,f)), @(u)(abu_Df(u,per_u,t,f)),10,v);
  
  

  hold on;
  xxx = bvp_u2x(u(1:(length(u)-2)),length(t) );
  plot3(xxx(1,:),xxx(2,:),xxx(3,:),'r*-');
 
end

