function eta = Estimate_Eigen_2(mesh,Base,lambda,solu,coeff_stiff,coeff_mass,varphi,diff_stiff)


% ESTIMATE_Eigen computes a posteriori error estimator on each element
% for second order elliptic eigenvalue problems with the residual type 
% USAGE
%    eta = Estimate_Eigen(mesh,Base, lambda, solu,coeff_stiff,coeff_mass,varphi,diff_stiff)
%
% INPUT 
%    mesh:  current mesh
%    Base: the base function defined on each element
%    solu:  finite element solution
%    coeff_stiff:  Coefficient function
%    coeff_mass: Mass function
%    varphi:   varphi function in the left hand side 
%
% OUTPUT
%    eta:  error estimator on each element

%--------compute the base value at the integral points---------------------
Nr_elem = size(mesh.elem,1);
%----------compute the freedom degree vector-------------------------------
[Free_deg_3,Dim] = DOF_Mangement(mesh,Base.DOF_3);
%[Direction_3,Direction_4] = Direction_Edge(mesh);
edgenormal = Normal4edge(mesh);
edgelength = Length4edge(mesh);
%full(edgenormal)
Radius=Radius4element(mesh);
%----get the number of triangles-------------------------------------------
II3 = 1:Nr_elem;
% II3 = find(mesh.elem(:,4)<0);
% Nr_elem3 = length(II3);
% II4 = setdiff([1:Nr_elem]',II3);
% Nr_elem4 = length(II4);
%solu = [(mesh.node(:,2)).^3;...
%    ((mesh.node(mesh.node4edge(:,1),2)+mesh.node(mesh.node4edge(:,2),2))/2).^3];

% if Nr_elem4 > 0
%     fprintf('There exist quadratic elements, this is forbidden!\n');
%     return;
% end
Nr_Edge = size(mesh.node4edge,1);
Jump_Edge = zeros(Nr_Edge,1);
eta = zeros(Nr_elem,1);

coord3 = quad_triangle(13);
Nr_int3 = size(coord3,1);
       
coord_edge = [0.5 0.5 1;
              0   0.5 1;
              0.5 0   1];

          
  
 if Nr_elem > 0
     
    m = size(Base.b_3([0,0]),2);
    
    %node4edge = [2,3,1];
    %------compute the edge jump of the a posteriori error estimate--------
    [area,inv_transf] = Transform_information(mesh,1,1);
    
    for kk=1:3
        
        for II = 1: 1
            
            index = kk;
            
            B = Base.gb_3([0 0]*(1-coord_edge(index,1)-coord_edge(index,2))+[1 0]*coord_edge(index,1)+[0 1]*coord_edge(index,2));
            %---------compute the gradient of the approximation------------
            coeff_matrix = coeff_stiff((1-coord_edge(index,1)-coord_edge(index,2))*mesh.node(mesh.elem(II3,1),1:2)+...
                                        coord_edge(index,1)*mesh.node(mesh.elem(II3,2),1:2)+...
                                        coord_edge(index,2)*mesh.node(mesh.elem(II3,3),1:2));
            v = zeros(Nr_elem,2);
            
            for i=1:m                
                v = v + [solu(Free_deg_3(II3,i)),solu(Free_deg_3(II3,i))].*[(inv_transf(II3,1,1)*B(1,i)+inv_transf(II3,2,1)*B(2,i)),...
                                                                            (inv_transf(II3,1,2)*B(1,i)+inv_transf(II3,2,2)*B(2,i))];
  
            end
            v = [v(:,1).*coeff_matrix(:,1,1)+v(:,2).*coeff_matrix(:,1,2),...
                 v(:,1).*coeff_matrix(:,2,1)+v(:,2).*coeff_matrix(:,2,2)];             
            % Time the normal vector
            
            v(:,1) = v(:,1).*edgenormal(II3,2*(kk-1)+1) + v(:,2).*edgenormal(II3,2*(kk-1)+2);
       
            Jump_Edge = Jump_Edge + sparse(mesh.edge(II3,kk),ones(Nr_elem,1),v(:,1),Nr_Edge,1);
            
        end
    end
    
    Jump_Edge = Jump_Edge.*Jump_Edge;  
    
    Jump_Edge(mesh.Dirichlet_edge,:) = 0;

    Jump_Edge =  Jump_Edge.*edgelength.^2;
                   
    %--------delete the Dirichlet boundary egde jump-----------------------
    Jump_Edge(mesh.Dirichlet_edge,1) = 0;
    

    %----compute the residual terms on each element -----------------------
    for kk=1:Nr_int3       
 
        B = Base.b_3([0 0]*(1-coord3(kk,1)-coord3(kk,2))+[1 0]*coord3(kk,1)+[0 1]*coord3(kk,2));
        Grad_B  = Base.gb_3([0 0]*(1-coord3(kk,1)-coord3(kk,2))+[1 0]*coord3(kk,1)+[0 1]*coord3(kk,2));
        Hessian = Base.Hessian_3([0 0]*(1-coord3(kk,1)-coord3(kk,2))+[1 0]*coord3(kk,1)+[0 1]*coord3(kk,2));
        
        v  = zeros(Nr_elem,1);
        W  = zeros(Nr_elem,2,2); 
        HH = zeros(Nr_elem,2,2);
        
        mass_matrix = coeff_mass((1-coord3(kk,1)-coord3(kk,2))*mesh.node(mesh.elem(II3,1),1:2)+...
                coord3(kk,1)*mesh.node(mesh.elem(II3,2),1:2)+...
                coord3(kk,2)*mesh.node(mesh.elem(II3,3),1:2));
            
        varphi_matrix = varphi((1-coord3(kk,1)-coord3(kk,2))*mesh.node(mesh.elem(II3,1),1:2)+...
                coord3(kk,1)*mesh.node(mesh.elem(II3,2),1:2)+...
                coord3(kk,2)*mesh.node(mesh.elem(II3,3),1:2));
            
        coeff_matrix = coeff_stiff((1-coord3(kk,1)-coord3(kk,2))*mesh.node(mesh.elem(II3,1),1:2)+...
                coord3(kk,1)*mesh.node(mesh.elem(II3,2),1:2)+...
                coord3(kk,2)*mesh.node(mesh.elem(II3,3),1:2));
            
         if nargin > 7
             
             diff_stiff_matrix = diff_stiff((1-coord3(kk,1)-coord3(kk,2))*mesh.node(mesh.elem(II3,1),1:2)+...
                                            coord3(kk,1)*mesh.node(mesh.elem(II3,2),1:2)+...
                                            coord3(kk,2)*mesh.node(mesh.elem(II3,3),1:2));
         end
         
         for i = 1 : m
             v = v + (lambda*B(i)*mass_matrix-varphi_matrix*B(i)).*solu(Free_deg_3(:,i));
             %---compute the Hessian matrix in each element ---------------
             %----K^{-T}*H-------------------------------------------------
             w(:,1,1) = [Hessian(1,i)*inv_transf(:,1,1)+Hessian(3,i)*inv_transf(:,2,1)];
             w(:,1,2) = [Hessian(2,i)*inv_transf(:,1,1)+Hessian(4,i)*inv_transf(:,2,1)];
             w(:,2,1) = [Hessian(1,i)*inv_transf(:,1,2)+Hessian(3,i)*inv_transf(:,2,2)];
             w(:,2,2) = [Hessian(2,i)*inv_transf(:,1,2)+Hessian(4,i)*inv_transf(:,2,2)];
             %----K^{-T}*H*K^{-1}------------------------------------------
             HH(:,1,1) = w(:,1,1).*inv_transf(:,1,1)+w(:,1,2).*inv_transf(:,2,1);
             HH(:,1,2) = w(:,1,1).*inv_transf(:,1,2)+w(:,1,2).*inv_transf(:,2,2);
             HH(:,2,1) = w(:,2,1).*inv_transf(:,1,1)+w(:,2,2).*inv_transf(:,2,1);
             HH(:,2,2) = w(:,2,1).*inv_transf(:,1,2)+w(:,2,2).*inv_transf(:,2,2);
             %-------------------------------------------------------------       
             v = v + (coeff_matrix(:,1,1).*HH(:,1,1)+coeff_matrix(:,1,2).*HH(:,1,2)+...
                      coeff_matrix(:,2,1).*HH(:,2,1)+coeff_matrix(:,2,2).*HH(:,2,2)).*solu(Free_deg_3(:,i));
                 
             if nargin > 7
                 %----compute the differential of matrix stiff coefficients
                 % [(A11)_x+(A21)_y, (A12)_x+(A22)_y]*grad(u)
                 w(:,1:2,1) = [solu(Free_deg_3(:,i)),solu(Free_deg_3(:,i))].*[(inv_transf(:,1,1)*Grad_B(1,i)+inv_transf(:,2,1)*Grad_B(2,i)),...
                                                                            (inv_transf(:,1,2)*Grad_B(1,i)+inv_transf(:,2,2)*Grad_B(2,i))];
                 v = v + w(:,1,1).*diff_stiff_matrix(:,1) + w(:,2,1).*diff_stiff_matrix(:,2);                 
             end           
         end
         
        %----compute the residual at the point coord(kk,1:2)            
        eta =  eta + coord3(kk,3)*v.^2.*area.*Radius.^2;       
    end
    %-----add the residual terms and edge jump terms
    eta = eta + Jump_Edge(mesh.edge(II3,1),1)+Jump_Edge(mesh.edge(II3,2),1)+ Jump_Edge(mesh.edge(II3,3),1);
 end