Comment calculez-vous le rayon de la Terre à une latitude géodésique donnée?

(Je vois qu'il y a une équation sur wikipedia qui fait exactement ce que je demande mais il n'y a pas de références. Je n'ai aucun moyen de confirmer la validité de cette équation!)

Je comprends déjà la différence entre Latitude géocentrique et Latitude géodésique.

En supposant que les rayons semi-majeurs aet semi-mineurs connus bsont donnés. Comment calculez-vous le rayon à une latitude géodésique donnée?

J'ai besoin d'une sorte de confirmation d'expert (dérivation, lien vers la dérivation, confirmation de l'expert, explication, etc.).

Cette question suppose un modèle ellipsoïdal de la terre. Sa surface de référence est obtenue en faisant tourner une ellipse autour de son petit axe (tracée verticalement par convention). Une telle ellipse n'est qu'un cercle qui a été étiré horizontalement d'un facteur a et verticalement d'un facteur b . En utilisant le paramétrage standard du cercle unitaire,

t --> (cos(t), sin(t))

(qui définit le cosinus et le sinus), on obtient un paramétrage

t --> (a cos(t), b sin(t)).

(Les deux composantes de cette paramétrisation décrivent un parcours autour de la courbe: elles précisent, en coordonnées cartésiennes, notre position au "temps" t .)

La latitude géodésique , f , de tout point est l'angle que "up" fait avec le plan équatorial. Lorsque a diffère de b , la valeur de f diffère de celle de t (sauf le long de l'équateur et aux pôles).

Dans cette image, la courbe bleue est un quadrant d'une telle ellipse (considérablement exagérée par rapport à l'excentricité de la Terre). Le point rouge dans le coin inférieur gauche est son centre. La ligne en pointillés désigne le rayon à un point de la surface. Sa direction "vers le haut" y est représentée par un segment noir: elle est, par définition, perpendiculaire à l'ellipse en ce point. En raison de l'excentricité exagérée, il est facile de voir que "up" n'est pas parallèle au rayon.

Dans notre terminologie, t est lié à l'angle fait par le rayon à l'horizontale et f est l'angle fait par ce segment noir. (Notez que n'importe quel point de la surface peut être vu de cette perspective. Cela nous permet de limiter à la fois t et f entre 0 et 90 degrés; leurs cosinus et sinus seront positifs, donc nous n'avons pas à nous soucier des négatifs racines carrées dans les formules.)

L'astuce consiste à convertir la paramétrisation t en une unité en termes de f , car en termes de t le rayon R est facile à calculer (via le théorème de Pythagore). Son carré est la somme des carrés des composantes du point,

R(t)^2 = a^2 cos(t)^2 + b^2 sin(t)^2.

Pour effectuer cette conversion, nous devons relier la direction "ascendante" f au paramètre t . Cette direction est perpendiculaire à la tangente de l'ellipse. Par définition, une tangente à une courbe (exprimée en vecteur) est obtenue en différenciant sa paramétrisation:

Tangent(t) = d/dt (a cos(t), b sin(t)) = (-a sin(t), b cos(t)).

(La différenciation calcule le taux de changement. Le taux de changement de notre position lorsque nous parcourons la courbe est, bien sûr, notre vitesse , et cela pointe toujours le long de la courbe.)

Faites-le pivoter de 90 degrés dans le sens des aiguilles d'une montre pour obtenir la perpendiculaire, appelée vecteur "normal":

Normal(t) = (b cos(t), a sin(t)).

La pente de ce vecteur normal, égale à (a sin (t)) / (b cos (t)) ("montée sur course"), est aussi la tangente de l'angle qu'il fait à l'horizontale, d'où

tan(f) = (a sin(t)) / (b cos(t)).

De manière équivalente,

(b/a) tan(f) = sin(t) / cos(t) = tan(t).

(Si vous avez une bonne idée de la géométrie euclidienne, vous pouvez obtenir cette relation directement à partir de la définition d'une ellipse sans passer par un trig ou un calcul, simplement en reconnaissant que les expansions horizontales et verticales combinées de a et b respectivement ont pour effet de changer toutes les pentes par ce facteur b / a .)

Regardez à nouveau la formule pour R (t) ^ 2: nous connaissons a et b - ils déterminent la forme et la taille de l'ellipse - nous n'avons donc qu'à trouver cos (t) ^ 2 et sin (t) ^ 2 en termes de f , ce que l'équation précédente nous permet de faire facilement:

cos(t)^2 = 1/(1 + tan(t)^2) 
         = 1 / (1 + (b/a)^2 tan(f)^2) 
         = a^2 / (a^2 + b^2 tan(f)^2);
sin(t)^2 = 1 - cos(t)^2 
         = b^2 tan(f)^2 / (a^2 + b^2 tan(f)^2).

(Lorsque tan (f) est infini, nous sommes au pôle, il suffit donc de définir f = t dans ce cas.)

C'est la connexion dont nous avons besoin. Substituez ces valeurs pour cos (t) ^ 2 et sin (t) ^ 2 dans l'expression de R (t) ^ 2 et simplifiez pour obtenir

R(f)^2 = ( a^4 cos(f)^2 + b^4 sin(f)^2 ) / ( a^2 cos(f)^2 + b^2 sin(f)^2 ).

Une simple transformation montre que cette équation est la même que celle trouvée sur Wikipédia. Parce que a ^ 2 b ^ 2 = (ab) ^ 2 et (a ^ 2) ^ 2 = a ^ 4,

R(f)^2 = ( (a^2 cos(f))^2 + (b^2 sin(f))^2 ) / ( (a cos(f))^2 + (b sin(f))^2 )

Il est intéressant de constater que ma solution analphabète en mathématiques a fait le travail avec 5 minutes de réflexion et de codage, le facteur d'aplatissement ne devrait-il pas être pris en compte plutôt qu'un modèle elliptique parfait?

        double pRad = 6356.7523142;
        double EqRad = 6378.137;                      
        return pRad + (90 - Math.Abs(siteLatitude)) / 90 * (EqRad - pRad); 

C'est du moins la formule que j'ai trouvée au Centre d'analyse et d'évaluation des données des États-Unis (DAAC) pour le wiki du programme de modernisation du calcul haute performance (HPCMP) du ministère de la Défense (DoD) . Cela dit qu'ils ont beaucoup emprunté à l' entrée de Wikipédia . Pourtant, le fait qu'ils aient conservé cette formule devrait compter pour quelque chose.