Obliczanie wschodów i zachodów Słońca

ebvalaim 2015-12-22 1 Comment

Dziś przesilenie zimowe - najkrótszy dzień w roku, a najdłuższa noc. Nie każdemu jednak znany jest fakt, że mimo iż dni będą już coraz dłuższe, to jednak Słońce nadal będzie wschodziło coraz później. Jest to spowodowane kształtem orbity Ziemi, który nie jest idealnie okrągły, a eliptyczny. Teraz zbliża się moment, w którym Ziemia będzie najbliżej Słońca w ciągu swojej rocznej wędrówki, a przez to również porusza się szybciej niż zwykle. To z kolei opóźnia południe słoneczne każdego dnia o kilkanaście sekund, powodując, że godziny wschodów są wciąż coraz późniejsze.

Kiedy zatem mamy najpóźniejszy wschód i najwcześniejszy zachód, jeśli nie w przesilenie? Mógłbym pewnie sprawdzić gdzieś w internecie, ale po co, skoro mogę policzyć sam przy pomocy komputera ;)

Do realizacji zadania wybrałem Haskella - głównie dlatego, że wciąż słabo go ogarniam, a jest to dla mnie bardzo ciekawy i zmieniający sposób myślenia język. Zdecydowałem więc, że warto go poćwiczyć.

Etap pierwszy - znalezienie odpowiednich równań

Jakkolwiek sporo zależności da się wyprowadzić z podstawowych zasad, wiele wynika z faktów obserwacyjnych, więc bez internetu się nie obejdzie. Np. każdy wie, że równonoc wiosenna wypada 21 marca (co, nawiasem mówiąc, przestało być prawdą ok. roku 2000 - teraz przypada ona 20 marca, wróci na 21 ok. roku 2100), jednak mniej oczywiste jest już, jak daleko od Słońca znajduje się wtedy Ziemia i jak szybko się porusza. Tyle wystarczyłoby, żeby liczyć długość dnia z dokładnością do kilku minut, jednak skoro chcemy liczyć takie szczegóły jak przesuwanie się południa słonecznego, trzeba to zrobić porządniej.

Znalezienie danych nie jest trudne - praktycznie od razu Google odsyła do artykułu na Wikipedii. Jest tam wypisane praktycznie krok po kroku co trzeba policzyć, choć oczywiście warto wiedzieć, o co w tym wszystkim chodzi ;) W szczególności dość ważne są pojęcia rektascensji i deklinacji.

Już dawno temu ludzie wpadli na pomysł, że przydałby im się system opisywania położeń gwiazd i planet na niebie. Problem w tym, że niebo pozornie cały czas się obraca (co, jak wiemy, jest efektem ruchu obrotowego Ziemi). Pojawił się wobec tego pomysł wprowadzenia układu współrzędnych niebieskich niezależnych od tego ruchu, nieruchomych względem gwiazd. Idea jest bardzo prosta.

Po pierwsze, wprowadzamy na niebie odpowiednik ziemskiej szerokości geograficznej, zwany deklinacją. Deklinacja +90° odpowiada północnemu biegunowi niebieskiemu, tj. punktowi, w który celuje oś Ziemi "wychodząc" z bieguna północnego. Deklinacja -90° to południowy biegun niebieski, 0° to rzut równika ziemskiego na sferę niebieską. To załatwia jedną współrzędną.

Druga współrzędna to rektascensja, odpowiednik długości geograficznej. Tak, jak na Ziemi mamy arbitralnie wybrany południk zerowy w Greenwich, tak na niebie rektascensji 0 odpowiada tzw. punkt Barana. Jest to punkt przecięcia ekliptyki (płaszczyzny orbity Ziemi, czy też - względem Ziemi - rocznego toru Słońca po niebie) z równikiem niebieskim. Przez ten punkt Słońce przechodzi w momencie równonocy wiosennej. Po przeciwnej stronie nieba mamy drugie przecięcie ekliptyki z równikiem - punkt Wagi.

Dla odmiany, rektascensję liczy się nie w stopniach, ale w godzinach, minutach i sekundach i przyjmuje ona wartości od 0h do 24h. Ma to swoje uzasadnienie - rektascensja oznacza bowiem, ile czasu mija, nim dany punkt na równiku znajdzie się w tym samym punkcie nieba, w którym w danym momencie jest punkt Barana.

