Well, I know Kepler's first law:
$$\text{r}=\frac{\text{p}}{1+\epsilon\cos\left(\theta\right)}\tag1$$
Where $\text{p}$ is the semi-latus rectum, $\epsilon$ is the eccentricity of the ellipse, $\text{r}$ is the distance from the Sun to the planet, and $\theta$ is the angle to the planet's current position from its closest approach, as seen from the Sun.
Now, in order to find $\theta$ I need to solve:
$$\frac{\text{n}}{\text{A}}\int_0^{2\pi}\frac{1}{\left(1+\epsilon\cos\left(x\right)\right)^2}\space\text{d}x=\int_0^\theta\frac{1}{\left(1+\epsilon\cos\left(x\right)\right)^2}\space\text{d}x\tag2$$
Where $\epsilon$ is the eccentricity of the ellipse, $\text{A}$ is the number of days in a year (roundabout $365.25636$) and $\text{n}$ is a day in a year (so $1\le\text{n}\le\text{A}$).
Now, I found the integrals and I came out with this expression:
$$\color{red}{\frac{\text{n}}{\text{A}}\cdot\frac{2\pi}{\sqrt{\left(1-\epsilon^2\right)^3}}=\frac{\epsilon}{\epsilon^2-1}\cdot\frac{\sin\left(\theta\right)}{1+\epsilon\cos\left(\theta\right)}-\frac{2}{\sqrt{\left(1-\epsilon^2\right)^3}}\cdot\arctan\left(\frac{\epsilon-1}{\sqrt{1-\epsilon^2}}\cdot\tan\left(\frac{\theta}{2}\right)\right)}\tag3$$
Now I can add in the upper and lower bound because I know the upper and lower bound of $\text{n}$:
$$0.0172093\approx\frac{1}{\text{A}}\cdot\frac{2\pi}{\sqrt{\left(1-\epsilon^2\right)^3}}\le\int_0^\theta\frac{1}{\left(1+\epsilon\cos\left(x\right)\right)^2}\space\text{d}x\le\frac{2\pi}{\sqrt{\left(1-\epsilon^2\right)^3}}\approx6.28582\tag4$$
But now we've got a problem because when I plot the RHS of equation $\left(2\right)$ I got a graph (see the picture below) that never reaches a value of $6.28582$ so where am I wrong? For the values of $\epsilon$ and $\text{p}$ I've: $\epsilon\approx0.0167086$ and $\text{p}\approx1.49558\times10^{11}$
First, a simplification (using $0<\epsilon<1$): $$\arctan\left(\frac{\epsilon-1}{\sqrt{1-\epsilon^2}}\tan\frac\theta2\right) =\arctan\left(-\frac{\sqrt{1-\epsilon}\sqrt{1-\epsilon}}{\sqrt{1-\epsilon}\sqrt{1+\epsilon}}\tan\frac\theta2\right)= -\arctan\left(\sqrt{\frac{{1-\epsilon}}{{1+\epsilon}}}\tan\frac\theta2\right)$$ But the main problem is: As $\theta$ grows from $0$ to $2\pi$, the angle $\frac\theta2$ grows from $0$ to $\pi$, hence crosses $\frac\pi2$ (when $\theta=\pi$). At that point, $\tan\frac\theta2$ goes to $+\infty$ and jumps to $-\infty$ (loosely speaking; more precisely, the expresion is undefined at $\theta=\pi$), and consequently, the $\arctan$ jumps from $\frac\pi2$ to $-\frac\pi2$.
The easiest way to get out of this is to replace $\arctan x$ with a function we might call $$ \widetilde{\arctan}\,x:=\begin{cases}\frac\pi2&\text{if }x=\pm\infty,\\\arctan x&\text{if }x\ge 0,\\\pi+\arctan x&\text{if }x<0.\end{cases}$$ This gets rid of both jump incontinuities at once. To understand this more deeply, recall that $\arctan x$ is just one of all possible angles with tangens equal to $x$. More precisely, the convention is that we pick the angle that is between $-\frac\pi2$ and $\frac\pi 2$. But that may not be a good choice in all circumstances.