Algorytmy astronomiczne czyli Python w rachunkach astro. - Cz. VII Współrzędne Ziemi czyli gdzie jest Słońce.

Wiesiek1952 · 16 Czerwca 2022

Tytuł brzmi nieco przewrotnie ale taka prawda. Do dokładnego wyliczenia położenia słońca na naszym niebie liczymy najpierw położenie Ziemi w układzie heliocentrycznym i dopiero z tego wyliczamy gdzie jest słońce. Niestety tory wszystkich planet podlegają wzajemnym perturbacjom, Jowisz zakłóca orbity wszystkich innych planet a te w rewanżu zakłócają drogę Jowisza.

w 1982 r opracowano metodę zwaną VSOP82 (Variations Seculiares des Orbiter Planetaires). Pełne tablice współczynników dla poszczególnych składowych harmonicznych ruchu planet pozwalające na wyliczenie dość dokładnych orbit liczą sobie po kilka tysięcy wierszy. W 1987 r Bretagnon i Franco uprościli nieco obliczenia skracając i odrzucając elementy, które dokładność poprawiają w tak znikomy sposób, że po prostu szkoda czasu na ich liczenie.

Poniżej dane dla orbity Ziemi i funkcja wyliczająca parametry położenia naszej planety, która jest adaptacją Jeana Meeus'a metody VSOP87 i jest wystarczająco dokładna dla prawie wszystkich praktycznych zastosowań. Funkcja heliocentric_coordinates_earth() wylicza to co ma w nazwie. W celu otrzymania współrzędnych równikowych trzeba skorzystać z funkcji z rozdziału VI niniejszej serii)

Na stronie http://neoprogrammics.com/vsop87/source_code_generator_tool/ jest generator kodu wyliczającego orbity wg VSOP87 dla wszystkich planet i dla różnych języków programowania - niestety nie ma Pythona, jest za to CPP. Kod obejmuje okres: Ephemeris Date Span : 2000 BC < Date < 6000 AD. (w przyciskach u góry strony wpisujemy np Venus, czy Earth i podobnie w dwóch pozostałych - np. dla CPP wystarczy wpisać C)

# ******************************************************************************************
#    PLIK earth_data.py 
#    based on Meeus books and publications
#    using VSOP82 methods (Variations Seculiares des Orbites Planetaires)
#    The VSOP82 consists of long series of periodic term for each major planet.
#    Using this method, one obtains the values of the following quantities
#    for the osculating orbit.
#                   a = semimajor axis ot the orbit
#                   lambda = mean longitude of the planet
#                   h = e * sin(pi)
#                   k = e * cos(pi)
#                   p = sin(1/2 * i) * sin(omega)
#                   q = sin(1/2 * i) * cos)omega
#           where   e - orbital eccentricity,
#                   pi - longitude of the perihelion
#                   i - inclination
#                   omega - longitude of ascending node
#    Bretagnon and Francou constructed VSOP87  - used here.
#                   L - ecliptical longitude
#                   B - ecliptical latitude
#                   R - radius vector (= distance to the Sun)
# ******************************************************************************************
#
# ******************************************************************************************


