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
Under construction...
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 as well as the iterative 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(\frac{y}{x}),$$
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= $$
under construction, stay tuned...