% Antti Hannukainen
%
% 
% Follow the path and look for Hopf / Fold Bif.
%
% 
% [U,Fold_point,Hopf_point]=path_follow(fu,Dfu,u0)
%
% fu  = function handle
% Dfu = function Jacobian handle
% u0  = initial point
%
% U          = the arc
% Fold_point = points, where conditions of Fold are met 
% Hopf_point = points, where conditions of Hopf are met 
%
%

function [U,Fold_point,Fold_indicator,Hopf_point,Hopf_indicator]=path_follow(fu,Dfu,Dxfu,u0,delta,deltaMax,dir,v0)

if( nargin < 4)
  delta = 0.01;
  deltaMax = 1;
 
end

if(nargin < 7)
   dir = 1;
end


U = u0;
u = u0;

Hopf_indicator = [det(bialt(Dxfu(u))) ];
Fold_indicator = [det(Dxfu(u))];

Hopf_point = [];
Fold_point = [];



for i=2:1:ceil(deltaMax/delta)

  uold = u;
  
  if(i==2 & nargin < 8)
    % use dir to determine direction of following !!
    [u,v] = parc(u,fu,Dfu,sign(dir)*delta)
  else
    [u,v] = parc(u,fu,Dfu,delta,v)
  end
  

  % update indicators
  
  Hopf_indicator = [Hopf_indicator det(bialt(Dxfu(u))) ];
  Fold_indicator = [Fold_indicator det(Dxfu(u)) ];
  
  % DETECT HOPF POINT  
  if( Hopf_indicator(i)*Hopf_indicator(i-1) < 0)
    
    % use bisection search
    
    bi_a = uold;
    bi_b = u;
    bi_delta = delta;
    
    for(k=1:50)
      bi_delta = bi_delta/2;
      bi_v = (bi_b - bi_a)/norm(bi_b - bi_a);
      bi_c = parc(bi_a,fu,Dfu,bi_delta,bi_v);
      
      if(  det(bialt(Dxfu(bi_a)))*det(bialt(Dxfu(bi_c))) < 0) 
	bi_b = bi_c;	
	disp('change b');
      else
	
	bi_a = bi_c;
	disp('change a');
      end
     
      if( abs(det(bialt(Dxfu(bi_a)))) < 1e-12)
	Hopf_point = [Hopf_point bi_a ];
	break;
      end
    end
  
  end
  
  % DETECT FOLD POINT (check condition d/ds a = 0)
  
  if( Fold_indicator(i)*Fold_indicator(i-1) < 0)      
     % use bisection search
    
    bi_a = uold;
    bi_b = u;
    bi_delta = delta;
    
    for(k=1:50)
      
      bi_delta = bi_delta/2;
      bi_v = (bi_b - bi_a)/norm(bi_b - bi_a);
      bi_c = parc(bi_a,fu,Dfu,bi_delta,bi_v);
      
      if(  det(Dxfu(bi_a))*det(Dxfu(bi_c)) < 0) 
	bi_b = bi_c;	
	disp('change b');
      else
	
	bi_a = bi_c;
	disp('change a');
      end
     
      if( abs(det(Dxfu(bi_a))) < 1e-12)
	Fold_point = [Fold_point bi_a];
	break;
      end
    end
  end
  

  % update the predictor direction ( use the secant method)
  v = (u-uold)/norm(u-uold);
  U = [U u];  
end