function W=lambertw(t)
% Lambert's W function satisfies t=W*exp(W) when W=lambertw(t)
% This implementation requires that t>=0.

if any(t<0 || ~isreal(t)),
  disp('lambertw: argument must be real non-negative');
  return;
  end
  
% Initial values
W=log(1+t); 
% Improve initial values for large W
k=find(W>1); W(k)=W(k)-log(W(k)); 

% Newton's method is quadratically convergent; can use coarse tolerance
tol=1e-10;
while 1,
  dW=(W-t.*exp(-W))./(1+W); W=W-dW;
  if all(abs(dW)<=tol*W) break; end
end