The plane through the ray and the origin intersects your ellipsoid to form a planar ellipse. To find this planar ellipse one must establish a basis vector on the plane such that any point on the plane can be represented by two parameters ($\alpha$ and $\beta$). $$ \vec{r} = \alpha \vec{p} + \beta \vec{d} $$
The ellipse on the plane is now given by the equation
$$ \alpha^2 C_{11} + 2 \alpha\beta C_{12} + \beta^2 C_{22} = 1 $$
$$ \pmatrix{\alpha \\ \beta \\ 1}^\top \pmatrix{C_{11} & C_{12} & \\ C_{12} & C_{22} & \\ & & -1} \pmatrix{\alpha \\ \beta \\ 1} =0 $$
where $$\begin{align} C_{11} &= \frac{x_p^2}{a^2} + \frac{y_p^2}{a^2} + \frac{z_p^2}{b^2} \\ C_{22} &= \frac{x_d^2}{a^2} + \frac{y_d^2}{a^2} + \frac{z_d^2}{b^2} \\ C_{12} &= \frac{x_p x_d}{a^2} + \frac{y_p y_d}{a^2} + \frac{z_p z_d}{b^2} \end{align} $$
Now for the fun part. To be perpendicular to the ray the point $\vec{r}$ must obey the projection $$ \vec{d}^\top (\vec{r}-\vec{p}) =0 $$
From the expansion of $\vec{r}$ in terms of the planar coordinates $\alpha$ and $\beta$ this projection is used to solve for
$$ \alpha = 1- \frac{\vec{d}^\top \vec{d}}{\vec{d}^\top \vec{p}} \beta = 1-\frac{x_d^2+y_d^2+z_d^2}{x_d x_p + y_d y_p + z_d z_p} \beta = 1- \lambda \, \beta $$
So given any distance along the ray $\beta$ the 3D position of the point is $\vec{r} = (1-\lambda \,\beta) \vec{p} + \beta \vec{d} $ where $\lambda$ is the fixed ratio $\lambda = \frac{x_d^2+y_d^2+z_d^2}{x_d x_p + y_d y_p + z_d z_p}$.
The points on the ellipse that are perpendicular to the ray solve the equation
$$ (1-\lambda\,\beta)^2 C_{11} + 2 (1-\lambda\,\beta)\beta C_{12} + \beta^2 C_{22} = 1 \\ (C_{11} \lambda^2 -2 C_{12} \lambda + C_{22}) \beta^2 + 2 (C_{12}-C_{11} \lambda) \beta + (C_{11-1}) = 0 $$
This has two solutions (as expected)
$$ \beta = \frac{C_{11} \lambda-C_{12} \pm \sqrt{ C_{11} (\lambda^2-C_{22}) + C_{12}^2 -2 C_{12} \lambda + C_{22}}}{C_11 \lambda^2 -2 C_{12} \lambda + C_{22}} $$
And the back substitution
$$ \alpha = 1 - \lambda\, \beta \\ \vec{r} = \alpha \vec{p} + \beta \vec{d} $$
Edit 1
To find the lines tangent to the planar ellipse, and passing through P you do the following:
- The planar coordinates of P are $(\alpha=1,\beta=0)$, or in homogeneous coordinates $$P=\pmatrix{1,& 0,& 1}$$
- The homogeneous coordinates of the planar ellipse are $$C=\left| \matrix{C_{11} & C_{12} & \\ C_{12} & C_{22} & \\ & & -1} \right|$$
- The polar line of P is $L=C\,P$ $$L = \left[ \matrix{C_{11}, & C_{12},& -1} \right]$$
The polar line intersects the ellipse at the tangent points $Q_1$ and $Q_2$. To find these points form the equation of the line $C_{11} \alpha + C_{12} \beta - 1 =0$ and plug it into the ellipse equation to find $$Q_1 = \pmatrix{ 1 + \tfrac{C_{12} \sqrt{C_{11}-1}}{K}, & - \tfrac{C_{11} \sqrt{ C_{11}-1 }}{K}, & C_{11} }$$ $$Q_2 = \pmatrix{ 1 - \tfrac{C_{12} \sqrt{C_{11}-1}}{K}, & \tfrac{C_{11} \sqrt{ C_{11}-1 }}{K}, & C_{11} }$$
where $K=\sqrt{C_{11} C_{22}-C_{12}^2}$.
The tangent line coordinates are $T_1 = C \,Q_1$ and $T_2 = C \, Q_2$
$$ T_1 = \left[ \matrix{ C_{11}, & C_{12} -K \sqrt{C_{11}-1} , & -C_{11} } \right] $$
$$ T_2 = \left[ \matrix{ C_{11}, & C_{12} +K \sqrt{C_{11}-1} , & -C_{11} } \right] $$
In equation form the tangent lines are
$$ (C_{11}) \alpha + (C_{12} - K \sqrt{C_{11}-1}) \beta - C_{11} = 0 $$
$$ (C_{11}) \alpha + (C_{12} + K \sqrt{C_{11}-1}) \beta - C_{11} = 0 $$
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.
Best Answer
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>$$