function value = heaviside( z, epsilon )
% HEAVISIDE Heaviside function of z
%    HEAVISIDE( z, epsilon ) Heaviside function of z, regularized
%    based on epsilon.

value = 0.5 .* ( 1 + ( 2 ./ pi ) .* atan( z ./ epsilon ) );