function [l1,l2,v1,v2,v3,v4] = Hessianf(k1,k2,k3,k4,k5,k6,k7,k8,k9,k10)
%
% Hessianf: compute eigenvalues and eigenvectors of Hessian matrix using
%           bicubic coefficients k1, k2, ... , k10.
% On input:
%     k1..k10: bicubic coefficients
% On output:
%     l1,l2: eigenvalues
%     v1..v4: eigenvectors
% Call:
%     [l1,l2,v1,v2,v3,v4] = Hessianf(k1,k2,k3,k4,k5,k6,k7,k8,k9,k10)
%

[m,n] = size(k1);
l1 = zeros(m,n);
l2 = zeros(m,n);
v1 = zeros(m,n);
v2 = zeros(m,n);
v3 = zeros(m,n);
v4 = zeros(m,n);

for i1 = 1:m
  for j1 = 1:n
    a = [k4(i1,j1) k5(i1,j1); k5(i1,j1) k6(i1,j1)];
    [vec,lamda] = eig(a);
    l1(i1,j1) = lamda(1,1);
    l2(i1,j1) = lamda(2,2);
    v1(i1,j1) = vec(1,1);
    v2(i1,j1) = vec(2,1);
    v3(i1,j1) = vec(1,2);
    v4(i1,j1) = vec(2,2);
  end
end
