Skocz do zawartości

Algorytmy astronomiczne czyli Python w rachunkach astro. - Cz. IV Precesja i nutacja

Wiesiek1952

Przez Wiesiek1952
w Oprogramowanie

 Udostępnij

Rekomendowane odpowiedzi

Wiesiek1952

Wiesiek1952

Opublikowano
Opublikowano

 

 

Oś obrotu Ziemi wykonuje w przestrzeni oscylacje. Oś zatacza stożek w okresie około 26 000 lat. Aktualnie oś "celuje" w okolice gwiazdy nazywanej przez nas Polarna (Polaris). Za około 12-13 tyś lat będzie w pobliżu gwiazdy Vega. Ten ruch po obwiedni stożka nazywamy precesją. W skutek tego punkt Barana (Aries) czyli punkt przecięcia równika ziemskiego i ekliptyki zmienia swoje położenie o około 50" rocznie.

 

Drugim ruchem osi ziemskiej jest nutacja. To takie małe kółeczka zataczane w trakcie głównego ruchu precesyjnego. Nutacja jest wynikiem (głównie) wpływu Księżyca jej największy składnik ma okres 6798.4 dnia (18.6 lat) co jest równe pełnemu okresowi w którym orbita Księżyca wypełnia swój torus wokół Ziemii. Są jednakże składowe mające bardzo krótki okres nawet poniżej 10 dni. Formuła jak widać poniżej jest dość długa i uwzględnia większość istotnych składowych. Nutacja została podzielona na dwie składowe - równoległą i prostopadłą do płaszczyzny ekliptyki. Mamy więc nutację "w długości" dotyczącą rektascencji i nutację w "nachyleniu" (obliquity) wpływającą na zmianę nachylenia równika w stosunku do ekliptyki.

 

https://pl.wikipedia.org/wiki/Precesja

https://pl.wikipedia.org/wiki/Nutacja

 

Na skutek tych ruchów położenia "stałych" gwiazd na niebie powoli się zmienia zarówno w rektascencji jak i w deklinacji. To właśnie dlatego katalogi gwiazd są tworzone dla danej epoki. Przyjęło się określanie epoki (epoch) w pełnych 50 lat czyli mamy epoch 1950.0, 2000.0, 2050.0.

 

 

Poniższe procedury wyliczają poprawki do położenia ciał niebieskich związane z precesją i nutacją. Obie procedury wykorzystują wcześniej opisane funkcje i klasy Pythona.

 

 

 

  

# ******************************************************************************************

def precession(current_date: day_in_calendar, current_time: present_time, CN: GeoPoint) -> GeoPoint:
    CN2 = GeoPoint()
    godziny = current_time.hour \
              + current_time.minute / 60.0 \
              + current_time.seconds / 3600.0

    local_date = day_in_calendar()
    local_date.day = current_date.day + godziny / 24.0
    local_date.year = current_date.year
    local_date.month = current_date.month
    t = (julian_day(local_date) - 2451545.0) / 36525.0

    zeta = 2306.2181 * t + 0.30188 * t * t + 0.017998 * t * t * t
    zet = 2306.2181 * t + 1.09468 * t * t + 0.018203 * t * t * t
    theta = 2004.3109 * t - 0.42665 * t * t - 0.041833 * t * t * t

    zeta = zeta / 3600.0 * TO_RAD
    zet = zet / 3600.0 * TO_RAD
    theta = theta / 3600.0 * TO_RAD

    # WPCoord CN2 ;

    A = math.cos(CN.latitude) * math.sin(CN.longitude + zeta)
    B = math.cos(theta) * math.cos(CN.latitude) * math.cos(CN.longitude + zeta) \
        - math.sin(theta) * math.sin(CN.latitude)
    C = math.sin(theta) * math.cos(CN.latitude) * math.cos(CN.longitude + zeta) \
        + math.cos(theta) * math.sin(CN.latitude)

    alfa = math.atan2(A, B)
    if CN.latitude < (math.pi / 2.0 - 0.1):
        declination = math.asin(C)
    else:
        declination = math.acos(math.sqrt(A * A + B * B))

    alfa += zet
    CN2.latitude = declination
    CN2.longitude = alfa

    return CN2