Earth_L0 = [[175347046, 0, 0],
            [3341656, 4.6692568, 6283.0758500],
            [34894, 4.62610, 12566.15170],
            [3497, 2.7441, 5753.3849],
            [3418, 2.8289, 3.5231],
            [3136, 3.6277, 77713.7715],
            [2676, 4.4181, 7860.4194],
            [2343, 6.1352, 3930.2097],
            [1324, 0.7425, 11506.7698],
            [1273, 2.0371, 529.691],
            [1199, 1.1096, 1577.3435],
            [990, 5.233, 5884.927],
            [902, 2.045, 26.298],
            [857, 3.508, 398.149],
            [780, 1.179, 5223.694],
            [753, 2.533, 5507.553],
            [505, 4.583, 18849.228],
            [492, 4.205, 775.523],
            [357, 2.920, 0.067],
            [317, 5.849, 11790.629],
            [284, 1.899, 796.298],
            [271, 0.315, 10977.079],
            [243, 0.345, 5486.778],
            [206, 4.806, 2544.314],
            [205, 1.869, 5573.143],
            [202, 2.458, 6069.777],
            [156, 0.833, 213.299],
            [132, 3.411, 2942.463],
            [126, 1.083, 20.775],
            [115, 0.645, 0.980],
            [103, 0.636, 4694.003],
            [102, 0.976, 15720.839],
            [102, 4.267, 7.114],
            [99, 6.21, 2146.17],
            [98, 0.68, 155.42],
            [86, 5.98, 161000.69],
            [85, 1.30, 6275.96],
            [85, 3.67, 71430.70],
            [80, 1.81, 17260.15],
            [79, 3.04, 12036.46],
            [75, 1.76, 5088.63],
            [74, 3.50, 3154.69],
            [74, 4.68, 801.82],
            [70, 0.83, 9437.76],
            [62, 3.98, 8827.39],
            [61, 1.82, 7084.90],
            [57, 2.78, 6286.60],
            [56, 4.39, 14143.50],
            [56, 3.47, 6279.55],
            [52, 0.19, 12139.55],
            [52, 1.33, 1748.02],
            [51, 0.28, 5856.48],
            [49, 0.49, 1194.45],
            [41, 5.37, 8429.24],
            [41, 2.40, 19651.05],
            [39, 6.17, 10447.39],
            [37, 6.04, 10213.29],
            [37, 2.57, 1059.38],
            [36, 1.71, 2352.87],
            [36, 1.78, 6812.77],
            [33, 0.59, 17789.85],
            [30, 0.44, 83996.85],
            [30, 2.74, 1349.87],
            [25, 3.16, 4690.48]
            ]

Earth_L1 = [[628331966747, 0, 0],
            [206059, 2.678235, 6283.075850],
            [4303, 2635.1, 12566.1517],
            [425, 1.590, 3.523],
            [119, 5.796, 26.298],
            [109, 2.966, 1577.344],
            [93, 2.59, 18849.23],
            [72, 1.14, 529.69],
            [68, 1.87, 398.15],
            [67, 4.41, 5507.55],
            [59, 2.89, 5223.69],
            [56, 2.17, 155.42],
            [45, 0.40, 796.30],
            [36, 0.47, 775.52],
            [29, 2.65, 7.11],
            [21, 5.34, 0.98],
            [19, 1.85, 5486.78],
            [19, 4.97, 213.30],
            [17, 2.99, 6275.96],
            [16, 0.03, 2544.31],
            [16, 1.43, 2146.17],
            [15, 1.21, 10977.08],
            [12, 2.83, 1748.02],
            [12, 3.26, 5088.63],
            [12, 5.27, 1194.45],
            [12, 2.08, 4694.00],
            [11, 0.77, 553.57],
            [10, 1.30, 6286.60],
            [10, 4.24, 1349.87],
            [9, 2.70, 242.73],
            [9, 5.64, 951.72],
            [8, 5.30, 2352.87],
            [6, 2.65, 9437.76],
            [6, 4.67, 4690.48]
            ]

Earth_L2 = [[52919, 0, 0],
            [8720, 1.0721, 6283.0758],
            [309, 0.867, 12566.152],
            [27, 0.05, 3.52],
            [16, 5.19, 26.30],
            [16, 3.68, 155.42],
            [10, 0.76, 18849.23],
            [9, 2.06, 77713.77],
            [7, 0.83, 775.52],
            [5, 4.66, 1577.34],
            [4, 1.03, 7.11],
            [4, 3.44, 5573.14],
            [3, 5.14, 796.30],
            [3, 6.05, 5507.55],
            [3, 1.19, 242.73],
            [3, 6.12, 529.69],
            [3, 0.31, 398.15],
            [3, 2.38, 553.57],
            [2, 4.38, 5223.69],
            [2, 3.75, 0.98]
            ]

Earth_L3 = [[289, 5.844, 6283.076],
            [35, 0, 0],
            [17, 5.49, 12566.15],
            [3, 5.20, 155.42],
            [1, 4.72, 3.52],
            [1, 5.30, 18849.23],
            [1, 5.97, 242.73]
            ]

Earth_L4 = [[114, 3.142, 0],
            [8, 4.13, 6283.08],
            [1, 3.84, 12566.15]
            ]

