Algorytmy astronomiczne czyli Python w rachunkach astro. - Cz. II czyli Słońce i stary dobry kalkulator HP 97

[Wiesiek1952]

Coś zupełnie spoza Meeusa i poważnej astronomii. Gdzieś w latach 90 modne były kalkulatory programowalne HP pracujące w notacji Łukasiewicza zwanej też polską lub odwrotną. Pomysł był super i szkoda, że go zarzucono choć nie do końca bo kompilatory czy interpretery dalej zmieniają nasz zapis w stylu x = (3 + 5.1) / 13.0 na zapis odwrotny czyli 3 ENTER 5.1 + 13 /      - nie są potrzebne nawiasy i zapis jest całkowicie jednoznaczny.   

 

Korzystając z notacji łukasiewicza całe (uproszczone) obliczenie efemeryd słońca dało się zmieścić w 256 krokach. Prawie niewiarygodne przy współczesnych "smart" i innych fonach.

 

Poniższy algorytm oblicza GHA słońca (Greenwich Hour Angle), który jest odpowiednikiem długości geograficznej na jakiej znajduje się słońce w danym momencie. Drugi parametr to deklinacja słońca czyli odpowiednik szerokości geograficznej nad jaką słońce się znajduje.

 

Dodałem kolejne klasy podobnie jak w cz I, kolejne funkcje ułatwiające życie jak signum i fractional (tak wiem, że można korzystać z gotowych ale moje liczą tak jak jest tu potrzebne) i są łatwo do analizy. Słynne signum zapisane w najzwięźlejszej postaci w C/C++ wygląda tak: return (x > 0) - (x < 0) prawda, że wygląda bardzo mądrze? I jest równie niezrozumiałe 

 

Dodałem też funkcje formatujące wydruk do strawnego dla oka zapisu oraz stałe typu TO_DEG co = 180.0 / pi. To dlatego, że te rachunki występują dziesiątki razy i łatwiej i szybciej jest wykorzystać stałą niż milion razy wykonywać działanie 180.0/pi

 

DOKŁADNOŚĆ: Algorytm wylicza GHA i Dec z dokładnością do 0.1-0.2' kątowej - wystarczająco dokładnie i w obserwacjach astro i w nawigacji (tej prawdziwej na statkach czy jachtach na oceanie).

 

W procedurze występuje zmienna GHA0 - tu niewykorzystywana choć można ją "wyjąć" minimalnie modyfikując kod. GHA0 to GHA punktu Barana (GHA Aries) czyli Greenwich Hour Angle punktu równonocy wiosennej, która jest odpowiednikiem południka zerowego dla współrzędnych ciał niebieskich. Na Ziemii pozycję określamy podając szerokość i długość geograficzną a na niebie Rektascencję (RA) i deklinację. Rektascencja liczona jest właśnie od punktu Barana i występuje w dwóch postaciach RA i SHA gdzie SHA to Sideral Hour Angle. Różnią się tym że jedno jest dopełnieniem do 360° drugiego. Sideral Time at Greenwich to właśnie GHA Aries.

 

Dla danych 2014 Oct 21, godz: 12:30:29 wychodzą takie wartości:

 

 

 

 

 

 

 

no i kod programu:

 

# This is a sample Python script.

# Press ⌃R to execute it or replace it with your code.
# Press Double ⇧ to search everywhere for classes, files, tool windows, actions, and settings.

import math

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

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


def signum(x: float) -> float:
    # THE BEST IN C + +
    # (x > 0) - (x < 0)
    # returning SIGNUM of double as double
    # return x > 0.0  ?  1.0: (x == 0.0  ?  0.0: -1.0);
    if x > 0:
        return 1.0
    elif x < 0:
        return -1.0
    else:
        return 0.0

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

def fraction(x: float) -> float:
    # returning FLOAT fractional part of float
    # --------------------------------------------

    return abs(x) - int(abs(x))

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


def string_latitude(fi: float) -> str:
    szer = " "
    znak = int(signum(fi))
    stop = abs(fi) * TO_DEG
    minuty = fraction(stop) * 60.0
    degrees = int(stop)
    if degrees < 10:
        szer += "0"
    szer += str(degrees) + "°"
    if minuty < 9.95:
        szer += "0"
    szer += f'{minuty:3.1f}'
    if znak > 0:
        szer += "\' N"
    if znak < 0:
        szer += "\' S"
    if znak == 0:
        szer += "   "
    return szer

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

