Jak obliczyć promień ziemi na danej szerokości geograficznej?

(Widzę, że na wikipedii jest równanie, które robi dokładnie to, o co proszę, ale nie ma żadnych odniesień. Nie mam sposobu na potwierdzenie ważności tego równania!)

Rozumiem już różnicę między geocentryczną szerokością geograficzną a geodezyjną.

Zakładając, że znane są pół-duże ai pół-mniejsze, bpodano promienie. Jak obliczyć promień na danej szerokości geograficznej?

Potrzebuję jakiegoś potwierdzenia eksperta (wyprowadzenie, link do wyprowadzenia, potwierdzenie od eksperta, wyjaśnienie itp.).

To pytanie zakłada elipsoidalny model ziemi. Jego powierzchnię odniesienia uzyskuje się poprzez obrócenie elipsy wokół jej mniejszej osi (wykreślonej pionowo zgodnie z konwencją). Taka elipsa jest po prostu okręgiem, który został rozciągnięty poziomo o współczynnik a i pionowo o współczynnik b . Używając standardowej parametryzacji koła jednostki,

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

(która definiuje cosinus i sinus), uzyskujemy parametryzację

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

(Dwa elementy tej parametryzacji opisują podróż wokół krzywej: określają, we współrzędnych kartezjańskich, naszą lokalizację w „czasie” t .)

Szerokości geodezyjnej , F , dla dowolnego punktu, to kąt, że „w górę” kieruje do płaszczyzny równikowej. Gdy a różni się od b , wartość f różni się od wartości t (z wyjątkiem wzdłuż równika i na biegunach).

Na tym zdjęciu niebieska krzywa jest jedną ćwiartką takiej elipsy (znacznie przesadzona w porównaniu z ekscentrycznością Ziemi). Czerwona kropka w lewym dolnym rogu to jej środek. Linia przerywana oznacza promień do jednego punktu na powierzchni. Jego „górny” kierunek jest pokazany za pomocą czarnego segmentu: z definicji jest on prostopadły do ​​elipsy w tym punkcie. Ze względu na przesadną mimośrodowość łatwo zauważyć, że „góra” nie jest równoległa do promienia.

W naszej terminologii t jest związany z kątem wykonanym przez promień do poziomu, a f jest kątem wykonanym przez ten czarny segment. (Zauważ, że z tej perspektywy można patrzeć na dowolny punkt na powierzchni. To pozwala nam ograniczyć zarówno t, jak i f, do wartości od 0 do 90 stopni; ich cosinus i sinus będą dodatnie, więc nie musimy się martwić o negatywne pierwiastki kwadratowe we wzorach).

Sztuką jest konwersja z parametryzacji t na jeden w kategoriach f , ponieważ w kategoriach t promień R jest łatwy do obliczenia (za pomocą twierdzenia Pitagorasa). Jego kwadrat jest sumą kwadratów składników punktu,

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

Aby dokonać tej konwersji, musimy powiązać kierunek „w górę” f z parametrem t . Ten kierunek jest prostopadły do ​​stycznej elipsy. Z definicji styczną do krzywej (wyrażoną jako wektor) uzyskuje się poprzez różnicowanie jej parametryzacji:

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

(Różnicowanie oblicza szybkość zmian. Szybkość zmiany naszej pozycji podczas podróży po krzywej jest oczywiście naszą prędkością i zawsze wskazuje na krzywą.)

Obróć to w prawo o 90 stopni, aby uzyskać prostopadłą, zwaną wektorem „normalnym”:

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

Nachylenie tego wektora normalnego, równe (sin (t)) / (b cos (t)) („wzniesienie po przebiegu”), jest również styczną kąta, jaki robi do poziomu, skąd

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

Równoważnie

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

(Jeśli masz dobry wgląd w geometrię euklidesową, możesz uzyskać ten związek bezpośrednio z definicji elipsy bez przechodzenia przez dowolny wyzwalacz lub rachunek różniczkowy, po prostu rozpoznając, że połączone poziome i pionowe rozwinięcia odpowiednio przez a i b powodują zmianę wszystkie zbocza tego współczynnika b / a .)

Ponownie przyjrzeć wzoru na R (t) ^ 2: wiemy, jak i B - oni określają kształt i wielkość elipsy - więc musimy tylko znaleźć cos (t) ^ 2 i sin (t) ^ 2 w kategoriach f , które powyższe równanie pozwala nam łatwo zrobić:

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

(Gdy tan (f) jest nieskończony, jesteśmy na biegunie, więc w takim przypadku po prostu ustaw f = t .)

To jest połączenie, którego potrzebujemy. Zamień te wartości na cos (t) ^ 2 i sin (t) ^ 2 w wyrażenie na R (t) ^ 2 i uproś, aby uzyskać

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

Prosta transformacja pokazuje, że to równanie jest takie samo jak w Wikipedii. Ponieważ a ^ 2 b ^ 2 = (ab) ^ 2 i (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 )

Ciekawe, że moje niepiśmienne rozwiązanie matematyczne wykonało zadanie po 5 minutach namysłu i kodowania, czy nie należy brać pod uwagę czynnika spłaszczania zamiast idealnego modelu eliptycznego?

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

Przynajmniej taką formułę znalazłem w amerykańskim Centrum analizy i oceny danych (DAAC) dla wiki Departamentu Obrony (DoD) High Performance Computing Modernization Program (HPCMP) wiki . Mówi się, że mocno pożyczyli od wpisu Wikipedii . Mimo to fakt, że zachowali tę formułę, powinien na coś liczyć.