Es simple pero desordenado.
Porque estás trabajando en ECEF, presumiblemente tienes el origen del rayo (x, y, z) y el vector de dirección (u, v, w) en coordenadas ECEF también. Por el momento, asumamos que durante el viaje a la superficie terrestre, la tierra no se mueve apreciablemente. (La parte más rápida de la tierra rotante, el Ecuador, se mueve alrededor de 0.45 km/seg y la luz se mueve aproximadamente 300,000 km/seg, por lo que un rayo que se origine, por ejemplo, a 1000 km sobre la tierra y se dirija más o menos directamente hacia abajo hacia el Ecuador tardará 1/300 de segundo en alcanzarlo, durante el cual el Ecuador se habrá movido 1.5 metros: probablemente ese sea un error aceptable.)
Solo necesitamos calcular la intersección de la línea parametrizada
t --> (x, y, z) + t*(u, v, w)
con la superficie de la tierra, que puede considerarse como el conjunto cero de la función
(x/a)^2 + (y/a)^2 + (z/b)^2 - 1
donde a es el semieje mayor (6,378,137 metros) y b es el semieje menor del elipsoide WGS84 (6,356,752.3142 metros). Sustituye la primera fórmula en la segunda y resuelve para t en términos de x, y, z, u, v, w. Es una ecuación cuadrática, por lo que obtienes hasta dos soluciones: una para entrar a la tierra y otra para salir nuevamente (lo que ocurriría, por ejemplo, para un neutrino). Elige la solución para la cual la distancia es la más corta. Esto da
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)])
Sustituye este valor en la primera ecuación para obtener el punto de intersección.
Para un rayo que se origine lejos, pero no demasiado lejos (por ejemplo, desde el sol pero no desde fuera del sistema solar), comienza con una estimación aproximada del tiempo T que debería tomar llegar a la tierra (en segundos): podrías usar la distancia desde (x, y, z) hasta el centro de la tierra, por ejemplo. Modifica las coordenadas iniciales (x, y, z) para tener en cuenta la cantidad de rotación de la tierra durante este tiempo: esto cambiará las coordenadas iniciales a
(x*c + y*s, -x*s + y*c, z)
(el punto parecerá moverse hacia atrás) donde c y s son el seno y el coseno de 0.000072921150 * T radianes. Calcula la intersección para un rayo que comienza en esta ubicación actualizada. Puedes tener un error de hasta 10 metros más o menos debido al uso de un tiempo estimado. Si esto es importante, vuelve a estimar el tiempo transcurrido basado en este punto de intersección y repite el cálculo con el nuevo valor de T.