# ******************************************************************************************
#   NUTACJA
#   funkcja zwraca nutację w RADIANACH
#   GeoPoint.longitude zwraca nutację w długości
#   GeoPoint.latitude zwraca nutację w szerokosci
# ******************************************************************************************

def nutation(current_date: day_in_calendar, current_time: present_time) -> GeoPoint:
    # **************************************************************************************
    #  D - Mean elongation of the Moon from the Sun
    #  M - Mean Anomaly of the Sun (Earth)
    #  M1 - Mean Anomaly of tha Moon
    #  F - Moon's argument of latitude
    #  OM - Longitude of Moon's ascending node
    #  JD - Julian Day
    #  JDE - Julian Day Ephemeris time
    #
    #  nutdl - nutacja w długości
    #  nuteps - nutacja w szerokości

    nutacja_corr = GeoPoint()
    godziny = current_time.hour \
              + current_time.minute / 60.0 \
              + current_time.seconds / 3600.0

    local_date = day_in_calendar()
    local_date.day = current_date.day + godziny / 24.0
    local_date.year = current_date.year
    local_date.month = current_date.month
    JD = 0.0
    JD = julian_day(local_date)
    t = (JD - 2451545.0) / 36525.0
    t2 = t * t
    t3 = t2 * t

    D = (297.85036 + 445267.111480 * t - 0.0019142 * t2 + t3 / 189474.0) * TO_RAD
    M = (357.52772 + 35999.050340 * t - 0.0001603 * t2 - t3 / 300000.0) * TO_RAD
    M1 = (134.96298 + 477198.867398 * t + 0.0086972 * t2 + t3 / 56250.0) * TO_RAD
    F = (93.27191 + 483202.017538 * t - 0.0036825 * t2 + t3 / 327270.0) * TO_RAD
    OM = (125.04452 - 1934.136261 * t + 0.0020708 * t2 + t3 / 450000.0) * TO_RAD

    nutdl = (-17.1996 - 0.01742 * t) * math.sin(OM) \
            + (-1.3187 - 0.00016 * t) * math.sin(-2.0 * D + 2.0 * F + 2.0 * OM) \
            + (-0.2274 - 0.00002 * t) * math.sin(2.0 * F + 2.0 * OM) \
            + (0.2062 + 0.00002 * t) * math.sin(2.0 * OM) \
            + (0.1426 - 0.00034 * t) * math.sin(M) \
            + (0.0712 + 0.00001 * t) * math.sin(M1) \
            + (-0.0517 + 0.00012 * t) * math.sin(-2.0 * D + M + 2.0 * F + 2.0 * OM) \
            + (-0.0386 - 0.00004 * t) * math.sin(2.0 * F + OM) \
            - 0.0301 * math.sin(M1 + 2.0 * F + 2.0 * OM) \
            + (0.0217 - 0.00005 * t) * math.sin(-2.0 * D - M + 2.0 * F + 2.0 * OM) \
            - 0.0158 * math.sin(-2.0 * D + M1) \
            + (0.0129 + 0.00001 * t) * math.sin(-2.0 * D + 2.0 * F + OM) \
            + 0.0123 * math.sin(-M1 + 2.0 * F + 2.0 * OM) \
            + 0.0063 * math.sin(2.0 * D) \
            + (0.0063 + 0.00001 * t) * math.sin(M1 + OM) \
            + (-0.0059) * math.sin(2.0 * D - M1 + 2.0 * F + 2.0 * OM) \
            + (-0.0058 - 0.00001 * t) * math.sin(- M1 + OM) \
            + (-0.0051) * math.sin(M1 + 2.0 * F + OM) \
            + 0.0048 * math.sin(-2.0 * D + 2.0 * M1) \
            + 0.0046 * math.sin(- 2.0 * M1 + 2.0 * F + OM) \
            + (-0.0038) * math.sin(2.0 * D + 2.0 * F + 2.0 * OM) \
            + (-0.0031) * math.sin(2.0 * M1 + 2.0 * F + 2.0 * OM) \
            + 0.0029 * math.sin(2.0 * M1) \
            + 0.0029 * math.sin(-2.0 * D + M1 + 2.0 * F + 2.0 * OM) \
            + 0.0026 * math.sin(2.0 * F) \
            + (-0.0022) * math.sin(-2.0 * D + 2.0 * F) \
            + 0.0021 * math.sin(-M1 + 2.0 * F + OM) \
            + (0.0017 - 0.00001 * t) * math.sin(2.0 * M) \
            + 0.0016 * math.sin(2.0 * D - M1 + OM) \
            + (-0.0016 + 0.00001 * t) * math.sin(-2.0 * D + 2.0 * M + 2.0 * F + 2.0 * OM) \
            + (-0.0015) * math.sin(M + OM) \
            + (-0.0013) * math.sin(-2.0 * D + M1 + OM) \
            + (-0.0012) * math.sin(-M + OM) \
            + 0.0011 * math.sin(2.0 * M1 - 2.0 * F) \
            + (-0.0010) * math.sin(2.0 * D - M1 + 2.0 * F + OM) \
            + (-0.0008) * math.sin(2.0 * D + M1 + 2.0 * F + 2.0 * OM) \
            + 0.0007 * math.sin(M + 2.0 * F + 2.0 * OM) \
            + (-0.0007) * math.sin(-2.0 * D + M + M1) \
            + (-0.0007) * math.sin(- M + 2.0 * F + 2.0 * OM) \
            + (-0.0007) * math.sin(2.0 * D + 2.0 * F + OM) \
            + 0.0006 * math.sin(2.0 * D + M1) \
            + 0.0006 * math.sin(-2.0 * D + 2.0 * M1 + 2.0 * F + 2.0 * OM) \
            + 0.0006 * math.sin(-2.0 * D + M1 + 2.0 * F + OM) \
            - 0.0006 * math.sin(2.0 * D - 2.0 * M1 + OM) \
            - 0.0006 * math.sin(2.0 * D + OM) \
            + 0.0005 * math.sin(-M + M1) \
            - 0.0005 * math.sin(-2.0 * D - M + 2.0 * F + OM) \
            - 0.0005 * math.sin(-2.0 * D + OM) \
            - 0.0005 * math.sin(2.0 * M1 + 2.0 * F + OM) \
            + 0.0004 * math.sin(-2.0 * D + 2.0 * M1 + OM) \
            + 0.0004 * math.sin(-2.0 * D + M + 2.0 * F + OM) \
            + 0.0004 * math.sin(M1 - 2.0 * F) \
            - 0.0004 * math.sin(-D + M1) \
            - 0.0004 * math.sin(-2.0 * D + M) \
            - 0.0004 * math.sin(D) \
            + 0.0003 * math.sin(M1 + 2.0 * F) \
            + (-0.0003) * math.sin(-2.0 * M1 + 2.0 * F + 2.0 * OM) \
            + (-0.0003) * math.sin(-D - M + M1) \
            + (-0.0003) * math.sin(M + M1) \
            + (-0.0003) * math.sin(-M + M1 + 2.0 * F + 2.0 * OM) \
            + (-0.0003) * math.sin(2.0 * D - M - M1 + 2.0 * F + 2.0 * OM) \
            + (-0.0003) * math.sin(3.0 * M1 + 2.0 * F + 2.0 * OM) \
            + (-0.0003) * math.sin(2.0 * D - M + 2.0 * F + 2.0 * OM)

    #     -------

    nuteps = (9.2025 + 0.00089 * t) * math.cos(OM) \
             + (0.5736 - 0.00031 * t) * math.cos(-2.0 * D + 2.0 * F + 2.0 * OM) \
             + (0.0977 - 0.00005 * t) * math.cos(2.0 * F + 2.0 * OM) \
             + (-0.0895 + 0.00005 * t) * math.cos(2.0 * OM) \
             + (0.0054 - 0.00001 * t) * math.cos(M) \
             - 0.0007 * math.cos(M1) \
             + (0.0224 - 0.00006 * t) * math.cos(-2.0 * D + M + 2.0 * F + 2.0 * OM) \
             + 0.0200 * math.cos(2.0 * F + OM) \
             + (0.0129 - 0.00001 * t) * math.cos(M1 + 2.0 * F + 2.0 * OM) \
             + (-0.0095 + 0.00003 * t) * math.cos(-2.0 * D - M + 2.0 * F + 2.0 * OM) \
             - 0.0070 * math.cos(-2.0 * D + 2.0 * F + OM) \
             - 0.0053 * math.cos(-M1 + 2.0 * F + 2.0 * OM) \
             - 0.0033 * math.cos(M1 + OM) \
             + 0.0026 * math.cos(2.0 * D - M1 + 2.0 * F + 2.0 * OM) \
             + 0.0032 * math.cos(- M1 + OM) \
             + 0.0027 * math.cos(M1 + 2.0 * F + OM) \
             - 0.0024 * math.cos(- 2.0 * M1 + 2.0 * F + OM) \
             + 0.0016 * math.cos(2.0 * D + 2.0 * F + 2.0 * OM) \
             + 0.0013 * math.cos(2.0 * M1 + 2.0 * F + 2.0 * OM) \
             - 0.0012 * math.cos(-2.0 * D + M1 + 2.0 * F + 2.0 * OM) \
             - 0.0010 * math.cos(-M1 + 2.0 * F + OM) \
             - 0.0008 * math.cos(2.0 * D - M1 + OM) \
             + 0.0007 * math.cos(-2.0 * D + 2.0 * M + 2.0 * F + 2.0 * OM) \
             + 0.0009 * math.cos(M + OM) \
             + 0.0007 * math.cos(-2.0 * D + M1 + OM) \
             + 0.0006 * math.cos(-M + OM) \
             + 0.0005 * math.cos(2.0 * D - M1 + 2.0 * F + OM) \
             + 0.0003 * math.cos(2.0 * D + M1 + 2.0 * F + 2.0 * OM) \
             - 0.0003 * math.cos(M + 2.0 * F + 2.0 * OM) \
             + 0.0003 * math.cos(- M + 2.0 * F + 2.0 * OM) \
             + 0.0003 * math.cos(2.0 * D + 2.0 * F + OM) \
             - 0.0003 * math.cos(-2.0 * D + 2.0 * M1 + 2.0 * F + 2.0 * OM) \
             - 0.0003 * math.cos(-2.0 * D + M1 + 2.0 * F + OM) \
             + 0.0003 * math.cos(2.0 * D - 2.0 * M1 + OM) \
             + 0.0003 * math.cos(2.0 * D + OM) \
             + 0.0003 * math.cos(-2.0 * D - M + 2.0 * F + OM) \
             + 0.0003 * math.cos(-2.0 * D + OM) \
             + 0.0003 * math.cos(2.0 * M1 + 2.0 * F + OM)
    nutdl = nutdl / 3600.0 * TO_RAD
    nuteps = nuteps / 3600.0 * TO_RAD

    nutacja_corr.latitude = nuteps
    nutacja_corr.longitude = nutdl

    return nutacja_corr


