LambertW函数

if abs(real(z))+abs(imag(z))<1e-10
      wx=z;
    return ;
    end
    lnz=log(z);
    zx=real(lnz);
    zy=imag(lnz);
 
    Fy=zeros(2);
    F=Fy;
    Fx=Fy;
    temp=0;
    x=zx;
    y=zy;
    x0=0;
    y0=0;
    Fy(1)=log(x^2+y^2)/2+x-zx;
    Fy(2)=y-zy+atan2(y,x);
    erro=abs(Fy(1))+abs(Fy(2));
    loopn=1000;
    w=1;
    while loopn>0 && w>-1.05 && erro>1e-10
        w=x^2+y^2;
        F(1)=x/w+1;
        F(2)=y/w;
        w=F(1)^2+F(2)^2;
        Fx(1)=(Fy(1)*F(1)-Fy(2)*F(2))/w;
        Fx(2)=(Fy(2)*F(1)+Fy(1)*F(2))/w;
        w=1;
        while w>=-1
            x0=x-w*Fx(1);
            y0=y-w*Fx(2);
            Fy(1)=log(x0^2+y0^2)/2+x0-zx;
            Fy(2)=y0-zy+atan2(y0,x0);
            temp=abs(Fy(1))+abs(Fy(2));
            if temp<erro
                erro=temp;
                x=x0;
                y=y0;
                w=-1.03;
            else
                 w=w-0.1;
            end
        end
        loopn=loopn-1;
    end
    wx=x+i*y;
    end