function [W,iter]=lambertw(t,tol,display)
% function W=lambertw(t)
%
% delivers Lambert's W function defined by
% W*exp(W)=t for t>=zero.
%
% t         vector of arguments
% tol       relative tolerance to stop iteration
% display   'on': plot some interesting values
%
% W         vector of solutions
% iter      number of necessary newton steps


% positivity check
if any(t<0) || any(~isreal(t))
    disp('lambertw: argument must be real and non-negative');
    return;
end

% initial values
W=log(1+t);

% improve initial values for large W using fixpoint iteration
% one step is sufficient
k=find(W>1);
W(k)=log(t(k))-log(W(k));

% apply Newton's method
iter=0;
while 1
    dW=(W-t.*exp(-W))./(1+W);
    W=W-dW;
    iter=iter+1;
    if all(abs(dW)<=tol*W) 
        break;
    end
end

% display
if isequal(display,'on')
    fprintf('lambertw: zero W=%20.16e\n',W);
    fprintf('lambertw: newton steps=%d\n',iter);
end

return