Solution 1:

Preliminaries

$\newcommand{\Res}{\operatorname*{Res}}$ $\Res\limits_{w=0}\left(w^k\,e^{\lambda/4(w+1/w)}\right)$ is a modified Bessel function of the first kind: $$ \begin{align} \Res_{w=0}\left(w^k\,e^{\lambda/4(w+1/w)}\right) &=\left[w^{-k-1}\right]\sum_{n=0}^\infty\frac{(\lambda/4)^n}{n!}\sum_{j=0}^n\binom{n}{j}w^{-n+2j}\tag{1a}\\ &=\sum_{j=0}^\infty\sum_{n=j}^\infty\frac{(\lambda/4)^n}{n!}\binom{n}{j}[n=2j+k+1]\tag{1b}\\ &=\sum_{j=0}^\infty\frac{(\lambda/4)^{2j+k+1}}{(2j+k+1)!}\binom{2j+k+1}{j}\tag{1c}\\ &=\sum_{n=0}^\infty\frac{(\lambda/4)^{2n+k+1}}{n!(n+k+1)!}\tag{1d}\\[9pt] &=I_{k+1}(\lambda/2)\tag{1e} \end{align} $$ Explanation:
$\text{(1a)}$: residue is the coefficient of $1/w$; apply the Binomial Theorem
$\text{(1b)}$: change order of summation and apply coefficient-of operator
$\text{(1c)}$: apply Iverson brackets
$\text{(1d)}$: simplify
$\text{(1e)}$: compare series


Evaluation of the Integral

For $\boldsymbol{|z|\gt1}$

Using $(1)$, we can evaluate the integral as a sum of modified Bessel functions of the first kind: $$ \begin{align} S(z) &=\frac1{2\pi}\int_0^{2\pi}\frac{e^{i\phi}+z}{e^{i\phi}-z}e^{-\lambda\sin^2(\phi/2)}\,\mathrm{d}\phi\tag{2a}\\ &=\frac1{2\pi i}\oint_{\partial B_1}\frac{w+z}{w-z}e^{-\lambda/4(2-w-1/w)}\frac{\mathrm{d}w}w\tag{2b}\\ &=\frac{e^{-\lambda/2}}{2\pi i}\oint_{\partial B_1}\left(\frac2{w-z}-\frac1w\right)e^{\lambda/4(w+1/w)}\,\mathrm{d}w\tag{2c}\\[3pt] &=-e^{-\lambda/2}\Res_{w=0}\left(\left(\frac2z\frac1{1-w/z}+\frac1w\right)e^{\lambda/4(w+1/w)}\right)\tag{2d}\\ &=-e^{-\lambda/2}\left(I_0(\lambda/2)+\sum_{k=0}^\infty\frac2{z^{k+1}}I_{k+1}(\lambda/2)\right)\tag{2e} \end{align} $$ Explanation:
$\text{(2a)}$: definition
$\text{(2b)}$: $w=e^{i\phi}$
$\text{(2c)}$: partial fractions
$\text{(2d)}$: since $|z|\gt1$, ignore the residue at $w=z$
$\text{(2e)}$: apply $(1)$

For $\boldsymbol{|z|\lt1}$

If we substitute $w\mapsto1/w$ in $(2b)$ we see that $$ S(z)=-S(1/z)\tag3 $$ Therefore, $$ S(z)=e^{-\lambda/2}\left(I_0(\lambda/2)+\sum_{k=0}^\infty2z^{k+1}I_{k+1}(\lambda/2)\right)\tag4 $$


The Real Part near $\boldsymbol{|z|=1}$

Since all the coefficients in the Laurent expansion are real, we have that $$ S\left(\bar z\right)=\overline{S(z)}\tag5 $$ Combining $(3)$ and $(5)$ gives $$ \overline{S(z)}=-S(1/\bar z)\tag6 $$ If $z$ is just inside the unit circle, then $1/\bar z$ is just on the other side of the circle and $$\newcommand{\Re}{\operatorname{Re}} \begin{align} 2\Re(S(z)) &=S(z)+\overline{S(z)}\tag{7a}\\ &=S(z)-S(1/\bar z)\tag{7b}\\ &=2e^{-\lambda/2}e^{\lambda/4(z+1/z)}\tag{7c}\\ \Re(S(z)) &=e^{\lambda/2\Re(z-1)}\tag{7d}\\ \Re(S(1/\bar z)) &=-e^{\lambda/2\Re(z-1)}\tag{7e} \end{align} $$ Explanation:
$\text{(7a)}$: formula for the real part of a complex number
$\text{(7b)}$: apply $(6)$
$\text{(7c)}$: $\text{(2c)}$ says that we include the residue at $w=z$
$\phantom{\text{(7c):}}$ for $|z|\lt1$, but not for $|z|\gt1$
$\text{(7d)}$: for $|z|=1$, $z+1/z=2\Re(z)$
$\text{(7e)}$: apply the real part of $(6)$


Graph Along the Real Axis

Along the real axis, for $\lambda=1$:

graph along the real axis

The function is discontinuous at $|x|=1$, which is to be expected since $(7)$ says $$ S(1^{\pm})=\mp1\qquad\text{and}\qquad S(-1^{\pm})=\pm e^{-\lambda}\tag8 $$ However, $S(x)$ is continuous at $x=0$, and $$ S(0)=e^{-\lambda/2}I_0(\lambda/2)\tag9 $$