# ******************************************************************************************
  • Lubię 3
Odnośnik do komentarza
Udostępnij na innych stronach
WielkiAtraktor

WielkiAtraktor

Opublikowano
Opublikowano (edytowane)
4 minuty temu, AviatorL napisał:

Czy nie brakuje klasy GeoPoint?

Pewnie była w poprzednich częściach wątku. @Wiesiek1952 ponawiam sugestię wrzucenia tego na GitHuba, byłoby jedno repozytorium ze wszystkim, do którego można co jakiś czas dodawać odkurzony i sprawdzony kod ;)

Edytowane przez WielkiAtraktor
Odnośnik do komentarza
Udostępnij na innych stronach
Wiesiek1952

Wiesiek1952

Opublikowano
  • Autor
Opublikowano (edytowane)
48 minut temu, AviatorL napisał:

Czy nie brakuje klasy GeoPoint?

:-) 

 

Faktycznie brakuje. Wszystkie klasy używane do teraz i w najbliższych odcinkach poniżej. One są wszystkie bardzo prymitywne, właściwie to kontenery ale python nie używa (w zasadzie) niczego innego poza klasami. Tym co lubią ponarzekać i zrobiliby to lepiej - tak wiem, że kod jest nieoptymalny, można to zrobić inaczej, można to zrobić lepiej itp. itd:

 

 

  

