If I use symmetry take 2 x integral from 0 to Pi and define the constraints on the parameters, Mathematica gives me an analytical answer to the integral:
In[34]:= Integrate[2 Cos[p b]/ (1 + e Cos[b]), {b, 0, Pi},
Assumptions -> e > 0 && p > 0 && p \[Element] Integers && e <= 1]
Out[34]= -((
2 \[Pi] HypergeometricPFQRegularized[{1/2, 1, 1}, {1 - p, 1 + p}, (
2 e)/(-1 + e)])/(-1 + e))