def string_stopnie(rad: float) -> str:
    stopnie = ""
    znak = int(signum(rad))
    if znak < 0:
        stopnie += "-"
    stop = abs(rad) * TO_DEG
    minuty = fraction(stop) * 60.0
    deg = int(stop)
    if deg < 100:
        stopnie += "0"
    if deg < 10:
        stopnie += "0"
    stopnie += str(deg) + "°"
    if minuty < 10:
        stopnie += "0"
    stopnie += f'{minuty:3.1f}\''
    return stopnie


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



M_FACTOR = 10800.0 / math.pi  # STAŁA ZAMIANY Radianów na Minuty
REV_FACTOR = math.pi / 10800.0
TWO_PI = 2.0 * math.pi
TO_DEG = 180.0 / math.pi
TO_RAD = math.pi / 180.0


# ******************************************************************************************
#    Uproszczony sposob obliczania efemeryd slonca
#    Algorytm na podstawie programu na kalkulator HP - 97
#    Oryginalnie napisany w latach 80 XX w
# ******************************************************************************************

def sun_efemeries_HP97(current_date: day_in_calendar, current_time: present_time) -> sun_efemeries:
    GMT = current_time.hour + current_time.minute / 60.0 + current_time.seconds / 3600.0
    A = math.trunc(30.56 * current_date.month) + current_date.day

    c = A / 365.25 + current_date.year - 7
    b = 20 * c
    cc = 1
    if current_date.year / 4 == math.trunc(current_date.year / 4):
        cc = 2
    if current_date.month >= 3:
        cc = 3
    f = current_date.year / 4 - math.trunc(current_date.year / 4) + cc

    xx = math.sin(b * TO_RAD)
    g = (A - f + 4 * xx / 896.0 + 56.853795 + c / 128.0) * 360.0

    t = 15.0 * GMT
    r = (g + t) / 365.25
    v = 118.1 - c / 48.4 + r
    u = 0.2 * math.cos(v * TO_RAD) - 9.58
    w = 0.2 * math.sin(v * TO_RAD) * u + r
    fi = (math.cos(b * TO_RAD) + 8531.5 - b / 427.0) / 360.0
    w *= TO_RAD
    fi *= TO_RAD

    DECS = math.asin(-math.sin(fi) * math.sin(w))
    GHA0 = (r + t) * TO_RAD
    while GHA0 > TWO_PI:
        GHA0 = GHA0 - TWO_PI

    l1 = -math.sin(w) * math.cos(fi)
    m1 = math.sqrt(math.cos(w) * math.cos(w) + math.sin(w) * math.cos(fi) * math.sin(w) * math.cos(fi))
    sdf = l1 / m1
    GHA1 = math.asin(sdf)
    if math.cos(w) > 0:
        GHA1 = math.pi - GHA1
    GHAS = GHA0 - GHA1
    while GHAS > TWO_PI:
        GHAS = GHAS - TWO_PI
    while GHAS <= 0:
        GHAS = GHAS + TWO_PI

    efem = sun_efemeries()
    efem.sha = GHAS
    efem.declination = DECS
    efem.semi_diameter = 0

    return efem


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

# Press the green button in the gutter to run the script.
if __name__ == '__main__':
    day_a = day_in_calendar(2014, 10, 21)
    time_a = present_time(12, 30, 29)

    sun_data = sun_efemeries_HP97(day_a, time_a)

    print('SUN  - HP 97')
    print('--------------------------------')
    print('Sun  GHA     = ', string_stopnie(sun_data.sha))
    print('Sun  DEC     = ', string_latitude(sun_data.declination))
    print('--------------------------------')

# See PyCharm help at https://www.jetbrains.com/help/pycharm/

[lkosz]

39 minut temu, Wiesiek1952 napisał:

w C/C++ wygląda tak: return (x > 0) - (x < 0)

i jest to wadą C/C++. Widząc tak napisany kod mam ochotę strzelać

 

a w kwestii funkcji signum: numpy.sign() robi to, czego potrzebujesz. Efemerydy policzysz pakietem astropy. Próg wejścia w niego jest ciut wysoki, ale nie jest wyższy niż pisanie na piechotę. Generalnie jeśli jest na coś gotowiec który jest standardem, to dobrą praktyką jest użycie gotowca, nawet jeśli się to wydaje mało wygodne. Potem przy aktualizacji się to opłaca, bo nie ma się na głowie jednorazowych adminów/jednorazowych deweloperów. Dodatkowo jest taki bonus, że ktoś za nas wykonuje optymalizację kodu.

[Wiesiek1952]

6 minut temu, lkosz napisał:

