function[mas,ev]=maslov(a,ainf,dx)
% This function compute the Maslov index of the unstable space of:
%
%   u_x(x)=A(x)u(x)
%
% where A(i*dx) is given by matrixA(a(i,:)) and 
% lim_{x\to\infty}A(x)=matrixA(ainf)
%
% The routine also needs the following parameters (computed with 
% exterior_algebra.m) to be defined: 
% matrixA         : A function of a which returns A(i*dx) when given as 
%                   a(:,i) as argument.
% compound_matrix : A function of a which returns the n^{th}-compound 
%                   matrix of A(i*dx) when given a(i,:) as argument.
% sigma_matrix    : A matrix such that U\wedge V=U^T sigma_matrix V for U,V 
%                   two n-forms written in the standard basis.
% maslov_form     : A line matrix such that maslov_form*U gives the Maslov 
%                   angle associated to U.
% maslov_matrix   : A matrix such that maslov_matrix*U is the polynomial 
%                   defined by det((V-iW)-lambda(V+iW)), where V and W are
%                             ( V )
%                   such that (   ) has its minor determinants equal to U.
%                             ( W )

global sigma_matrix
global maslov_form
global maslov_matrix

N=size(compound_matrix(ainf),1);
AINF=matrixA(ainf);

skip=1;
dx=skip*dx;
if (length(find(abs(real(eig(AINF)))<10^(-10)))==0)


    n=size(a,1)-1;
    BB=symbase();

    e=BB*[1;zeros(N-1,1)];
    d=BB*[zeros(N-1,1);1];
    d=d/(transpose(e)*(sigma_matrix*d));
    s=(e'*compound_matrix(ainf)*e)/(e'*e);

    maslov_form2=maslov_form*BB^(-1);
    maslov_matrix2=maslov_matrix*BB^(-1);
    mas=0;
    

    for ii=1:2*skip:(n-1)
        ep=e;
        inter=transpose(e)*sigma_matrix*d;

        % Fourth-order Runge-Kutta method integration step.
        AA=compound_matrix(a(ii,:));
        f1=AA*e-s*e;

        AA=compound_matrix(a(ii+skip,:));
        aux=e+dx*f1;
        f2=AA*aux-s*aux;

        AA=compound_matrix(a(ii+skip,:));
        aux=e+dx*f2;
        f3=AA*aux-s*aux;

        AA=compound_matrix(a(ii+2*skip,:));
        aux=e+2*dx*f3;
        f4=AA*aux-s*aux;

        e=e+(2*dx/6)*(f1+2*f2+2*f3+f4);

        % Computation of the Maslov angle.
        deltamas=angle(maslov_form2*e*exp(-i*mas));
        if abs(deltamas)>.1
            disp('Warning: the lift might not be computed accurately')
        end
        mas=mas+deltamas;
        
%            roots(maslov_matrix2*e)
%            err=[norm(abs(roots(maslov_matrix*e))-1),err];
%            prod(roots(maslov_matrix2*e))


    end
%   semilogy(err);drawnow;
    ev=(transpose(e)*(sigma_matrix*d));
    mas=mas/(pi);
else
    disp('This point is in the continuous spectrum')
    mas=NaN;
    ev=NaN;
end


%Computation of a basis adpated to the stable and the unstable space.
    function[y]=symbase()
        nn=size(AINF,1);

        J=zeros(nn,nn);
        J(nn/2+1:nn,1:nn/2)=speye(nn/2,nn/2);
        J(1:nn/2,nn/2+1:nn)=-speye(nn/2,nn/2);

        [v,e]=eig(AINF);
        e=diag(e);
        I1=find(real(e)>0 & imag(e)>0);
        I2=find(real(e)>0 & imag(e)==0);
        U=[real(v(:,I2)),real(v(:,I1)),imag(v(:,I1))];

        I1=find(real(e)<0 & imag(e)>0);
        I2=find(real(e)<0 & imag(e)==0);
        S=[real(v(:,I2)),real(v(:,I1)),imag(v(:,I1))];



        No=(U'*J*S)^(-1);
        V=[U*No',S];
        y=zeros(N,N);
        subsets=nchoosek(1:nn,nn/2);

        for r=1:N
            for s=1:N
                y(r,s)=det(V(subsets(r,:),subsets(s,:)));
            end
        end
    end

end