Jak już wspomniałem, odzwierciedleniem ruchu orbitalnego Ziemi jest pozorny ruch Słońca po ekliptyce. Znając współrzędne niebieskie (rektascensję i deklinację) Słońca, współrzędne geograficzne miejsca obserwacji i wiedząc, który południk niebieski góruje w danym momencie możemy obliczyć, kiedy Słońce wzejdzie, a kiedy zajdzie.

Dokładnie tych danych dostarczają nam równania z Wikipedii. Po przeliczeniu daty i godziny na tzw. datę juliańską (Julian Date), otrzymujemy współrzędne ekliptyczne Słońca ($\lambda$ i $\beta$ na Wiki), które następnie przeliczamy na rektascensję ($\alpha$) i deklinację ($\delta$), korzystając z kąta nachylenia osi ziemskiej do płaszczyzny orbity ($\epsilon$). Niejako za darmo otrzymujemy przy tym odległość Ziemi od Słońca w danym momencie ;)

To, co bardzo się przyda do obliczenia wschodów i zachodów, to wysokość Słońca nad horyzontem (która jest 0 podczas wschodu i zachodu). Po te wartości Wikipedia odsyła nas do artykułów dot. kąta zenitalnego oraz azymutu. Najprościej rzecz ujmując, te dwa kąty tworzą układ współrzędnych sferycznych, w którym to horyzont pełni rolę równika, zenit jest biegunem północnym, nadir - południowym, a kierunek północny wskazuje południk 0.

Oczywiście to wszystko jest nie do obliczenia, jeśli nie znamy orientacji Ziemi w przestrzeni. Pojawiają się tu takie pojęcia, jak kąt godzinny i czas gwiazdowy. Czas gwiazdowy to nic innego, jak rektascensja aktualnie górującego południka - gdy góruje akurat punkt Barana, mamy czas gwiazdowy 0. Mamy więc czas gwiazdowy Greenwich (czyli czas gwiazdowy na południku Greenwich), który można obliczyć na podstawie daty i godziny, oraz lokalny czas gwiazdowy, który jest po prostu czasem Greenwich skorygowanym o przesunięcie wynikające z długości geograficznej punktu obserwacji. Gdy do tego wszystkiego doliczyć rektascensję interesującego nas punktu (w tym wypadku - Słońca), otrzymujemy kąt godzinny dla tego punktu.

Kąt godzinny określa, jak daleko południk przechodzący przez dany obiekt jest od górowania. Znając go i deklinację obiektu, możemy już obliczyć azymut i wysokość (co sprowadza się do przeliczenia kilku obrotów na sferze niebieskiej). Gotowe równania są podane na Wiki (choć w momencie, gdy je znalazłem, był w nich drobny błąd wynikający z pomieszania terminologii - pochwalę się tu, że go wykryłem i poprawiłem ;)).

Mogąc obliczać wysokość Słońca nad horyzontem dla danego momentu i lokalizacji jesteśmy już tylko o krok od wschodów i zachodów - wystarczy znaleźć momenty, w których wysokość jest 0.

Etap drugi - kodowanie

Mając ogarniętą teorię, możemy wziąć się do pisania kodu. Gotowy program znajduje się tutaj.

W pierwszym podejściu zacząłem od utworzenia własnych typów danych do przechowywania daty i godziny. Dobrzy programiści wiedzą jednak, że obsługa dat, obok kryptografii, jest jednym z przysłowiowych trudnych problemów (nieregularność długości miesięcy to tylko pierwsza przeszkoda - gdy wchodzą różne standardy mierzenia czasu, dopiero robi się ciekawie), których nie należy rozwiązywać samodzielnie. Tak więc postanowiłem poszukać gotowych rozwiązań i odnalazłem moduł Data.Time.

Pierwszy krok, czyli obliczenie daty juliańskiej, moduł ten załatwia praktycznie za nas. Co prawda oblicza tylko część całkowitą (tj. numer dnia) i w wersji tzw. zmodyfikowanej (liczonej od innego dnia i od północy, zamiast od południa), ale doliczenie reszty jest proste - wystarczy przeliczyć godzinę na sekundy, które minęły od północy i podzielić przez 86400, po czym przesunąć wynik o pół dnia. Dalej należy przeliczyć to na tzw. epokę J2000, co sprowadza się do odjęcia pewnej stałej.

Kolejne funkcje obliczają parametry ruchu Ziemi po orbicie i długość ekliptyczną (dla Słońca ekliptyczna szerokość jest zawsze 0, co trochę ułatwia):