class GeoPoint:
    def __init__(self, latitude: float = 0.0, longitude: float = 0.0) -> object:
        self.latitude = latitude
        self.longitude = longitude


class GeoPointPolar:
    def __init__(self, direction: float = 0.0, distance: float = 0.0):
        self.direction = direction
        self.distance = distance

class day_in_calendar:
    def __init__(self, year: int = 0, month: int = 0, day: float = 0):
        self.year = year
        self.month = month
        self.day = day


class present_time:
    def __init__(self, hour: int = 0, minute: int = 0, seconds: float = 0):
        self.hour = hour
        self.minute = minute
        self.seconds = seconds

class sun_efemeries:
    def __init__(self, sha: float = 0, declination: float = 0, semi_diameter: float = 0):
        self.sha = sha
        self.declination = declination
        self.semi_diameter = semi_diameter

class star_ephemeries:
    def __init__(self, rektascencja: float = 0.0, deklinacja: float = 0.0, nazwa: str =''):
        self.nazwa = nazwa
        self.rektascencja = rektascencja
        self.deklinacja = deklinacja


class moon_ephemeries:
    def __init__(self, rektascencja: float = 0.0, deklinacja: float = 0.0, paralaksa: float = 0, sd: float = 0):
        self.rektascencja = rektascencja
        self.deklinacja = deklinacja
        self.paralaksa = paralaksa
        self.sd = sd
