Jump to content

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


Wiesiek1952
 Share

Recommended Posts

 

 

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


# ******************************************************************************************
  • Like 3
Link to comment
Share on other sites

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

Edited by WielkiAtraktor
Link to comment
Share on other sites

Posted (edited)
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

Edited by Wiesiek1952
Link to comment
Share on other sites

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ż

 

Link to comment
Share on other sites

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

Link to comment
Share on other sites

  • 1 month later...

Join the conversation

You can post now and register later. If you have an account, sign in now to post with your account.

Guest
Reply to this topic...

×   Pasted as rich text.   Paste as plain text instead

  Only 75 emoji are allowed.

×   Your link has been automatically embedded.   Display as a link instead

×   Your previous content has been restored.   Clear editor

×   You cannot paste images directly. Upload or insert images from URL.

 Share

  • Recently Browsing   0 members

    • No registered users viewing this page.
  • Our picks

    • Migracja Astropolis na nowy serwer - opinie
      Kilka dni temu mogliście przeczytać komunikat o wyłączeniu forum na dobę, co miało związek z migracją na nowy serwer. Tym razem nie przenosiłem Astropolis na większy i szybszy serwer - celem była redukcja dosyć wysokich kosztów (ok 17 tys rocznie za dedykowany serwer z administracją). Biorąc pod uwagę fakt, że płacę z własnej kieszeni, a forum jest organizacją w 100% non profit (nie przynosi żadnego dochodu), nie znalazłem w sobie uzasadnienia na dalsze akceptowanie tych kosztów.
        • Thanks
        • Like
      • 58 replies
    • Droga Mleczna w dwóch gigapikselach
      Zdjęcie jest mozaiką 110 kadrów, każdy po 4 minuty ekspozycji na ISO 400. Wykorzystałem dwa teleskopy Takahashi Epsilon 130D i dwa aparaty Nikon D810A zamocowane na montażu Losmandy G11 wynajętym na miejscu. Teleskopy były ustawione względem siebie pod lekkim kątem, aby umożliwić fotografowanie dwóch fragmentów mozaiki za jednym razem.
        • Love
        • Thanks
        • Like
      • 48 replies
    • Przelot ISS z ogniskowej 2350 mm
      Cześć, po kilku podejściach w końcu udało mi się odpowiednio przygotować cały sprzęt i nadążyć za ISS bez stracenia jej ani razu z pola widzenia. Wykorzystałem do tego montaż Rainbow RST-135, który posiada sprzętową możliwość śledzenia satelitów.
      Celestron Edge 9,25" + ZWO ASI183MM. Czas ekspozycji 6 ms na klatkę, końcowy film składa się z grup 40 klatek stackowanych, wyostrzanych i powiększonych 250%.
      W przyszłości chciałbym wrócić do tematu z kamerką ASI174MM, która z barlowem 2x da mi podobną skalę, ale 5-6 razy większą liczbę klatek na sekundę.
      Poniżej film z przelotu, na dole najlepsza klatka.
        • Love
        • Thanks
        • Like
      • 72 replies
    • Big Bang remnant - Ursa Major Arc or UMa Arc
      Tytuł nieco przekorny bo nie chodzi tu oczywiście o Wielki Wybuch ale ... zacznijmy od początku.
       
      W roku 1997 Peter McCullough używając eksperymentalnej kamery nagrał w paśmie Ha długą na 2 stopnie prostą linie przecinajacą niebo.
       
      Peter McCullough na konferencji pokazał fotografię Robertowi Benjamin i obaj byli pod wrażeniem - padło nawet stwierdzenie: “In astronomy, you never see perfectly straight lines in the sky,”
        • Love
        • Thanks
        • Like
      • 17 replies
    • Jeśli coś jest głupie, ale działa, to nie jest głupie - o nietypowych rozwiązaniach sprzętowych
      Sformułowanie, które można znaleźć w internetach jako jedno z "praw Murphy'ego" przyszło mi na myśl, gdy kolejny raz przeglądałem zdjęcia na telefonie z ostatniego zlotu, mając z tyłu głowy najgłośniejszy marsjański temat na forum. Do rzeczy - jakie macie (bardzo) nietypowe patenty na usprawnienie sprzętu astronomicznego bądź jakieś kreatywne improwizacje w razie awarii czy niezabrania jakiegoś elementu sprzętu  Obstawiam, że @HAMAL mógłby samodzielnie wypełnić treścią taki wątek.
        • Haha
        • Like
      • 43 replies
×
×
  • Create New...

Important Information

We have placed cookies on your device to help make this website better. You can adjust your cookie settings, otherwise we'll assume you're okay to continue.