Point d'intersection pour un rayon et la surface de la Terre

11

Disons que j'ai un vecteur rayon provenant de quelque part dans l'espace. Un exemple pourrait être un rayon de lumière du soleil. Comment puis-je calculer le point d'intersection (s'il existe) entre le rayon et la surface de la Terre? J'utilise des coordonnées cartésiennes (ECEF) et j'aimerais que la géométrie ellipsoïde de la Terre soit prise en compte dans le calcul.

Pris
la source

Réponses:

12

C'est simple mais désordonné.

Parce que vous travaillez dans ECEF, vous avez sans doute aussi l'origine du rayon (x, y, z) et le vecteur de direction (u, v, w) en coordonnées ECEF. Supposons pour le moment que pendant le temps de voyage à la surface de la Terre, la Terre ne bouge pas sensiblement. (La partie la plus rapide de la terre en rotation, l'équateur, se déplace d'environ 0,45 km / s et la lumière se déplace d'environ 300 000 km / s, donc un rayon provenant, disons, de 1000 km au-dessus de la terre et dirigé plus ou moins directement vers l'équateur prendra 1/300 seconde pour l'atteindre, pendant laquelle l'équateur aura déplacé 1,5 mètre: c'est probablement une erreur acceptable.)

Il suffit de calculer l'intersection de la ligne paramétrée

t --> (x,y,z) + t*(u,v,w)

avec la surface de la terre, qui peut être considérée comme le zéro de la fonction

(x/a)^2 + (y/a)^2 + (z/b)^2 - 1

a est le demi-grand axe (6 378 137 mètres) et b est le demi-petit axe de l' ellipsoïde WGS84 (6 356 752 3142 mètres). Branchez la première formule dans la seconde et résolvez pour t en termes de x, y, z, u, v, w . C'est une équation quadratique, donc vous obtenez jusqu'à deux solutions: une pour entrer dans la terre et une autre pour la quitter à nouveau (ce qui arriverait, par exemple, pour un neutrino). Choisissez la solution pour laquelle la distance est la plus courte. Cela donne

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)])

Branchez cette valeur dans la première équation pour obtenir le point d'intersection.

Pour un rayon originaire de loin, mais pas terriblement loin ( par exemple, du soleil mais pas de l'extérieur du système solaire), commencez par une estimation grossière du temps T qu'il faudrait pour atteindre la terre (en secondes): vous pourriez utilisez la distance de (x, y, z) au centre de la terre, par exemple. Modifiez les coordonnées de départ (x, y, z) pour tenir compte de la rotation de la Terre pendant cette période: cela changera les coordonnées de départ en

(x*c + y*s, -x*s + y*c, z)

(le point semblera reculer ) où c et s sont le sinus et le cosinus de 0,000072921150 * T radians . Calculez l'intersection d'un rayon à partir de cet emplacement mis à jour. Vous pouvez être éloigné de près de 10 mètres en raison de l'utilisation d'un temps estimé. Si cette question est importante, réestimer le temps écoulé en fonction de ce point d'intersection et répéter le calcul avec la nouvelle valeur de T .

whuber
la source
Sensationnel. Merci beaucoup pour la réponse incroyablement détaillée!
Pris
Cette question semble trop proche de celle-ci: gis.stackexchange.com/questions/86031/… Cependant je n'utilise pas ECEF. Peut être résolu de manière similaire @whuber? Thx
alvarolb
1
@alvarolb Une référence montrant comment convertir entre ECEF et (lon, lat, altitude) est donnée sur gis.stackexchange.com/questions/20714 .
whuber
Ne peut-on pas calculer directement sans passer de la géodésique à l'EFEC et les renvoyer à la géodésique?
alvarolb