Say I have a ray vector originating from somewhere out in space. An example might be a ray of light from the sun. How can I calculate the point of intersection (if it exists) between the ray and the Earth's surface? I'm using Cartesian coordinates (ECEF) and I'd like the Earth's ellipsoid Geometry to factor into the calculation.
It's simple but messy.
Because you're working in ECEF, presumably you have the ray's origin (x,y,z) and direction vector (u,v,w) in ECEF coordinates, too. For the moment let's assume that during the time of travel to the earth's surface, the earth does not appreciably move. (The fastest part of the rotating earth, the Equator, moves about 0.45 km/sec and light moves about 300,000 km/sec, so a ray originating, say, 1000 km above earth and headed more or less straight down towards the Equator will take 1/300 second to reach it, during which the Equator will have moved 1.5 meters: that's probably an acceptable error.)
We only need to compute the intersection of the parameterized line
t --> (x,y,z) + t*(u,v,w)
with the earth's surface, which can be considered the zero set of the function
(x/a)^2 + (y/a)^2 + (z/b)^2 - 1
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). Plug the first formula into the second and solve for t in terms of x, y, z, u, v, w. It's a quadratic equation, so you get up to two solutions: one for entering the earth and another for leaving it again (which would happen, for instance, for a neutrino). Choose the solution for which the distance is shortest. This gives
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)])
Plug this value into the first equation to obtain the point of intersection.
For a ray originating far away, but not terribly far away (e.g., from the sun but not from outside the solar system), start with a crude estimate of the time T it should take to reach the earth (in seconds): you could use the distance from (x,y,z) to the earth's center, for instance. Modify the starting coordinates (x,y,z) to account for the amount of the earth's rotation during this time: this will change the starting coordinates to
(x*c + y*s, -x*s + y*c, z)
(the point will seem to move backwards) where c and s are the sine and cosine of 0.000072921150 * T radians. Compute the intersection for a ray starting at this updated location. You may be off by as much as 10 meters or so due to the use of an estimated time. If this matters, re-estimate the elapsed time based on this point of intersection and repeat the calculation with the new value of T.