One possible solution can be as follows:
Let $\rm P{(x_1,y_1)} \ Q{(x_2,y_2)}, and\ R{(h,k)}, $ be the points.
Let ellipse ($\rm S$) be a special ellipse with centre (0,0)
$$\begin{align}\rm S &=\rm \dfrac{x^2}{a^2}+ \dfrac{y^2}{b^2} = 1 \\
\rm Equation \ of \ tangent\ at\ R &:\rm \dfrac{xh}{a^2}+ \dfrac{yk}{b^2}-1 = 0 \ \ \ \ \ \ ... (i)\end{align}$$
Since the line segment $\rm PR$ (slope $\rm m_2$) is equally inclined to the tangent ( slope $\rm m_1$) as the segment $\rm QR$ ( slope $\rm m_3$), we can equalise the angles.
$$\rm\tan\theta=\dfrac{m_1 - m_2}{1+m_1m_2}$$
$$\rm\tan\phi=\dfrac{m_1 - m_3}{1+m_1m_3}$$
$$\rm \implies \tan\theta = \tan\phi \ \ \ \ \ \ ... (ii)$$
Find $\rm m_1, m_2, m_3$
You'll get $2$ equations with $2$ unknowns $\rm h$ and $\rm k$.
Note In my answer, I used eqn of standard ellipse with centre at origin. But this method will work for any ellipse. For finding tangent eqn, just put $T=0$
EDIT
Here is a graph which shows tangent, and a line inclined equally to it (it tangent is the angle bisector of the lines.) : Desmos
The $\nabla$ operator is known as the gradient operator. It's a vector of all of the partial derivatives of the function with respect to all of its variables, i.e., $$\nabla f=\bigg<\frac{\partial f}{\partial x},\frac{\partial f}{\partial y},\frac{\partial f}{\partial z}\bigg>$$The gradient vector will give you your desired normal vector.
Now, all of that may sound like gibberish to you, so let's break down the intuition behind it.
Suppose our surface involved 1 variable, instead of 3, such as $f=\frac{x^2}{a^2}$. How would we find the normal vector at a given point? What we would want to do, is find the tangent line, and then find the vector perpendicular to it.
Now what does the tangent line look like? We know that it must have a slope equal to the derivative of the function at the given point! In other words, it would be of the form $$\frac{df}{dx}\cdot x+c=0$$With three variables, we can draw an analogy with partial derivatives instead of full ones (since we have more than two variables) and likewise create an equation of a tangent plane, which would look something like$$\frac{\partial f}{\partial x}\cdot x+\frac{\partial f}{\partial y}\cdot y+\frac{\partial f}{\partial z}\cdot z + C=0$$So, now all we want is the normal vector to this plane. However, we know through the magic of the cross product, that the normal vector to a plane $$ax+by+cz+d=0$$ is $$<a,b,c>$$So, you can see why the gradient gives you the desired result.
So, your normal vector at a given point $(x,y,z)$ would be $$2\cdot\bigg<\frac x{a^2},\frac y{b^2},\frac z{c^2}\bigg>$$
Best Answer
Reflection angle $\theta$ can be found from the sine rule:
$$ {r\over\sin\alpha}={d\over\sin\theta}, $$ where $d$ is the distance from the light source to the centre of the sphere, while $\alpha$ is the angle the incident ray makes with the line joining the source and the centre.
You can obtain the limiting value of $\alpha$ by setting $\sin\theta=1$ in the above formula.