An integral arising from Kepler's problem $\frac{1}{2\pi}\int_0^{2\pi}\ \frac{\sin^2\theta\ \text{d}\theta}{(1 + \epsilon\cos\theta)^2}$ I'm dealing with this integral in my spare time, since days and days, and it's really interesting. I'll provide to write what I tried until now, and I would really appreciate some help in understanding how to continue.
$$\Phi(\epsilon) = \frac{1}{2\pi}\int_0^{2\pi}\ \frac{\sin^2\theta\ \text{d}\theta}{(1 + \epsilon\cos\theta)^2}$$
Note: The integral arises in studying the Kepler's problem, which is why $\epsilon\in [0;\ 1]$.
I started with the passages to the complex plane:
$$\sin\theta = \frac{1}{2i}(z - z^{-1}) ~~~~~~~ \cos\theta = \frac{1}{2}(z + z^{-1}) ~~~~~~~ \text{d}\theta = \frac{\text{d}z}{iz}$$
thence
$$\begin{align*}
\Phi(\epsilon) & = \frac{1}{2\pi}\oint_{|z| = 1}\ \frac{-\frac{1}{4}(z - z^{-1})^2\ \text{d}z}{iz\left[1 + \epsilon\frac{(z + z^{-1}}{2}\right]^2}
\\\\
& = -\frac{1}{8\pi i}\oint_{|z| = 1}\ \frac{(z - z^{-1})^2\ \text{d}z}{z\left[\frac{2 + \epsilon(z + z^{-1})}{2}\right]^2}
\\\\
& = -\frac{1}{4\pi i}\oint_{|z| = 1}\ \frac{(z - z^{-1})^2\ \text{d}z}{z(2 + \epsilon(z + z^{-1})^2}
\\\\
& = -\frac{1}{4\pi i}\oint_{|z| = 1}\ \frac{\left(\frac{z^2 - 1}{z}\right)^2\ \text{d}z}{z\left(2 + \epsilon\frac{z^2+1}{z}\right)^2}
\\\\
& = -\frac{1}{4\pi i}\oint_{|z| = 1}\ \frac{(z^2-1)^2\ \text{d}z}{z(2z + \epsilon(z^2+1))^2}
\end{align*}$$ 
Poles et cetera
Poles are at $z_0 = 0$ and when $2z + \epsilon(z^2+1) = 0$, namely 
$$z_1 = \frac{1+\sqrt{1 - \epsilon^2}}{\epsilon} ~~~~~~~~~~~ z_2 = \frac{1-\sqrt{1 - \epsilon^2}}{\epsilon}$$
Checking that $z_1\cdot z_2 = 1$ is ok, the question is now:
Question 1.: which one, aside $z_0$ lies into the unitary circle?
Question 2.: How to proceed? 

I am assuming $0 < \epsilon < 1$ in the following. The cases
$\epsilon = 0$ and $\epsilon = 1$ have to be handled separately,
see below.
There is a small error in your calculation, the root of
$2z + \epsilon(z^2+1) = 0$ are 
$$
z_1 = \frac{-1-\sqrt{1 - \epsilon^2}}{\epsilon} \, , \quad z_2 = \frac{-1+\sqrt{1 - \epsilon^2}}{\epsilon}
$$

Question 1.: which one, aside $z_0$ lies into the unitary circle?

From
$$
 1+\sqrt{1 - \epsilon^2} > 1 > \epsilon
$$
it follows that $z_1 < -1$ and consequently, $-1 < z_2 < 0$, i.e. $z_2$
is inside the unit disk and $z_1$ outside.

Question 2.: How to proceed?

You already have that
$$
\Phi(\epsilon) = -\frac{1}{4\pi i}\oint_{|z| = 1} \frac{(z^2-1)^2\ }{z(2z + \epsilon(z^2+1))^2} \, dz =
\frac{1}{2\pi i} \oint_{|z| = 1} f(z) \, dz
$$
with
$$
 f(z) := -\frac{1}{2 \epsilon^2} \frac{(z^2-1)^2\ }{z(z-z_1)^2(z-z_2)^2} 
$$
From the residue theorem it follows that
$$
\Phi(\epsilon) = \text{Res}(f, 0) + \text{Res}(f, z_2) \, .
$$
$f$ has a simple pole at $z = 0$, therefore
$$
  \text{Res}(f, 0) = \lim_{z \to 0} z f(z) = -\frac{1}{2 \epsilon^2} \, .
$$
$f$ has a double pole at $z = z_2$.
One possible method to compute the residue is the
limit formula for higher order poles:
$$
 \text{Res}(f, z_2) = \lim_{z \to z_2} \frac{d}{dz} \bigl((z-z_2)^2 f(z) \bigr) \, .
$$

Special cases: For $\epsilon  = 0$ the integral becomes
$$
 \Phi(0) = -\frac{1}{4\pi i}\oint_{|z| = 1} \frac{(z^2-1)^2}{2z^3} \, dz
= -\frac{1}{8\pi i}\oint_{|z| = 1} \bigl( z - \frac 2z + \frac{1}{z^3}\bigr) \, dz
$$
which can easily be computed using the residue theorem. 
Alternatively,
$$
  \Phi(0) = \frac{1}{2\pi}\int_0^{2\pi} \sin^2\theta \, d\theta
$$
can be computed directly.
For $\epsilon = 1$, the integral becomes
$$
\Phi(1) = -\frac{1}{4\pi i}\oint_{|z| = 1} \frac{(z-1)^2}{z(z+1)^2} \, dz
$$
which is infinite due to the singularity at $z=-1$. Alternatively,
$$
 \Phi(1) = \frac{1}{2\pi}\int_0^{2\pi} \frac{\sin^2\theta}{(1 + \cos\theta)^2} \, d\theta 
 = \frac{1}{2\pi}\int_0^{2\pi} \frac{(2 \sin \frac{\theta}{2} \cos \frac{\theta}{2})^2 }{(2 \cos^2 \frac{\theta}{2})^2} \, d\theta 
 = \frac{1}{2\pi}\int_0^{2\pi} \tan^2 \frac{\theta}{2}\, d\theta
$$
is infinite.
