note: I believe the Earth coordinates should be written with capital letters $(X, Y, Z)$ but quotes from the other answer use lower case $(x, y, z)$. Considering the radical term barely fits with little letters also, I've kept lower case in most places. For this answer upper and lower cases represent the same coordinates.
tl;dr: Plug in your satellite's position ($x, y, z$) $$\mathbf{x_0} = x\mathbf{\hat{x}} + y\mathbf{\hat{y}} + z\mathbf{\hat{z}}$$ and the vector's normal ($u, v, w$) $$\mathbf{\hat{n}} = u\mathbf{\hat{x}} + v\mathbf{\hat{y}} + w\mathbf{\hat{z}}$$ to get $t$ the length of the vector from satellite to the ellipsoid, then use $$\mathbf{x_{intersect}} = \mathbf{x_0} + t\mathbf{\hat{n}}$$ to get the coordinates of the intersection with the reference ellipsoid.
Then use either the analytical solution in the second section below, or one of the iterative solutions the references there to convert those to latitude and longitude.
Intersection math to get (x, y, z)
From this handy answer in GIS SE:
The WGS84 reference ellipsoid is a biaxial (and oblate) ellipsoid. It's shorter at the poles than the equator, and the equator is a circle.
The equation for it is:
$$\frac{x^2}{a^2} + \frac{y^2}{a^2} + \frac{z^2}{b^2} - 1 = 0$$
where $x, y, z$ is a point on the ellipsoid's surface, and
where $a$ is the semi-major axis (6,378,137 meters) and $b$ is the semi-minor axis of the WGS84 ellipsoid (6,356,752.3142 meters).
I'll continue to use the math from this answer and reformat in MathJax and later Python:
Your satellite's position is $x,y, z$ and the direction normal is $u, v, w$. The length of the vector from the satellite to the first and closest intersection (assuming of course that your satellite is outside of the Earth:
t = -(1/(b^2 (u^2 + v^2) + a^2 w^2)) * (b^2 (u x + v y) + a^2 w z + 1/2 Sqrt[
4 (b^2 (u x + v y) + a^2 w z)^2 -
4 (b^2 (u^2 + v^2) + a^2 w^2) (b^2 (-a^2 + x^2 + y^2) + a^2 z^2)])
In MathJax:
$$A = -\frac{1}{b^2(u^2+v^2) + a^2w^2} $$
$$B = b^2(ux + vy) + a^2wz $$
$$C = \frac{1}{2}\sqrt{4(b^2(ux + vy) + a^2wz)^2 - 4(b^2(u^2 + v^2) + a^2*w^2) (b^2(-a^2 + x^2 + y^2) + a^2z^2)} $$
$$ t = A(B+C)$$
which is easier to read than the way MathJax displays it all in one line:
$$t = -\frac{1}{b^2(u^2+v^2) + a^2w^2} \left( b^2(ux + vy) + a^2wz + \frac{1}{2}\sqrt{4(b^2(ux + vy) + a^2wz)^2 - 4(b^2(u^2 + v^2) + a^2w^2) (b^2(-a^2 + x^2 + y^2) + a^2z^2)} \right)$$
Implemented in Python:
import numpy as np
a, b = 0.9, 1.1
asq, bsq = a**2, b**2
x, y, z = 5.0, 0.0, 0.0
u, v, w = -np.sqrt(1 - 0.1**2 - 0.1**2), 0.1, 0.1
xsq, ysq, zsq = x**2, y**2, z**2
usq, vsq, wsq = u**2, v**2, w**2
A = -(1/(bsq*(usq + vsq) + asq*wsq))
B = bsq*(u*x + v*y) + asq*w*z
C = 0.5*np.sqrt(4*(bsq*(u*x + v*y) + asq*w*z)**2 -
4*(bsq*(usq + vsq) + asq*wsq) *
(bsq*(-asq + xsq + ysq) + asq*zsq))
t = A * (B + C)
print "t: ", t
xyz, uvw = np.array([x, y, z]), np.array([u, v, w])
xyzi = xyz + t*uvw
xi, yi, zi = xyzi
print "point of intersection: ", xyzi
print "check, is it zero within roundoff? ", xi**2/asq + yi**2/asq + zi**2/bsq - 1
it yields:
t: 4.33961998769
point of intersection: [ 0.70399539 0.433962 0.433962 ]
check, is it zero? -8.54871728961e-15
The solution does indeed fall on the ellipsoid.
Convert (x, y, z) to latitude and longitude
We need to invert the equations (shown in this answer) in order to get the latitude and longitude on the WSG84 ellipsoid corresponding to these coordinates.
According to Wikipedia's [Geographic_coordinate_conversion#From_geodetic_to_ECEF_coordinates][1]
The 3D cartesian coordinates $X, Y, Z$ in Earth-centered, Earth-fixed coordinates assuming an ellipsoidal shape is given by:
$$X = \left(N(\phi) + h \right) \cos\phi \cos\lambda $$
$$Y = \left(N(\phi) + h \right) \cos\phi \sin\lambda $$
$$Z = \left(\frac{b^2}{a^2} N(\phi) + h \right) \sin\phi $$
where $\phi, \lambda, h$ are latitude, longitude, and altitude, and $a, b$ are the equatorial and polar radii of the ellipsoid used, and
$$N(\phi) = \frac{a^2}{\sqrt{a^2\cos^2\phi + b^2 \sin^2\phi}}. $$
The inversion is addressed in Datum Transformations of GPS Positions; Application Note found in this GIS question as well as in this link in comments.
From the Application Note I'll show the analytical solution.
Definitions:
$$a = 6378137 $$ $$b = a(1-f) = 6356752.31424518$$ $$f = \frac{1}{298.257223563}$$ $$e = \sqrt{\frac{a^2-b^2}{a^2}}$$ $$e' = \sqrt{\frac{a^2-b^2}{b^2}}.$$
Solving for longitude is trivial. The Application Note gives:
$$\lambda = \arctan\left(\frac{y}{x}\right),$$
however in practice you can't blindly use arctan because it can't distinguish all four quadrants. So instead use
$$\lambda = \arctan2(y, x).$$
Solving for latitude analytically, the Application Note gives
$$\phi= \arctan\left( \frac{z + e'^2b \sin^3(\theta)}{p-e^2a\cos^3(\theta)} \right)$$
where
$$p = \sqrt{x^2+y^2}$$ $$\theta = \arctan\left(\frac{za}{pb}\right).$$
Again you should probably use $arctan2()$ instead of minding the quadrants.
You can also find iterative solutions in this reference as well as the other (mentioned already above). Iterative solutions were probably significantly faster in the past (and might be now as well) which would be important if you were doing ray-tracing for millions of points, e.g. generating an image of the Earth's surface from the satellites point of view with an accurate representation.