% Numerically calculate local value of hx from average value of h 
function hx = hLoc(flag,hfun,x)
  span=[0.999; 1.001]*x(:)';     % numerical step to calculate slopes
  if flag=='T'
    hL=@(L)hfun(L).*L;
    hx=(hL(span(2,:))-hL(span(1,:)))./diff(span);
  elseif flag=='H'
    Lh=@(L)L./hfun(L);
    hx=diff(span)./(Lh(span(2,:))-Lh(span(1,:)));
  else
    error('Unknown flag, ''T'' or ''H'' required');
  end
end