Edytowane przez Wiesiek1952
Odnośnik do komentarza
Udostępnij na innych stronach
Wiesiek1952

Wiesiek1952

Opublikowano
  • Autor
Opublikowano
44 minuty temu, WielkiAtraktor napisał:

Pewnie była w poprzednich częściach wątku. @Wiesiek1952 ponawiam sugestię wrzucenia tego na GitHuba, byłoby jedno repozytorium ze wszystkim, do którego można co jakiś czas dodawać odkurzony i sprawdzony kod ;)

 

 

Hmmmm... nie wiem. Nigdy niczego nie umieszczałem na githubie - nawet nie wiem jak to zrobić. Może na koniec wrzucę w postaci zipa tutaj ale zobaczę na githubie też

 

Odnośnik do komentarza
Udostępnij na innych stronach
zbuffer

zbuffer

Opublikowano
Opublikowano

Moim zdaniem lepszym rozwiązaniem niż kod Python'a wrzucony tylko na GitHub byłoby skorzystanie z Jupyter Notebook's na jakimś darmowym serwisie. Nie używałem jeszcze ale dość popularny jest Colab od Google. Największa przewaga Jupyter Notebooks jest taka, że można napisać całą teorię, łącznie z formułami matematycznymi do każdej z funkcji, osadzone bezpośrednio w dokumencie. W ten sposób nie byłby to tylko kod ale baza wiedzy na temat obliczeń astronomicznych. Można do tego złączać obrazy no i generować wykresy osadzone w tym samym dokumencie. Dodatkowa przewaga jest taka, że używając Colab'a lub podobnych serwisów nie ma konieczności instalacji pythona ani żadnych bibliotek, bo wszystko działa przez przeglądarkę.

Odnośnik do komentarza
Udostępnij na innych stronach
  • 1 miesiąc temu...

Dołącz do dyskusji

Możesz dodać zawartość już teraz a zarejestrować się później. Jeśli posiadasz już konto, zaloguj się aby dodać zawartość za jego pomocą.

Gość
Dodaj odpowiedź do tematu...

×   Wklejono zawartość z formatowaniem.   Usuń formatowanie

  Dozwolonych jest tylko 75 emoji.

×   Odnośnik został automatycznie osadzony.   Przywróć wyświetlanie jako odnośnik

×   Przywrócono poprzednią zawartość.   Wyczyść edytor

×   Nie możesz bezpośrednio wkleić grafiki. Dodaj lub załącz grafiki z adresu URL.

×
 Udostępnij
Tematy
×
×
  • Dodaj nową pozycję...

Powiadomienie o plikach cookie

Umieściliśmy na Twoim urządzeniu pliki cookie, aby pomóc Ci usprawnić przeglądanie strony. Możesz dostosować ustawienia plików cookie, w przeciwnym wypadku zakładamy, że wyrażasz na to zgodę.