function out = homoclinicT
%
% homoclinic tangency 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)
global cds homTds

  [x,YS,YU,p] = rearr(arg);
  J=homTds.Niterations;
  n1=homTds.nphase;
  n2= homTds.npoints;
  nu=homTds.nu;
  ns=homTds.ns;
  b=homTds.b;
  c=homTds.c;
  K1=n1*(n2-1)+nu*(n1-nu)+ns*(n1-ns)+2*n1-nu-ns;
  jac=BVP_homT_jac(x,p,YS,YU,J);
  Bord=[jac b;c' 0];
  bunit=[zeros(K1,1);1];
  v=Bord\bunit;
  g=Bord'\bunit;
  f = BVP_homT(x,YS,YU,p);
  func = [f ; v(end)];   
  
%-----------------------------------------------------
function varargout = jacobian(varargin)
  global homTds  cds  
  
%-----------------------------------------------------

function varargout = hessians(varargin)

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

function varargout = defaultprocessor(varargin)
global homTds opt cds
  %1,size(varargin{1})
  %2,size(varargin{2})
 % [x,YS,YU,p] = rearr(varargin{1});
  %v = rearr(varargin{2});  
 
  %p1 = num2cell(p);
  %homTds.YS = YS;
  %homTds.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} = homTds.npoints';

%-------------------------------------------------------
  
function option = options
global homTds cds
  % Check for symbolic derivatives in mapfile
% homTds,pause, homTds.Jacobian,pause
     symjac  = ~isempty(homTds.Jacobian);
     symhes  = ~isempty(homTds.Hessians);
     symder  = ~isempty(homTds.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;
 option = contset(option, 'SymDerivative', symord);
 %option = contset(option, 'Workspace', 1);
%  option = contset(option, 'Locators', zeros(1,13));
  symjacp = ~isempty(homTds.JacobianP); 
%   symhes  = ~isempty(homTds.HessiansP);
%   symordp = 0;
   if symjacp, symordp = 1; end
   if symhes,  symordp = 2; end
   option = contset(option, 'SymDerivativeP', symordp);
  
  cds.symjac  = 0;
  cds.symhess = 1;
  
%------------------------------------------------------  
  
function [out, failed] = testf(id, x0, v)
global homTds 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  homTds cds
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(homTds.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 homTds
[x0,YS,YU,p] = rearr(x);
p = n2c(p);
ndim = cds.ndim; 
nphase=homTds.nphase;
n=homTds.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 homTds 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 homTds cds

res = []; % no re-evaluations needed
[x1,YS,YU,p] = rearr(x);

cds.adapted = 1;
% 
J=homTds.Niterations;
p=homTds.P0;
YS=homTds.YS;
YU=homTds.YU;
jac =BVP_homT_jac(x1,p,YS,YU,J);
[U,S,V]=svd(full(jac));
homTds.b=U(:,end);
homTds.c=V(:,end);

res = 1;


%homTds.Q0,homTds.Q1,pause
%----------------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ---------------------------------------------------------------
 
function [x,YS,YU,p] = rearr(x1)
% Rearranges x1 into all of its components
global homTds

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

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

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

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

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

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

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