function out = homoclinic
%
% homoclinic curve definition file for a problem in mapfile
% 
global homds cds
    out{1}  = @curve_func;
    out{2}  = @defaultprocessor;
    out{3}  = @options;
    out{4}  = @jacobian;
    out{5}  = [];%@hessians;
    out{6}  = @testf;
    out{7}  = [];%@userf;
    out{8}  = @process;
    out{9}  = @singmat;
    out{10} = [];%@locate;
    out{11} = @init;
    out{12} = @done;
    out{13} = @adapt;
return


%----------------------------------------------------
function func = curve_func(arg)

  [x,YS,YU,p] = rearr(arg);
  func = BVP_hom(x,YS,YU,p);

%-----------------------------------------------------
function varargout = jacobian(varargin)
  global homds  cds  
  J= homds.niteration;
  [x,YS,YU,p] = rearr(varargin{1});
  varargout{1} = BVP_Hom_jac(homds.func,x,p,YS,YU,J);

function varargout = hessians(varargin)

%------------------------------------------------------

function varargout = defaultprocessor(varargin)
global homds opt cds
    
 % [x,x0,p,T,eps0,eps1,YS,YU] = rearr(varargin{1});
  [x,YS,YU,p] = rearr(varargin{1});
  v = rearr(varargin{2});
  
  %homds.ndim = length(varargin{1});

  % update
%   if ~isempty(homds.ups)
%     homds.upold = homds.ups;
%   end
%  homds.ups = reshape(x,homds.nphase,homds.tps);
%  homds.vps = reshape(v,homds.nphase,homds.tps);
% figure
% plot(homds.ups(1,:),homds.ups(2,:))
  % update upoldp
  p1 = num2cell(p);
%   for i=1:homds.tps
%     homds.upoldp(:,i) = 2*T*feval(homds.func, 0, homds.ups(:,i), p1{:});
%   end
  % homds.eps0 = eps0;
  % homds.eps1 = eps1;
  homds.YS = YS;
  homds.YU = YU;
 

  if nargin > 2
    % set data in special point structure
    s = varargin{3};
    varargout{3} = s;
  end
%   % all done succesfully
   varargout{1} = 0;
   varargout{2} = homds.npoints';
%   
%   if (cds.options.Eigenvalues==1)
%       varargout{2} = [varargout{2}; D];
%   end

%-------------------------------------------------------
  
function option = options
global homds cds
  % Check for symbolic derivatives in odefile
  
  symjac  = ~isempty(homds.Jacobian);
  symhes  = ~isempty(homds.Hessians);
  symder  = ~isempty(homds.Der3);
  
  symord = 0; 
   if symjac, symord = 1; end
   if symhes, symord = 2; end
   if symder, symord = 3; end
  %if higher>2, symord = higher; end

  option = contset;
%   switch homds.nphase
%       case 1
%           option=contset(option,'IgnoreSingularity',[2 3 4]);
%       case 2
%           option=contset(option,'IgnoreSingularity',[4]);
%   end
  option = contset(option, 'SymDerivative', symord);
  option = contset(option, 'Workspace', 1);
  %option = contset(option, 'Locators', zeros(1,13));
   symjacp = ~isempty(homds.JacobianP); 
   symhes  = ~isempty(homds.HessiansP);
   symordp = 0;
   if symjacp, symordp = 1; end
   if symhes,  symordp = 2; end
   option = contset(option, 'SymDerivativeP', symordp);
  
  cds.symjac  = 1;%1
  cds.symhess = 0;
  
%------------------------------------------------------  
function [out, failed] = testf(id, x0, v)
global homds cds        

[x,YS,YU,p] = rearr(x0);
p = n2c(p);
ndim = cds.ndim;
J=contjac(x0);%eig((J(:,1:ndim-1))+eye(ndim-1)),
failed = [];
for i=id
  lastwarn('');
  
  switch i
     