Earth_L5 = [[1, 3.14, 0]]

Earth_B0 = [[280, 3.199, 84334.662],
            [102, 5.422, 5507.553],
            [80, 3.88, 5223.69],
            [44, 3.70, 2352.87],
            [32, 4.00, 1577.34]
            ]

Earth_B1 = [[9, 3.90, 5507.55],
            [6, 1.73, 5223.69]
            ]

Earth_R0 = [[100013989, 0, 0],
            [1670700, 3.0984635, 6283.0758500],
            [13956, 3.05525, 12566.15170],
            [3084, 5.1985, 77713.7715],
            [1628, 1.1739, 5753.3849],
            [1576, 2.8469, 7860.4194],
            [925, 5.453, 11506.770],
            [542, 4.564, 3930.210],
            [472, 3.661, 5884.927],
            [346, 0.964, 5507.553],
            [329, 5.900, 5223.694],
            [307, 0.299, 5573.143],
            [243, 4.273, 11790.629],
            [212, 5.847, 1577.344],
            [186, 5.022, 10977.079],
            [175, 3.012, 18849.228],
            [110, 5.055, 5486.778],
            [98, 0.89, 6069.78],
            [86, 5.69, 15720.84],
            [86, 1.27, 161000.69],
            [65, 0.27, 17260.15],
            [63, 0.92, 529.69],
            [57, 2.01, 83996.85],
            [56, 5.24, 71430.70],
            [49, 3.25, 2544.31],
            [47, 2.58, 775.52],
            [45, 5.54, 9437.76],
            [43, 6.01, 6275.96],
            [39, 5.36, 4694.00],
            [38, 2.39, 8827.39],
            [37, 0.83, 19651.05],
            [37, 4.90, 12139.55],
            [36, 1.67, 12036.46],
            [35, 1.84, 2942.46],
            [33, 0.24, 7084.90],
            [32, 0.18, 5088.63],
            [32, 1.78, 398.15],
            [28, 1.21, 6286.60],
            [28, 1.90, 6279.55],
            [26, 4.59, 10447.39]
            ]

Earth_R1 = [[103019, 1.107490, 6283.075850],
            [1721, 1.0644, 12566.1517],
            [702, 3.142, 0],
            [32, 1.02, 18849.23],
            [31, 2.84, 5507.55],
            [25, 1.32, 5223.69],
            [18, 1.42, 1577.34],
            [10, 5.91, 10977.08],
            [9, 1.42, 6275.96],
            [9, 0.27, 5486.78]
            ]

Earth_R2 = [[4359, 5.7846, 6283.0758],
            [124, 5.579, 12566.152],
            [12, 3.14, 0],
            [9, 3.63, 77713.77],
            [6, 1.87, 5573.14],
            [3, 5.47, 18849.23]
            ]

Earth_R3 = [[145, 4.273, 6283.076],
            [7, 3.92, 12566.15]
            ]

Earth_R4 = [[4, 2.56, 6283.08]]

# ******************************************************************************************
#                          END OF FILE
# ******************************************************************************************

# ******************************************************************************************
#    EARTH
# ******************************************************************************************

