function a=kermat(x,y,c)
f=exp(-(x.^2+y.^2)/c^2);
a= [ -2*f.*(-c^2+2*y.^2)/c^4,           4*x/c^4.*y.*f;...
        4*x/c^4.*y.*f,         -2*f.*(-c^2+2*x.^2)/c^4];
% multiquadric
% nu =[ -c^2-x.^2,     -x.*y; -x.*y, -c^2-y.^2];
% de =repmat((c^2+x.^2.+y.^2).*c^2.*((c^2+x.^2.+y.^2)/c^2).^(1/2),2,2);
% a=nu./de;
return