case 1 % LP
     out(1) = v(end);
 case 2 % BP
      B = [J; v'];
      out(2) = det(B);
   
  otherwise
    error('No such testfunction');
  end
  if ~isempty(lastwarn)
    msg = sprintf('Could not evaluate tf %d\n', i);
    failed = [failed i];
  end
  
end

%-------------------------------------------------------------

function [out, failed] = userf(userinf, id, x, v)
global  homds
dim =size(id,2);
failed = [];
[x,x0,p,T,eps0,eps1,YS,YU] = rearr(x); p = num2cell(p);
for i=1:dim
  lastwarn('');
  if (userinf(i).state==1)
      out(i)=feval(homds.user{id(i)},0,x0,p{:});
  else
      out(i)=0;
  end
  if ~isempty(lastwarn)
    msg = sprintf('Could not evaluate userfunction %s\n', id(i).name);
    failed = [failed i];
  end
end

%-----------------------------------------------------------------

function [failed,s] = process(id, x, v, s)
global  cds homds
[x0,YS,YU,p] = rearr(x);
p = n2c(p);
ndim = cds.ndim; 
nphase=homds.nphase;
n=homds.niteration;
 % WM: Removed SL array
fprintf('label = %s, x = ', s.label); printv(x)
p1=p;
 switch id     
  case 1 % LP      
     jac =homjac(x,p,n);
     [V,D]= eig(jac-eye(nphase));
     [Y,i]=min(abs(diag(D)));
     vext=real(V(:,i));
     vext=vext/norm(vext);
     [V,D]= eig(jac'-eye(nphase));
     [Y,i]=min(abs(diag(D)));
     wext=real(V(:,i));      
     wext=wext/(vext'*wext);
     s.msg=sprintf('Limit point\n');
       
  case 2 %BP
      s.msg=sprintf('Branch point\n');  
      s.data.v=v;
end

% Compute eigenvalues for every singularity
J=contjac(x);
if ~issparse(J)
  [v,d]=eig(J(:,1:ndim-1));
else
  opt.disp=0;
  % WM: fixed a bug (incompatability between MatLab 6.0 and 5.5?)
  [v,d]=eigs(J(:,1:ndim-1),min(6,ndim-1),'lm',opt);
end
%d=d+eye(nphase);
s.data.evec = v;
%s.data.eval = diag(d)';

failed = 0;
%-------------------------------------------------------------  

function [S,L] = singmat
global hetds cds
 
% 0: testfunction must vanish
% 1: testfunction must not vanish
% everything else: ignore this testfunction

  S = [  0 8 
         8 1 ] ;

  L = [  'LP  ';'BP  ' ];


  %elseif strcmp(arg, 'locate')


%--------------------------------------------------------

function [x,v] = locate(id, x1, v1, x2, v2)
msg = sprintf('No locator defined for singularity %d', id);
error(msg);
    
%----------------------------------------------------------

function varargout = init(varargin)

WorkspaceInit(varargin{1:2});
% all done succesfully
varargout{1} = 0;

%-----------------------------------------------------------

function varargout = done

%-----------------------------------------------------------

function [res,x,v] = adapt(x,v)
global homds cds

res = []; % no re-evaluations needed


 cds.adapted = 1;
% 
 YU = homds.YU;
 YS = homds.YS;

Q1S = homds.Q1;
QbS1 = Q1S(:,1:homds.nu);
QbS2 = Q1S(:,homds.nu+1:end);
if ~isempty(YS)
    [U1,S1,R1] = svd(QbS1 + QbS2*YS , 0);
    [U2,S2,R2] = svd(QbS2 - QbS1*YS', 0);
    Q1S = [U1*R1', U2*R2'];
else
    Q1S = Q1S;
end

Q0U = homds.Q0;
QbU1 = Q0U(:,1:homds.ns);
QbU2 = Q0U(:,homds.ns+1:end);
if ~isempty(YU)
    [U1,S1,R1] = svd(QbU1 + QbU2*YU , 0);
    [U2,S2,R2] = svd(QbU2 - QbU1*YU', 0);
    Q1U = [U1*R1', U2*R2']; 
else
    Q1U = Q0U;
end
homds.Q0= Q1U;
homds.Q1 = Q1S;
 res = 1;

%----------------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ---------------------------------------------------------------
 
function [x,YS,YU,p] = rearr(x1)
% Rearranges x1 into all of its components
global homds

x = x1(1:homds.nphase*homds.npoints,1);
p = homds.P0;
% eps0 = homds.eps0;
% eps1 = homds.eps1;
idx=homds.npoints*homds.nphase;%+homds.nu.homds.ns;
ju=homds.nphase-homds.nu;
js=homds.nphase-homds.ns;
YU = reshape(x1(idx+1:idx+ju*homds.nu,1),homds.nphase-homds.nu,homds.nu);
idx = idx + ju*homds.nu;
YS = reshape(x1(idx+1:idx+js*homds.ns,1),homds.nphase-homds.ns,homds.ns);
idx = idx + js*homds.ns;
p(homds.ActiveParams) = x1(end,1);

%size(x),pause,YS,YU,eps0,eps1,p,pause
    
% -------------------------------------------------------------

% ---------------------------------------------------------------

function WorkspaceInit(x,v)
global cds homds
% homds.cols_p1 = 1:(homds.ncol+1);
% homds.cols_p1_coords = 1:(homds.ncol+1)*homds.nphase;
% homds.ncol_coord = homds.ncol*homds.nphase;
% homds.col_coords = 1:homds.ncol*homds.nphase;
% homds.coords = 1:homds.ncoords;
% homds.pars = homds.ncoords+(1:2);
% homds.tsts = 1:homds.ntst;
% homds.cols = 1:homds.ncol;
% homds.phases = 1:homds.nphase;
% homds.ntstcol = homds.ntst*homds.ncol;
% 
% homds.idxmat = reshape(fix((1:((homds.ncol+1)*homds.ntst))/(1+1/homds.ncol))+1,homds.ncol+1,homds.ntst);
% homds.dt = homds.msh(homds.tsts+1)-homds.msh(homds.tsts);
% 
% homds.wp = kron(homds.wpvec',eye(homds.nphase));
% homds.pwwt = kron(homds.wt',eye(homds.nphase));
% homds.pwi = homds.wi(ones(1,homds.nphase),:);
% 
% homds.wi = nc_weight(homds.ncol)';
% 
% [homds.bialt_M1,homds.bialt_M2,homds.bialt_M3,homds.bialt_M4]=bialtaa(homds.nphase);

% ------------------------------------------------------

function [x,v,s] = WorkspaceDone(x,v,s)

%------------------------------------------------------------

function K = fastkron(c,p,A,B)
t = p:((c+2)*p-1);
K = A(ones(1,p),fix(t/p)).*B(:,rem(t,p)+1);