def heliocentric_coordinates_earth(current_date: day_in_calendar, current_time: present_time)-> sun_efemeries:

    czas = current_time.hour + current_time.minute / 60 + current_time.seconds / 3600
    actual_day = current_date.day
    current_date.day += czas / 24.0
    JDE = julian_day(current_date)
    tau = (JDE - 2451545.0) / 365250.0

    a_sum1 = a_sum2 = a_sum3 = a_sum4 = a_sum5 = a_sum6 = 0.0
    b_sum1 = b_sum2 = b_sum3 = b_sum4 = b_sum5 = 0.0
    r_sum1 = r_sum2 = r_sum3 = r_sum4 = r_sum5 = r_sum6 = 0.0
    tau2 = tau * tau
    tau3 = tau2 * tau
    tau4 = tau3 * tau
    tau5 = tau4 * tau
    for k in range(len(Earth_L0)):
        a_sum1 += Earth_L0[k][0] * cos(Earth_L0[k][1] + Earth_L0[k][2] * tau)
    for k in range(len(Earth_L1)):
        a_sum2 += Earth_L1[k][0] * cos(Earth_L1[k][1] + Earth_L1[k][2] * tau)
    for k in range(len(Earth_L2)):
        a_sum3 += Earth_L2[k][0] * cos(Earth_L2[k][1] + Earth_L2[k][2] * tau)
    for k in range(len(Earth_L3)):
        a_sum4 += Earth_L3[k][0] * cos(Earth_L3[k][1] + Earth_L3[k][2] * tau)
    for k in range(len(Earth_L4)):
        a_sum5 += Earth_L4[k][0] * cos(Earth_L4[k][1] + Earth_L4[k][2] * tau)
    for k in range(len(Earth_L5)):
        a_sum6 += Earth_L5[k][0] * cos(Earth_L5[k][1] + Earth_L5[k][2] * tau)

    for k in range(len(Earth_B0)):
        b_sum1 += Earth_B0[k][0] * cos(Earth_B0[k][1] + Earth_B0[k][2] * tau)
    for k in range(len(Earth_B1)):
        b_sum2 += Earth_B1[k][0] * cos(Earth_B1[k][1] + Earth_B1[k][2] * tau)



    for k in range(len(Earth_R0)):
        r_sum1 += Earth_R0[k][0] * cos(Earth_R0[k][1] + Earth_R0[k][2] * tau)
    for k in range(len(Earth_R1)):
        r_sum2 += Earth_R1[k][0] * cos(Earth_R1[k][1] + Earth_R1[k][2] * tau)
    for k in range(len(Earth_R2)):
        r_sum3 += Earth_R2[k][0] * cos(Earth_R2[k][1] + Earth_R2[k][2] * tau)
    for k in range(len(Earth_R3)):
        r_sum4 += Earth_R3[k][0] * cos(Earth_R3[k][1] + Earth_R3[k][2] * tau)
    for k in range(len(Earth_R4)):
        r_sum5 += Earth_R4[k][0] * cos(Earth_R4[k][1] + Earth_R4[k][2] * tau)

    # heliocentric eliptical longitude L
    L = (round(a_sum1) + round(a_sum2) * tau
         + round(a_sum3) * tau2 + round(a_sum4) * tau3
         + round(a_sum5) * tau4 + round(a_sum6 * tau5)) / 10e+7

    # heliocentric eliptical latitude B
    suma_b = (round(b_sum1) + round(b_sum2) * tau) / 10e+7
    B = -suma_b

    suma_r = (round(r_sum1) + round(r_sum2) * tau
              + round(r_sum3) * tau2 + round(r_sum4) * tau3
              + round(r_sum5) * tau4) / 10e+7
    R = suma_r

    while L < 0:
        L += TWO_PI
    while L > TWO_PI:
        L -= TWO_PI

    T = tau * 10

    theta = L + math.pi
    theta -= 0.09033/3600.0 * TO_RAD
    # conversion to FK5 system
    L1 = (theta * TO_DEG - 1.397 * T - 0.00031 * T * T) * TO_RAD
    while L1 < 0:
        L1 += TWO_PI
    while L1 > TWO_PI:
        L1 -= TWO_PI

    delta_b  =  (0.03916 * (math.cos(L1) - math.sin(L1))) * TO_RAD / 3600.0
    aberation = 20.4898 / R / 3600.0 * TO_RAD
    nutacja = nutation(current_date, current_time)
    sun_apparent_longitude = theta + nutacja.longitude - aberation

    sun_apparent_delination = B + delta_b

    epsilon = ecliptic_obliquity(current_date, current_time) + nutacja.latitude

    efem_sun_value = sun_efemeries()
    efem_sun_value.sha = sun_apparent_longitude
    efem_sun_value.declination = sun_apparent_delination
    efem_sun_value.semi_diameter = R

    current_date.day = actual_day
    return efem_sun_value

Edytowane 16 Czerwca 2022 przez Wiesiek1952

zbuffer · 17 Czerwca 2022

Nie myślałeś o tym żeby do obliczeń używać biblioteki Numpy? Jak wszyscy….

Zamiast tych wszystkich pętli można po prostu użyć odpowiednich funkcji a indeksowanie „samo się robi”, jak w Matlabie.
Nie podoba mi się to mieszanie języka polskiego z angielskim.