... to dobrą praktyką jest użycie gotowca, nawet jeśli się to wydaje mało wygodne...

 

Tylko ma to zerową wartość edukacyjną no i jak to się ma do powyższego algorytmu z HP-97, którego zresztą nigdzie nie trafiłem w necie a po prostu miałem kiedyś HP-97 i na pasku magnetycznym był programik liczący efemerydy słońca w taki właśnie sposób. Wystarczyło poświęcić trochę czasu i zrobić "reverse engineering" i dzięki temu mam funkcję, którą mogę użyć np. w excelu do liczenia fotowoltaiki bez nadmiernych komplikacji. Ale tak, gdybym pisał program użytkowy z całym środowiskiem i poważnym UI to pewnie sięgnął bym do gotowca. W czasach kiedy to pisałem gotowców nie było, pythona zresztą też nie było...

[lkosz]

Wartość edukacyjna jest następująca: nie oszukiwać się, że się samemu zrobi coś lepiej, tylko użyć standardu. Po to są standardowe biblioteki, żeby ich używać, żeby inni potem umieli skorzystać z takiego kodu i nie musieli pół godziny główkować nad wielką funkcją pozbawioną dokumentacji szukając w niej błędu. Jednorazowe programowanie to dotkliwa pułapka. Kończy się porzuconym, nieczytelnym, nieoptymalnym kodem, którego nikomu nie opłaca się poprawiać. Jak przychodzi do aktualizacji, to się robi problem, w efekcie i tak trzeba jeszcze raz napisać na prędce to samo, tylko wówczas opóźnia to resztę zadań. Jedna osoba nie jest w stanie konkurować z zespołem pracującym nad biblioteką. Jak się ma jeden skrypcik w domu, to można nie poczuć skutków, ale jak się ma tysiąc maszyn, z tysiącem takich skrypcików robionych przez tysiąc osób, każda po swojemu nie trzymając się standardu, to się z tego robi horror.

 

Twój kod można mocno zoptymalizować, np. te działania modulo... oraz skrócić do postaci zrozumiałej od razu na pierwszy rzut oka.

[Wiesiek1952]

12 minut temu, lkosz napisał:

...

Twój kod można mocno zoptymalizować, np. te działania modulo... oraz skrócić do postaci zrozumiałej od razu na pierwszy rzut oka.

No przecież pisałem o tym w cz I. Wszystkie Twoje programy są oczywiście optymalne i nie można ich napisać lepiej. łechhhhhh....

 

dalej nie widzę związku z algorytmem z HP-97

[oicam]

@Wiesiek1952 wspaniały temat i bardzo edukacyjny. Gotowce są fajne, ale to czarna magiczna skrzynka.

 

W archiwach Sky&Telescope jest jeszcze pełno źródeł napisanych w Basic-u. Można by je było przetłumaczyć na Pythona.

 

Fajnie by było napisać kod do przewidywania tranzytów WCP

[astro4u.net]

@Wiesiek1952zainspirowałeś mnie, fajny projekt. Po tym tekście wyszukałem kalkulator, który leżał od studiów na półce to HP 42S z odwrotną notacją polską i zastanawiam się jak go można teraz użyć.

[oicam]

@lkosz przy projektach większych niż swój własny programik do przeliczania dat, albo wyznaczania momentów górowania Słońca czy innych  podobnych to masz rację. @Wiesiek1952 pokazuje krótkie edukacyjne programy dla tych, którzy chcieliby sami popróbować cokolwiek napisać, a przy okazji liznąć trochę astronomii obliczeniowej.

Jakby miał tutaj korzystać z gotowców, to w zasadzie mógłby dać link do githuba i dopisać "sami sobie sprawdźcie jak to liczy"...

[Wiesiek1952]

9 minut temu, astro4u.net napisał:

@Wiesiek1952zainspirowałeś mnie, fajny projekt. Po tym tekście wyszukałem kalkulator, który leżał od studiów na półce to HP 42S z odwrotną notacją polską i zastanawiam się jak go można teraz użyć.

Ta odwrotna notacja polska to wymysł amerykanów i ichnich "Polish jokes" według których Polacy wszystko robią odwrotnie :-) Notacja Lukasiewicza jakoś lepiej brzmi :-) W samolotach jest usterzenie V lub odwrócone V (coś jak A bez kreski) - u nas to usterzenie Rudlickiego , który to chyba pierwszy zastosował. :-)

 

Nie znam HP 42s ale one wszystkie podobne i bardzo szkoda, że umarły i zniknęły. W Apollo 13 wszystkie rachunki związane z wyliczeniem trajektorii powrotu były liczone na takich kalkulatorach - kurde nie pamiętam HP-65 albo coś w tym stylu.

