function out = heteroclinic
%
% heteroclinic curve definition file for a problem in mapfile
% 
global hetmds 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,eps1,p] = rearr(arg);
  func = BVP_Het(x,YS,YU,eps1,p);
 
  
%-----------------------------------------------------

function varargout = hessians(varargin)

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

function varargout = defaultprocessor(varargin)
global homds opt cds
    
%   [x,x0,p,T,eps0,eps1,YS,YU] = 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;
%   homds.T = T;
%   homds.x0 = x0;
%   
% % Update dimensions
% % -----------------
% p = num2cell(p);
% A = Hom_odejac(x0,p);
% D = eig(A);
% % nneg = dimension of stable subspace
% homds.nneg = sum(real(D) < 0);
% 
% % If one eigenvalue is (practically) zero, and the one of the subspaces has
% % zero dimension, change this dimension with 1.
% if (homds.nneg == homds.nphase)
%     if min(abs(real(D))) < 1e-2
%         homds.nneg = homds.nneg -1;
% %     else
% %         error('This is an HSN.');
%     end
% end
% if (homds.nneg == 0)
%     if min(abs(real(D))) < 1e-2
%         homds.nneg = homds.nneg +1;
% %     else
% %         error('This is an HSN.');
%     end
% end
% homds.npos = homds.nphase-homds.nneg;
% homds.Ysize = homds.nneg*homds.npos;
%   
%   if nargin > 2
%     % set data in special point structure
%     s = varargin{3};
%     s.data.timemesh = homds.msh;
%     s.data.ntst = homds.ntst;
%     s.data.ncol = homds.ncol;
%     s.data.parametervalues = p;
%     s.data.T = T;
%     varargout{3} = s;
%   end
%   % all done succesfully
%   varargout{1} = 0;
%   varargout{2} = homds.msh';
%   
%   if (cds.options.Eigenvalues==1)
%       varargout{2} = [varargout{2}; D];
%   end

%-------------------------------------------------------
  
function option = options
global hetds cds
  % Check for symbolic derivatives in odefile
  