Komentuję bo szkoda byłoby żeby Twoja praca i chęć dzielenia się wiedzą nie została doceniona przez sposób podania, który jest mało atrakcyjny i

pewnie odstrasza początkujących. Jak będę miał chwilę to wrzucę jakiś Twój kod do Jupytera i zobaczysz co mam na myśli jeśli go nie używasz.

Wiesiek1952 · 17 Czerwca 2022

Godzinę temu, zbuffer napisał:

Nie myślałeś o tym żeby do obliczeń używać biblioteki Numpy? Jak wszyscy….

Zamiast tych wszystkich pętli można po prostu użyć odpowiednich funkcji a indeksowanie „samo się robi”, jak w Matlabie.

Pole kwadratu można policzyć zwykłym mnożeniem i można użyć całek. Każde dodanie kolejnego poziomu komplikacji tylko zaciemnia obraz. Niestety Numpy dla początkującego to dodatkowa komplikacja i robi się z tego kolejny "black box" a jak ma być "black box" to są inne gotowce, lepsze i wydajniejsze choć zrozumienie tego co robią wcale łatwiejsze nie jest. Dla zaawansowanych programistów może się wydawać "numpy, jak wszyscy..." . - no nie wszyscy... Przerobiłem dość sporo różnych języków programowania od Algol, Basic, Fortran, Pascal, C, C++, Swift, Scratch po "wynalazki" typu Visual Basic i parę innych już całkiem zapomnianych. Jestem całkowicie przekonany, że każde zapakowanie kawałka kodu w jakąś funkcję to zbędna komplikacja - łatwo się pisze i programista może "się wykazać" jaki jest mądry i biegły w programowaniu ale... no właśnie jest poważne "ale". Żaden z języków programowania nie jest wieczny i po kilku latach robi się problem co dana funkcja biblioteczna robi i odtworzenie tego co się kryje w środku bywa pracochłonne - przerabiałem na sobie. Kod powstawał w latach 80 i 90 ubiegłego stulecia i Numpy ani pythona wtedy nie było. Wzorki się nie zmieniły a wszystkie te funkcje to nic innego jak najprostszy możliwy zapis tych wzorków. Zdecydowanie wolę podejście "patrzcie jakie to proste" niż, niestety częste obecnie "patrzcie jaki jestem mądry"... :-)

Edytowane 17 Czerwca 2022 przez Wiesiek1952

PiotrStopa · 17 Czerwca 2022

Myślę, że to zależy od zastosowania. Jeżeli potrzebujesz takie obliczenia od czasu do czasu to nie ma problemu. Jeżeli trzeba to liczyć kilka razy na sekundę to pewnie wspomaganie bibliotekami było wskazane (i pewnie przejście na coś innego niż Python).

Co do samego kodu ja bym się zastanowił co zwinąć do funkcji, co dokładnie to będzie wiadomo jak pomyśleć "co z tego kodu potrzebuje żeby wyliczyć położenie Jowisza". Na pierwszy rzut oka liczenie sum zwinąłbym do funkcji. To samo z przeliczenia mi między systemami koordynat.

zbuffer · 17 Czerwca 2022

@Wiesiek1952 Rozumiem co masz na myśli ale to zupełnie nietrafione w przypadku biblioteki, która jest wykorzystywana przez naprawdę wszystkich do robienia obliczeń. Idąc takim rozumowaniem programowanie obiektowe i biblioteki w ogóle nie powinny istnieć.

Kompletnie się nie zgadzam z podejściem, że wszystko trzeba pisać za pomącą czystego języka. Są pewne standardy i osoba początkująca musi je po prostu poznać, jeśli myśli o obliczeniach w Pythonie.

To tak jakbyś powiedział, że ktoś ma liczyć zawsze na kartce chociaż istnieją kalkulatory.

Trzeba najpierw nauczyć jakie są narzędzia i potem je wykorzystać.

Szczególnie, że jeśli dobrze wytłumaczysz te narzędzia albo odeślesz do dobrze przygotowanych materiałów na ich temat to taka osoba pozna coś dzięki czemu zrozumie nie tylko Twój kod, ale niezliczone przykłady, które znajdzie w internecie. Nie ma chyba bardziej typowego początku programu niż te dwie linie:

import numpy as np
import matplotlib.pyplot as plt

Edytowane 17 Czerwca 2022 przez zbuffer

szuu · 17 Czerwca 2022

1 godzinę temu, zbuffer napisał:

Idąc takim rozumowaniem programowanie obiektowe i biblioteki w ogóle nie powinny istnieć.

no właśnie, python to przecież również "black box" i on zaciemnia to co naprawde robi procesor

Wiesiek1952 · 17 Czerwca 2022

35 minut temu, szuu napisał:

no właśnie, python to przecież również "black box" i on zaciemnia to co naprawde robi procesor

"rence opadywują" - tego rodzaju przemyślenia są dla mnie za trudne.

oicam · 17 Czerwca 2022

Wiesiek, ja mam super pomysł.

Lepiej napisz taki program z własną biblioteką i z numpy i innymi i tu opublikuj. Chłopaki będą zachwyceni.

Będzie super edukacyjny:

import moja_biblioteka

PolicziZrobWszystko()

Panowie, facet robi fajną robotę. Jakie to ma znaczenie czy użyje numpy czy nie?

To są krótkie programy obliczeniowe dla astronomów amatorów. Albo argumenty, że python jest powolny. Proszę Was jakie to ma tutaj znaczenie? Żadne.

zbuffer · 17 Czerwca 2022

Użycie Numpy albo nie to drugorzędna sprawa. Dla mnie szkoda fajnej inicjatywy, która wg mnie jest kiepsko podana. Jako astro amator chciałbym zrozumieć teorię, która stoi za obliczeniami, na ile to możliwe i np zobaczyć wykres z orbitą czy z drogą Słońca na niebie itp.
Tutaj dostajemy kod z komentarzami lub nie, z mieszanką języków, i teorię wytłumaczoną w postach na forum albo w komentarzy do funkcji.

W przypadku tego typu materiału naprawdę widzę duży potencjał użycia Jupyter Notebooks, bo to najlepsze narzędzie do celów edukacyjnych. Mówię z doświadczenia bo uczę robotyki na uniwersytecie (w Hiszpanii) i tu prawie każdy przygotowuje materiały i ćwiczenia w Notebooks.

@Wiesiek1952 na jakiej licencji jest ten kod? Pytam bo w wolnym czasie mógłbym część przenieść do Colab i Ci pokazać co mam na myśli.

Wiesiek1952 · 17 Czerwca 2022

14 minut temu, zbuffer napisał:

@Wiesiek1952 na jakiej licencji jest ten kod? Pytam bo w wolnym czasie mógłbym część przenieść do Colab i Ci pokazać co mam na myśli.

Zupełnie nie myślałem nad jakimkolwiek "licencjonowaniem" tego oprogramowania więc najbliżej chyba do FLOSS no to niech i będzie FLOSS (co jest nawet ciut szerzej od "open source) - Free/Libre and Open Source Software.

zbuffer · 17 Czerwca 2022

Ok. Zawsze lepiej napisać żeby nie było żadnych problemów, bo teoretycznie jak coś wrzucasz bez licencji to nie wolno tego nikomu przesyłać dalej ani zmieniać.

zbuffer · 20 Czerwca 2022

Tutaj próbka tego co miałem na myśli. Wrzuciłem teorię, kody i dodałem kilka gadżetów. Możesz zobaczyć jak fajnie w dokumencie można osadzić równania i wykresy. Colab jest używany przez wiele osób pracujących nad deep-learning i big data, bo google udostępnia w tym celu GPU.

https://colab.research.google.com/drive/1jfG0PDYO3CoL4SINmCa-O2Aueiw7BlIB?usp=sharing

Jak otworzysz wszystko powinno być policzone ale możesz odpalić z menu "Runtime -> Run all".

Możesz dowolnie zmieniać kod, np testując sobie inne daty itp. bo ten link działa w opcji "viewer" - żadne zmiany nie zostaną zapisane w oryginalnym pliku. Jeśli chcesz zapisać zmiany musisz zapisać plik na swoim dysku google.

Edytowane 20 Czerwca 2022 przez zbuffer