[astro4u.net]

Ten mój to taki wynalazek, nawet wykopałem instrukcje

[lkosz]

@oicam Podstawowy problem w obliczeniach koordynatów to bałagan w jednostkach astronomicznych. Pamiętam jak sam przez to brnąłem 2 lata temu. Jak się TO załapie, to reszta jest nudnymi rachunkami na podstawie wzorów z książki, którą właśnie zajmują się biblioteki. Niestety żeby ich użyć, trzeba wiedzieć co się chce policzyć. Rozumiem entuzjazm, ale ten skrypt nie jest szczególnie dobrą demonstracją prostoty Pythona przeznaczoną dla nowicjuszy. Nie z klasami i obliczeniami robionymi bardzo naokoło. Można o wiele prościej, w kilku linijkach, składnią zrozumiałą na pierwszy rzut oka. Nie jest to wycieczka personalna, ot po prostu polemika z kodem.

Edytowane 2 Czerwca 2022 przez lkosz

[oicam]

@lkosz Każdy piszący swój kod ma swoje kung fu. Ja na ten przykład do obliczeń w ogóle nie wchodzę w programowanie obiektowe bo moim zdaniem to zbędny zabieg. Ale jak ktoś lubi.

[Wiesiek1952]

1 minutę temu, lkosz napisał:

@oicam Podstawowy problem w obliczeniach koordynatów to bałagan w jednostkach astronomicznych. Pamiętam jak sam przez to brnąłem 2 lata temu. Jak się TO załapie, to reszta jest nudnymi rachunkami na podstawie wzorów z książki, którą właśnie zajmują się biblioteki. Rozumiem entuzjazm, ale ten skrypt nie jest szczególnie dobrą demonstracją prostoty Pythona przeznaczoną dla nowicjuszy. Nie z klasami i obliczeniami robionymi bardzo naokoło. Można o wiele prościej, w kilku linijkach, składnią zrozumiałą na pierwszy rzut oka. Nie jest to wycieczka personalna, ot po prostu polemika z kodem.

 

 

Dalej czekam na wyjaśnienie związku z algorytmem HP-97.  

 

Tak poza tym - ostatnio w nauczaniu klasy wprowadza się bardzo wcześnie. W Pythonie typy podstawowe to też klasy. W tym moim kodzie klasa jest odpowiednikiem struktury czy kontenera danych, nie ma metod czy innych zaawansowanych rzeczy. Jest używana niejako "zamiast" niemodyfikowalnych tuples (tfu, po polsku krotka, tfu...)

[szuu]

Godzinę temu, Wiesiek1952 napisał:

Nie znam HP 42s ale one wszystkie podobne i bardzo szkoda, że umarły i zniknęły.

umarły i zniknęły bo nie mają sensu w swojej dawnej roli, zastąpiły je nowe lepsze narzędzia.

ale jednak ciągle żyją w sercach fanów

 

https://thomasokken.com/free42/ - softwarowy klon hp-42s

 

https://www.swissmicros.com/products - i są nawet do kupienia hardwarowe implementacje zachowujące ducha dawnych kalkulatorów ale już we współczesnej technologii!

 

 

 

[tunarob]

10 godzin temu, Wiesiek1952 napisał:

W tym moim kodzie klasa jest odpowiednikiem struktury czy kontenera danych, nie ma metod czy innych zaawansowanych rzeczy. Jest używana niejako "zamiast" niemodyfikowalnych tuples (tfu, po polsku krotka, tfu...)

To może warto użyć dataclass? https://docs.python.org/3/library/dataclasses.html lub namedtuple?

[smopi]

13 godzin temu, oicam napisał:

Gotowce są fajne, ale to czarna magiczna skrzynka.

Nie, jeśli jest to Open Source.

[oicam]

20 minut temu, smopi napisał:

Nie, jeśli jest to Open Source.

Nie chodzi mi oto, że nie można podejrzeć jak działa metoda, tylko program napisany przez Kolegę Wieśka składałby się z trzech linijek kodu.

I gdzie tu edukacja? Każdy kto chce teraz zacząć robić cokolwiek w pythonie, zerknie na ten wątek i może się zainspiruje.

Pomyśli "aaaa, fajnie, może sobie napiszę krótki programik, żeby mi przeliczał daty, albo przewidzi tranzyty księżyców galileuszowych"...

I sam powoli, małymi krokami napisze sobie program, zacznie szukać bibliotek itd, itd...