Wszystkie funkcje jako parametr przyjmują datę przeliczoną na epokę J2000.

Teraz jesteśmy gotowi do obliczenia rektascensji i deklinacji Słońca:

Ostatnią rzeczą, którą możemy obliczyć nie znając położenia obserwatora, jest czas gwiazdowy Greenwich:

Dalej potrzebujemy już lokalizacji, więc przyda się jakiś typ danych, który go będzie przechowywał:

Możemy teraz obliczyć lokalny czas gwiazdowy i kąt godzinny Słońca:

A stąd już prosta droga do kąta zenitalnego (i tym samym wysokości nad horyzontem) oraz azymutu:

Teraz pytanie, jak sprawdzić, kiedy wysokość wynosi 0.

W pierwszym podejściu zastosowałem naiwny algorytm - dla podanej daty wybieram południe i północ (godziny 0 i 12), po czym sprawdzam, czy w jednym z momentów Słońce jest nad horyzontem, a w drugim pod. Jeśli nie, to stwierdzam dzień/noc polarny/ą, jeśli tak, to dzielę interwał na 2 i wybieram te dwa punkty, dla których Słońce znów jest raz nad, raz pod horyzontem. Powtarzam procedurę dotąd, aż różnica czasu spadnie poniżej zadanego progu i przyjmuję, że jeden z otrzymanych momentów to szukany czas.

Procedura działała całkiem ładnie... dla lokalizacji niezbyt oddalonych od południka 0 i sytuacji bez dnia/nocy polarnego/ej. W innych przypadkach zaczynały się problemy.

Przede wszystkim, długość geograficzna dość mocno wpływa na moment południa słonecznego. Ponieważ wszystko liczymy w czasie UTC, w pobliżu południka Greenwich nie miało to wpływu, ale odpowiednio daleko południe mogło uciec tak, że godziny 0 i 12 obie wypadały w dzień lub w nocy. No i mamy błędnie rozpoznaną sytuację polarną.

To było łatwo skorygować (przesunąć godzinę o długość geograficzną / 15), lecz to nie wszystko. Drugi problem jest mniej oczywisty i ma związek ze wspomnianą na początku eliptycznością orbity. Otóż południe słoneczne nie zawsze wypada o 12 czasu lokalnego; potrafi się wahać o kilkanaście minut do przodu i do tyłu. W efekcie np. w pobliżu nocy polarnej godzina 12 potrafiła wypadać w nocy, mimo że Słońce tego dnia wschodziło. Program rozpoznawał to jednak jako noc polarną.

Na szczęście działające rozwiązanie okazało się niewiele bardziej skomplikowane. Wystarczy najpierw znaleźć południe słoneczne, czyli moment z największą wysokością Słońca. To można zrobić analogicznym algorytmem - zaczynamy np. od lokalnej 11:30 i 12:30 (bo moment południa nie waha się o więcej niż kilkanaście minut, więc nie warto przeszukiwać całej doby i ryzykować, że zahaczymy o południe sąsiedniego dnia). Następnie wybieramy środkowy punkt i zastępujemy nim ten z pary, w którym wysokość Słońca była mniejsza. Powtarzając to aż momenty się zbliżą, osiągniemy południe słoneczne.

Mając południe słoneczne, możemy wybrać momenty oddalone o 12 godzin w obie strony i tym sposobem wrócić do sytuacji z początkowego algorytmu, jednak bez ryzyka, że przegapimy dzień lub noc. W ten sposób znajdziemy wschód i zachód, kiedy to tylko możliwe.

Obliczenie długości dnia to już tylko kwestia odjęcia wschodu od zachodu - nie ma co się nad tym rozpisywać ;)

A oto gotowy algorytm:

Etap trzeci - aplikacja i odpowiedź na początkowe pytanie

Pozostało stworzyć prostą aplikację, która pozwoli na wykorzystanie napisanych funkcji. Skorzystałem w tym celu z optparse-applicative - biblioteka ta pozwala na łatwe utworzenie parsera opcji z linii poleceń. Korzystając z gotowej aplikacji, można już odpowiedzieć kiedy Słońce najpóźniej wschodzi i kiedy najwcześniej zachodzi:

Dla Warszawy, jak się okazuje, najpóźniejszy wschód przypada na 30 grudnia, a najwcześniejszy zachód - 13/14 grudnia :)