# Maple example to find normal vector on a height field with regularly spaced
# values in x and y (separated by d)

with(linalg);

# Compute the normal vector at point p; here is point p and its
# adjacent 4 vertices
#
#        pt
#         |
#   pl -- p -- pr
#         |
#        pb

p  := [px,   py,   pz];
pl := [px-d, py,   plz];
pr := [px+d, py,   prz];
pt := [px,   py+d, ptz];
pb := [px,   py-d, pbz];

# 4 vectors along edges at point p (pointing outward from p)
vl := evalm(pl-p);
vr := evalm(pr-p);
vt := evalm(pt-p);
vb := evalm(pb-p);

# compute smoothed normal vector by averaging normals at each corner
n := simplify(evalm(crossprod(vt,vl) + crossprod(vl,vb) + crossprod(vb,vr) + crossprod(vr,vt)));

# we can divide out a common factor as we're normalizing this...
simplify(evalm(n/(2*d)));

# take note how the answer doesn't depend on the location of p = (px,py,pz)
# but rather just the differences in height in each direction, and the
# separation distance betweeen the vertices