Is Lambert W Function in complex domain $(x+iy)e^{x+iy}=a+ib$ solvable? I am solving by method of approximation complex Lambert $W$ Function $(x+iy)e^{x+iy}=a+ib$ for $x$ and $y$, when real values of $a$ and $b$ are given and $i=\sqrt{-1}$. I want to know whether such type of function has been solved earlier. If answer is in affirmative, I pray, I may kindly be informed of its reference or any other source.
 A: 
I want to know whether such type of function has been solved earlier.

It's been done before, for example it's implemented in SciPy.  If installed, you can access it from Python as scipy.special.lambertw.

If answer is in affirmative, I pray, I may kindly be informed of its reference or any other source.

SciPy is available on GutHub.  The module providing Lambert $W$ is here: https://github.com/scipy/scipy/blob/main/scipy/special/_lambertw.py

Possible issues:
The evaluation can become inaccurate very close to the branch point
at $-1/e$. In some corner cases, lambertw might currently
fail to converge, or can end up on the wrong branch.

The references given are Wikipedia and the PDF already linked before, "On the Lambert W function" from Corless et al.
_lambertw.py is just a top-level Python interface, the very imlpementation is in Cython module _lambertw.pxd.
The implementation basically goes like this:

*

*Determine a starting value for Halley's method. This is accomplished by an asymptotic expansion of $W$:$$W_k(z)\approx \log z + 2\pi ik - \log(\log z + 2\pi ik)$$ There are cases where a Padé approximant is better, and thus used in such cases (branch $k=0$ around $z=0$). Around the branch point $z=-1/e$, yet another expansion is used.


*Apply Halley's method until the desired accuracy  is reached: $$w \mapsto w - 
\dfrac{we^w-z}{e^w(w+1) - \dfrac{(w+2)(we^w-z)}{2w+2}} $$Notice that the accuracy is an input to lambertw.
For details see the code resp. the document linked above and formulae (4.20), (4.22) and (5.9).
A: To stay with Lambert function.
The equations to be solved are
$$-a-e^x y \sin (y)+e^x x \cos (y)=0 \tag 1$$
$$-b+e^x x \sin (y)+e^x y \cos (y)=0 \tag 2$$
Using $(1)$
$$x=W\left(a e^{-y \tan (y)} \sec (y)\right)+y \tan (y)\tag 3$$
Plug in $(2)$ and you need to solve numerically for $y$
$$\frac{a y \sec ^2(y)}{W\left(a e^{-y \tan (y)} \sec (y)\right)}+a \tan (y)-b=0 \tag 4$$
Trying with $a=\pi$ and $b=e$, Newton method gives
$$y=0.394401045253167$$ from which
$$x=1.194608982218346$$
Checking, the lhs is
$$3.14159265358979 + 2.71828182845905\,i$$
Edit
Discarding Lambert function  and using
$$x=\frac{a \cos (y)+b \sin (y)}{b \cos (y)-a \sin (y)}\,y$$ we need to find the zero of function
$$f(y)=\frac{y}{b \cos (y)-a \sin (y)} \exp\Bigg[\frac{a \cos (y)+b \sin (y)}{b \cos (y)-a \sin (y)} \,y\Bigg]-1\tag 5$$   which will have a solution between $0$ and $\tan ^{-1}\left(\frac{b}{a}\right)$.
Because of this a priori bkown upper bound, it could be interesting to let $y=2 \tan ^{-1}(t)$ and to look for the zero of
$$G(t)={2 \left(t^2+1\right) \tan ^{-1}(t)}\,\,\exp\Bigg[2\frac{at^2-2bt-a }{bt^2+2at-b }\tan ^{-1}(t)\Bigg]+(b t^2+2 a t-b)\tag 6$$
and the solution will easily be bracketed before starting an iterative process.
$$\min \Bigg[0,\tan \left(\frac{1}{2} \tan
   ^{-1}\left(\frac{b}{a}\right)\right)\Bigg] < t_* <\max \Bigg[0,\tan \left(\frac{1}{2} \tan
   ^{-1}\left(\frac{b}{a}\right)\right)\Bigg] $$
A: $$(x+iy)e^x(\cos(y)+i\sin(y))=a+ib$$
$$\Rightarrow e^x(x\cos(y)-y\sin(y))=a ~~~\text{and} ~~~e^x(y\cos(y)+x\sin(y))=b$$
$$\Rightarrow x=\frac{(a+b)y\cos(y)}{b(\cos(y)-\sin(y))}$$
Plug $x$ into any of above two equations and numerically solve for $y$.
