function n = normalVec(T, i, j, m, n)

global HUGE

% compute surface normal to T at (i, j)
if i < m & T(i+1, j) ~= HUGE
    dYfrw = T(i+1, j) - T(i, j);
else
    dYfrw = 0;
end
if j < n & T(i, j+1) ~= HUGE
    dXfrw = T(i, j+1) - T(i, j);
else
    dXfrw = 0;
end
if i > 1 & T(i-1, j) ~= HUGE
    dYbck = T(i, j) - T(i-1, j);
else
    dYbck = 0;
end
if j > 1 & T(i, j-1) ~= HUGE
    dXbck = T(i, j) - T(i, j-1);
else
    dXbck = 0;
end
n = [0 0];
if dYfrw ~= 0
    if dXfrw ~= 0
        n1 = [dYfrw dXfrw];
        len = norm(n1);
        if len ~= 0 
            n = n1/len;
        end
    end
    if dXbck ~= 0
        n1 = [dYfrw dXbck];
        len = norm(n1);
        if len ~= 0 
            n = n + n1/len;
        end
    end
end
if dYbck ~= 0
    if dXfrw ~= 0
        n1 = [dYbck dXfrw];
        len = norm(n1);
        if len ~= 0 
            n = n + n1/len;
        end
    end
    if dXbck ~= 0
        n1 = [dYbck dXbck];
        len = norm(n1);
        if len ~= 0 
            n = n + n1/len;
        end
    else
        if n == [0 0]
            n = [0 1];  % [0 dXbck] normalized 
        end
    end
else
    if n == [0 0]
        if (dXfrw ~= 0) | (dXbck ~= 0)
            n = [1 0];
        end
    end
end
