function label = zcnt1(rho,theta,k,my_eps)
%
table2 = [7,8,9,0,0,0,9,10,11];

rzc = rho*sin(theta); czc = rho*cos(theta);
zcrr = 2*k(4)+6*k(7)*rzc+2*k(8)*czc;
zcrc = k(5)+2*k(8)*rzc+2*k(9)*czc;
zccc = 2*k(6)+6*k(10)*czc+2*k(9)*rzc;
Hzc = [zcrr zcrc; zcrc zccc];
[Wzc,Lzc] = eig(Hzc);
if abs(Lzc(1,1))<my_eps
  Lzc(1,1) = 0.0;
end;
if abs(Lzc(2,2))<my_eps
  Lzc(2,2) = 0.0;
end;
frzc = k(2)+2*k(4)*rzc+k(5)*czc+3*k(7)*(rzc^2)+2*k(8)*rzc*czc+k(9)*(czc^2);
fczc = k(3)+k(5)*rzc+2*k(6)*czc+k(8)*(rzc^2)+2*k(9)*rzc*czc+3*k(10)*(czc^2);
if norm([frzc,fczc],2)<my_eps
  lindex = 1+(1+sign(Lzc(2,2)))*3+(1+sign(Lzc(1,1)));
  label = table2(lindex);
else
  thetazc = atan2(czc,rzc); sinzc = sin(thetazc); coszc = cos(thetazc);
  Bzc = k(4)*(sinzc^2)+k(5)*sinzc*coszc+k(6)*(coszc^2);
  if abs(Bzc)<my_eps    % 2nd deriv in zc dir is 0
    fzc = k(1)+k(2)*rzc+k(3)*czc+k(4)*(rzc^2)+k(5)*rzc*czc+k(6)*(czc^2)+ ...
          k(7)*(rzc^3)+k(8)*(rzc^2)*czc+k(9)*rzc*(czc^2)+k(10)*(czc^3);
    frc = k(1);
    if fzc>frc
      label = 8;   % ridge
    else
      label = 10;   % ravine
    end;
  elseif Bzc < -my_eps
    label = 8;      % ridge
  else
    label = 10;     % ravine
  end;
end;
