Just for variety:
Note also outward normal to surface $z=f(x,y)$: $\nabla(z-f(x,y))=<-f_x,-f_y,1>$ (inward negative)...not same as "NormalVector" function listed (apologies if I have misinterpreted or made error).
I have rewritten taking into account:
f[x_, y_] := Sin[x y];
n[x_, y_] := {-D[f[a, b], a], -D[f[a, b], b], 1} /. {a -> x, b -> y};
va[x_, y_] := VectorAngle[n[x, y], {0, 0, 1}];
Manipulate[
Show[Plot3D[f[x, y], {x, 0, 3}, {y, 0, 3},
ColorFunction ->
Function[{x, y, z},
ColorData["Rainbow"][Rescale[va[x, y], {0, 1.4}]]],
ColorFunctionScaling -> False, PlotRange -> {-2, 2},
PerformanceGoal -> "Quality"],
Graphics3D[{Red, Thick,
Arrow[{{s[[1]], s[[2]],
f[s[[1]], s[[2]]]}, {s[[1]], s[[2]], f[s[[1]], s[[2]]]} +
Normalize@n[s[[1]], s[[2]]]}, 0.03]}],
ContourPlot3D[{x - s[[1]], y - s[[2]], z - f[s[[1]], s[[2]]]}.n[
s[[1]], s[[2]]] == 0, {x, s[[1]] - 0.2,
0.2 + s[[1]]}, {y, -0.2 + s[[2]], 0.2 + s[[2]]}, {z, -2, 2},
Mesh -> False]], {s, {0.1, 0.1}, {3, 3}}]
The graphic illustrates color of surface based on angle between normal and $\vec{k}$ (which seems to be intended color function). I show the unit normal and tangent plane.