%  symjac  = ~isempty(hetds.Jacobian);
%  symhes  = ~isempty(hetds.Hessians);
%  symder  = ~isempty(hetds.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(hetds.JacobianP); 
%   symhes  = ~isempty(hetds.HessiansP);
%   symordp = 0;
%   if symjacp, symordp = 1; end
%   if symhes,  symordp = 2; end
%   option = contset(option, 'SymDerivativeP', symordp);
  
  cds.symjac  = 0;
  cds.symhess = 0;
  
%------------------------------------------------------  
  
function [out, failed] = testf(id, x0, v)
global homds cds BigJac OldBigY

if isempty(BigJac)
    tmpt = jacobian(x0);
end
        
[x,x0,p,T,eps0,eps1,YS,YU] = rearr(x0);
ups = reshape(x,homds.nphase,homds.tps);

A = Hom_odejac(x0,num2cell(p));
D = eig(A);
nneg = sum(real(D) < 0);

[D1,indlist] = sort(real(D));
[val, ind] = min(abs(real(D1)));

if real(D1(ind)) > 0
    lambda1 = ind;
    lambda2 = ind+1;
    lambda3 = ind+2;
    mu1 = ind-1;
    mu2 = ind-2;
    mu3 = ind-3;
else
    lambda1 = ind+1;
    lambda2 = ind+2;
    lambda3 = ind+3;
    mu1 = ind;
    mu2 = ind-1;
    mu3 = ind-2;
end

D2 = D1;
D2(ind) = [];
[val2,ind2] = min(abs(real(D2)));
if ind2 >= ind
    ind2 = ind2 + 1;
end

res = zeros(14,1);
out = res;
failed = [];
% return

% at least 1 negative and 1 positive eigenvalue
if (mu1 > 0) & (lambda1 <= homds.nphase)
    % 1. Neutral saddle, saddle-focus or bi-focus
    res(1,1) = real(D1(mu1)) + real(D1(lambda1));

    % 2 negative eigenvales and 1 positive
    if (mu2 > 0)   
        % 2. Double real stable leading eigenvalues
        if imag(D1(mu1)) < cds.options.VarTolerance
            res(2,1) = (real(D1(mu1)) - real(D1(mu2)))^2;
        else
            res(2,1) = -(imag(D1(mu1)) - imag(D1(mu2)))^2;
        end
        
        % 4. Neutrally-divergent saddle-focus (stable)
        res(4,1) = real(mu1) + real(mu2) + real(lambda1);
        
        if (mu3 > 0)
            % 6. Three leading eigenvalues (stable)
            res(6,1) = real(mu1) - real(mu3);
        end
    end
    
    if (lambda2 <= homds.nphase)
        % 3. Double real unstable leading eigenvalues
        if imag(D1(lambda1)) < cds.options.VarTolerance
            res(3,1) = (real(D1(lambda1)) - real(D1(lambda2)))^2;
        else
            res(3,1) = -(imag(D1(lambda1)) - imag(D1(lambda2)))^2;
        end
        
        % 5. Neutrally-divergent saddle-focus (unstable)
        res(5,1) = real(mu1) + real(lambda2) + real(lambda1);
        
        if (lambda3 > 0)
            % 7. Three leading eigenvalues (unstable)
            res(7,1) = real(mu1) - real(mu3);
        end
    end
elseif mu1 == 0
    res(1,1) = D1(lambda1);
    res(4,1) = real(D1(lambda1));
else
    res(1,1) = D1(mu1);
    res(4,1) = real(D1(mu1));
end

% 8. Eigenvalue with zero real part => non-hyperbolic equilibrium
if nneg == homds.nneg
    % The signs of the eigenvalues are still all the same, if an
    % eigenvalue is small enough, we have a zero real part.    
    res(8,1) = sign(real(D1(ind)))*real(D1(ind)) - 10*cds.options.VarTolerance;
    
else
    % One of the negative eigenvalues has turned positive
    res(8,1) = real(D1(ind));
end

% 10. Non-central HSN: The eigenvalue with smallest real part has zero imaginary part
if abs(imag(D1(ind))) < cds.options.VarTolerance
    res(9,1) = res(8,1);
else
    res(9,1) = imag(D1(ind));
end

% 11. Bogdanov-Takens point: 2nd smallest eigenvalue must be zero
%prod = D1(mu1) * D1(lambda1);
%res(11,1) = sign(prod)*prod - cds.options.VarTolerance;
%res(11,1) = sign(real(D1(ind2)))*real(D1(ind2)) - cds.options.VarTolerance;
res(10,1) = norm(D1(ind2)) - cds.options.VarTolerance;

if (D1(1) > 0) || (D1(end) < 0)
    % This is a saddle-node, not a saddle!, so BT, NCH or something should
    % have been detected
    res(11:14,1) = 0;
    return;
end

V = [];
if (mu1 > 0)
if abs(imag(D1(mu1))) < cds.options.VarTolerance
    % 11. Orbit-flip with respect to stable manifold
    [V,D] = eig(A');
    % We still have the ordering of the eigenvalues of A (and thus A') from
    % before
    VN = V(:,indlist);
    w1s = VN(:,mu1);
    res(11,1) = exp(-D1(mu1) * homds.T) * (w1s' * (ups(:,end) - x0));
end
end

if (lambda1 <= homds.nphase)
if abs(imag(D1(lambda1))) < cds.options.VarTolerance
    if isempty(V)
        [V,D] = eig(A');
        VN = V(:,indlist);
    end
    w1u = VN(:,lambda1);
    % 12. Orbit-flip with respect to unstable manifold
    res(12,1) = exp(D1(lambda1) * homds.T) * (w1u' * (ups(:,1) - x0));
end

end

V = [];
if (mu1 > 0) & (~isempty(YS))
    if abs(imag(D1(mu1))) < cds.options.VarTolerance
        % 13. Inclination-flip with respect to stable manifold
        [V,D] = eig(A);
        VN = V(:,indlist);
        if homds.nneg == 0
            homds.nneg = homds.nneg + 1;
        elseif homds.nneg == homds.nphase
            homds.nneg = homds.nneg - 1;
        end
        
        BigJac2 = BigJac([1:homds.ncoords-homds.nphase end-1-homds.nphase:end-2],1:homds.ncoords);
        
        if (isempty(OldBigY)) | (size(OldBigY,2) < 2) | (size(OldBigY,1) ~= homds.ncoords)
            psiold = rand(homds.ncoords,1);
            phiold = rand(homds.ncoords,1);
            BigJacT = [BigJac2 psiold; phiold' 0];
            while condest(BigJacT) > 1e10
                psiold = rand(homds.ncoords,1);
                phiold = rand(homds.ncoords,1);
                BigJacT = [BigJac2 psiold; phiold' 0];
            end    
        else
            psiold = OldBigY(:,1);
            phiold = OldBigY(:,2);  
            BigJacT = [BigJac2 psiold; phiold' 0];