Jump to content

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


Wiesiek1952
 Share

Recommended Posts

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 :matrix:

 

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:

 

 

 

 

image.png.6146cffa730834ed491b3d9d5bc532be.png

 

 

 

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/
  • Like 6
Link to comment
Share on other sites

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ć :D

 

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.

Link to comment
Share on other sites

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

  • Like 1
Link to comment
Share on other sites

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.

Link to comment
Share on other sites

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

Link to comment
Share on other sites

@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:emotion-5:

  • Like 1
Link to comment
Share on other sites

@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"...

  • Like 2
Link to comment
Share on other sites

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.

  • Like 1
Link to comment
Share on other sites

Posted (edited)

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

Edited by lkosz
Link to comment
Share on other sites

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

Link to comment
Share on other sites

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

Link to comment
Share on other sites

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 :D

 

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!

 

 

 

Link to comment
Share on other sites

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?

Link to comment
Share on other sites

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

 

  • Like 1
Link to comment
Share on other sites

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
      • 16 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.