Number of normals from a point to an ellipse

From a typical point $P$ inside an ellipse, how many points $Q_i$ on the ellipse have $PQ_i$ normal to the ellipse? Someone asked me at school many years ago but I don't think I worked it out.

For every ellipse $\mathcal{E}$, there is a curve called ellipse evolute $^\color{blue}{[1]}$ associated with it. The ellipse evolute is the locus of centers of curvature $^\color{blue}{[2]}$ for $\mathcal{E}$. It is also the envelope of the normal lines of \mathcal{E}$$^\color{blue}{[2]}. For a point P inside \mathcal{E}, the number of points Q on \mathcal{E} where PQ normal to \mathcal{E} can be either 2, 3 or 4. It depends on whether P is lying outside, on or inside the ellipse evolute. The graph below illustrates what happens when P lies on the ellipse evolute (the black star shaped curve) and the three Q such that PQ normal to \mathcal{E}. Following is a brief analysis of the problem. For an alternate and more complete treatment, take a look at the paper Apollonius' ellipse and evolute revisited by Frederick Hartmann and Robert Jantzen. Let use choose a coordinate system such that \mathcal{E} is given by the equation$$\frac{x^2}{a^2} + \frac{y^2}{b^2} = 1$$where a > b > 0 are the semi-major and semi-minor axis of \mathcal{E}. Let c = \sqrt{a^2-b^2}, the two foci of \mathcal{E} are located at (\pm c,0). Let P = (p,q) be a point inside \mathcal{E} and Q = (x,y) be a point on \mathcal{E}. It is known that the normal of \mathcal{E} at Q is along the direction (\frac{x}{a^2},\frac{y}{b^2}). The condition for PQ normal to \mathcal{E} is given by$$\frac{x}{a^2} : \frac{y}{b^2} = x - p : y - q \quad\iff\quad \frac{x}{a^2}(y-q) - \frac{y}{b^2}(x-p) = 0$$This is the equation for a hyperbola \mathcal{H} ( the orange hyperbola in above graph ) with P lying on it. Since P is inside \mathcal{E} and the two arms of the branch of \mathcal{H} holding P extends to infinity. Each of the arm will intersect \mathcal{E} at least once. This means \mathcal{E} and \mathcal{H} intersected at least twice. Since five points determine a conic, the number of intersections between \mathcal{E} and \mathcal{H} is at most 4. To determine the actual number of intersections, let us introduce a new coordinate system$$(u,v) = (\frac{x}{a},\frac{y}{b})\quad\iff\quad (x,y) = (au, bv)$$In terms of (u,v), the ellipse becomes the unit circle u^2 + v^2 = 1 and the hyperbola \mathcal{H} becomes$$uv + \tilde{q} u - \tilde{p} v = 0\quad\text{ where }\quad \begin{cases}\tilde{p} = \frac{ap}{c^2}\\ \tilde{q} = \frac{bq}{c^2}\end{cases}$$Consider the rational parametrization of the circle$$t\quad \mapsto\quad (\frac{1-t^2}{1+t^2}, \frac{2t}{1+t^2}) $$The condition that a point on the circle intersect \mathcal{H} becomes$$\begin{align} & \left(\frac{1-t^2}{1+t^2}\right)\left(\frac{2t}{1+t^2}\right) + \tilde{q}\left(\frac{1-t^2}{1+t^2}\right) - \tilde{p}\left(\frac{1-t^2}{1+t^2}\right) = 0\\ \\ \iff & \tilde{q}t^4 + 2 (\tilde{p}+1) t^3 + 2(\tilde{p}-1) t - \tilde{q} = 0 \end{align}$$The last expression is a quartic equation which "usually" have either two or four real roots. When we changes the position of P, the number of roots jumps at those place where the discriminant of the quartic polynomial vanishes. By brute force computation, the discriminant of the quartic polynomial is given by$$\Delta(\tilde{p},\tilde{q}) = -256\left[ (\tilde{p}^2+\tilde{q}^2-1)^3 + 27\tilde{p}^2\tilde{q}^2\right]$$We have 3 possible cases: 1. When P is near O, \Delta(\tilde{p},\tilde{q}) > 0. The quartic equation has either 4 or 0 real roots. Since we know \mathcal{H} intersect \mathcal{E} at least twice, there are four Q on \mathcal{E} such that PQ is normal to \mathcal{E}. 2. When P is far away from O, \Delta(\tilde{p},\tilde{q}) < 0. The quartic equation has only 2 real roots and hence there are only two Q that make PQ normal to \mathcal{E}. 3. On the special case \Delta(\tilde{p},\tilde{q}) = 0, the quartic equation has 3 distinct real roots. One of them is a double root which corresponds to \mathcal{H} is tangent to \mathcal{E} at some points. There are three Q that make PQ normal to \mathcal{E}. As mentioned before, the picture above illustrates the 3^{rd} case. The black star shaped curve is the ellipse evolute where \Delta(\frac{ap}{c^2},\frac{bq}{c^2}) = 0. When P is lying on it, one branch of the hyperbola \mathcal{H} will be touching the ellipse \mathcal{E}. If you move P inside the ellipse evolute, the bottom branch of \mathcal{H} will move inwards too and \mathcal{H} will start to intersect \mathcal{E} at four places. To obtain a simpler expression for the ellipse evolute, let r = |\tilde{p}|^{\frac23}, s = |\tilde{q}|^{\frac23}, we have$$\begin{align} \Delta(\tilde{p},\tilde{q}) = 0 \iff & (r^3 + s^3 - 1)^3 + 27r^3s^3 = 0\\ \iff & r^3 + s^3 - 1 + 3rs = 0\\ \iff & (r+s-1)(s^2-rs+s+r^2+r+1) = 0\\ \iff & (r+s-1)\left((s+r+2)^2+3(s-r)^2\right) = 0\\ \iff & r + s = 1\\ \iff & |ap|^{\frac23} + |bq|^{\frac23} = c^{\frac43} \end{align}$Reproducing the equation for ellipse evolute commonly appear in the literature. Notes$\color{blue}{[1]}$ellipse evolute is a special case of a kind of curve called astroid.$\color{blue}{[2]}\$ The wiki page of evolute has the definition of center of curvature. It also has a nice animation showing the ellipse evolute as an envelop of the normals.