Postępy Astronomii
Tom XXXVI (1988), Zeszyt 1, 49–55.
Katedra Radioastronomii
Uniwersytetu M. Kopernika (Toruń)
(Otrzymano 22 października 1987 r.)
| S t r e s z c z e n i e — W pracy przedstawiono zestaw programów w języku FORTRAN dla komputerów osobistych pozwalających na obliczenie pozycji Słońca, Księżyca i planet z dokładnością 1' na dowolny moment w okresie 600 lat ześrodkowanym na latach 70. naszego stulecia. |
| SOLAR SYSTEM EPHEMERIDES FOR EVERYBODY. S u m m a r y - A set of FORTRAN routines suitable for personal computers of the IBM PC type that allow to find geocentric coordinates of the Sun, Moon and all planets with 1' accuracy within about 300 years of the present and based on Van Flandern and Pulkkinen (1979) work is described and examples of printouts are appended. |
|
|
Zjawiska związane z ruchem ciał niebieskich, Słońca, Księżyca i planet, zawsze wzbudzały żywe i powszechne zainteresowanie ludzi, nie tylko tych związanych zawodowo z astronomią. Od zarania dziejów próbowano rozszyfrować prawa nimi rządzące. I chociaż owe prawa ostatecznie okazały się względnie proste (prawa Keplera), to do dziś sztuka przewidywania wzajemnych położeń ciał niebieskich dostępna jest w zasadzie tylko wąskiej grupie specjalistów. Wynika to z faktu mnogości ciał wzajemnie na siebie oddziaływających grawitacyjnie, jak również z istnienia efektów relatywistycznych i innych niegrawitacyjnych. Czynniki te, chociaż w większości przypadków, albo nawet zgrupowane, zdają się być mało znaczące w danym momencie, kumulują się z upływem czasu do trudnych do oceny rozmiarów.
Ktoś zainteresowany efemerydami, czyli położeniami ciał, na wybraną chwilę w przeszłości lub w niedalekiej przyszłości, jest najprawdopodobniej skazany na odszukanie odpowiedniego rocznika astronomicznego, odczytanie zeń potrzebnych informacji i przeprowadzenie pewnych dodatkowych rachunków. Istnieje przy tym poważne niebezpieczeństwo niewłaściwej interpretacji odczytu albo pomyłek przy elementarnych redukcjach na daną chwilę i położenie geograficzne, dlatego w praktyce zawsze należałoby korzystać dodatkowo z pomocy zawodowego astronoma. Nie jest to z pewnością sytuacja zadowalająca. Jest jednak jeszcze gorzej jeśli w grę wchodzą epoki odległe o stulecia lub więcej od współczesności i wtedy można dopuścić, że zainteresowana osoba pod naporem trudności w ogóle zrezygnuje z poszukiwań.
Tak było od dawien i dziś jeszcze. Sytuacja zdaje się jednak radykalnie zmieniać na naszych oczach z powodu niebywałego rozwoju technik obliczeniowych. Wprawdzie na powszechną dostępność bardzo precyzyjnych efemeryd raczej nie można jeszcze liczyć, ale te najczęściej zadowalające przeciętnego użytkownika są już w zasięgu ręki, przynajmniej dla tej dość już licznej grupy, która ma dostęp do komputerów osobistych. Jedna typowa dyskietka z zapasem pomieści programy umożliwiające obliczanie od ręki i niezawodnie przeróżnych zjawisk związanych z ruchami ciał niebieskich widzianych z dowolnego miejsca na Ziemi i w dowolnej chwili z zakresu obejmującego całe stulecia. Podstawowe programy takiego zestawu są przedmiotem tego opisu.
Algorytmy obliczania położeń ciał Układu Słonecznego, które przedstawię zapewniają formalnie dokładność 1' na przestrzeni około 600 lat ześrodkowanej na latach 70. naszego stulecia (z wyjątkiem Plutona, dla którego niepewność położenia ocenia się na ok. 15'). Kilka lat temu zakodowałem je w języku FORTRAN IV na komputerze RIAD 32, a następnie (1986/87) zaadoptowałem na mniejszej jednostce, SM1420 (odpowiednik PDP 11/35). Wkrótce też pojawiła się możliwość przeniesienia programów na komputer osobisty typu IBM PC. Właśnie ten ostatni etap wraz z faktem, że algorytmy sprawdziły się bardzo dobrze w czasie kilkuletniej eksploatacji sprawił, że zdecydowałem się spopularyzować je w naszej społeczności astronomicznej.
Z praktycznych dotąd zrealizowanych zastosowań wymienię tutaj kalendarz astronomiczny zawierający m.in. codzienne współrzędne Słońca, Księżyca i planet i ich momenty wschodów i zachodów, momenty pór roku i faz Księżyca, program na znajdowanie momentów i obliczanie okoliczności zaćmień Słońca i Księżyca oraz horoskopy (bez interpretacji, naturalnie). Zastosowania te wymagają oczywiście dodatkowych algorytmów, których teraz nie będę opisywał w szczegółach, ale zainteresowany Czytelnik może niewielkim nakładem pracy sam je stworzyć bazując na przedstawionych tu efemerydach i własnej inwencji. Np., wyznaczenie momentu najbliższego nowiu Księżyca może być zrealizowane przez iteracyjny proces, w którym elongacja Księżyca (róznica długości ekliptycznej Księżyca i Słońca) zbiega do zera, a inkrementy czasu dobierane są w zależności od znaku i szybkości zmian elongacji (metoda Newtona znajdowania zera funkcji, którą tutaj jest elongacja). W tym przykładzie 3–5 iteracji da dokładność rzędu 1 minuty czasu, co jest na granicy dokładności efemerydy Księżyca (nasz naturalny satelita w ciągu minuty przemieszcza się po niebie o ok. 0,5').
Z naszego doświadczenia wynika, że efemerydy te są w istocie dokładniejsze niż wspomniana 1', chociaż nie możemy tego stwierdzenia poprzeć jakąkolwiek analizą statystyczną, ze względu na wyrywkowość przeprowadzonych obliczeń i porównań. Pewne potwierdzenie wspomnianej dokładności i wartości tych algorytmów zdaje się wynikać także z prac Todorovic-Juchniewicz (1985) publikowanych na tych łamach.
Pełny opis algorytmów przedstawili ich autorzy: Van Flandern i Pulkkinen (1979) (w dalszym ciągu będziemy używali skrótu VFP). Oczywiście, nie będziemy tutaj powtarzać całości tej publikacji, której lwią część stanowią tabele współczynników liczbowych do wyrażeń trygonometrycznych na obliczanie współrzędnych heliocentrycznych i (osobno) równikowych. Ponieważ sposób prowadzący do położeń heliocentrycznych wydaje się nieco prostszy, ograniczymy się tylko do niego pomijając całkowicie ten, który daje wprost interesujące nas geocentryczne współrzędne równikowe. Wprawdzie w VFP podano w zasadzie pełną informację potrzebną przy praktycznym zastosowaniu algorytmów, ale zdając sobie sprawę z trudności, jakie czekają nawet wprawnego w programowaniu czytelnika, któremu jednak mogą być obce meandry analitycznych teorii ruchu ciał Układu Słonecznego, podamy tu wystarczająco dużo materiału — tak, by nawet nie znający języka angielskiego mógł sam zaprogramować całość efemeryd posługując się tym opracowaniem i oryginałem VFP.
Podstawowym argumentem prezentowanej teorii VFP jest czas, t, liczony w dniach od południa 1 stycznia 2000 r. Zmienną tę można wyrazić wzorem
| t = JD – 2451545.5 + UT/24, | (1) |
Algorytmy VFP stanowią w istocie pewne uproszczenie dokładnych analitycznych teorii. Kilka podstawowych informacji o konstrukcji takich teorii zawiera wspomniana już praca Todorovic-Juchniewicz (1985). W naszym przypadku każdą ze współrzędnych planety oblicza się z wyrażenia typu:
| w(t) = wo(t) + Σk aksin[Σi Jik Aik(t)], | (2) |
Tak więc trzonem algorytmów VFP są tabele współczynników a i J dla każdej planety i każdej ze współrzędnych oraz zestaw elementów średnich A w postaci prostych, liniowych funkcji czasu. Mamy nadzieję, że to skrótowe objaśnienie wystarczy zainteresowanemu Czytelnikowi przy jednoczesnej inspekcji podanych w Dodatku konkretnych przykładów.
W naszych efemerydach szczególną rolę odgrywa Słońce, gdyż jego współrzędne (a właściwie współrzędne Ziemi) stanowią punkt odniesienia dla obliczeń geocentrycznych położeń wszystkich planet. Drugim wyjątkiem jest Księżyc, którego pozycję dostaje się wprost w układzie geocentrycznym i z odległością (od Ziemi) wyrażoną nie w jednostkach astronomicznych, jak to jest u wszystkich pozostałych obiektów, lecz w promieniach Ziemi (1 j.a. = 23454,8 promieni Ziemi). Te fakty uzasadniają pełne przedstawienie podprogramów na obliczanie pozycji tych dwóch ciał (w Dodatku procedury SLONCE i KSIEZY). Uzupełnia to reprezentatywna procedura dla planet (WENUS). Wybraliśmy w to miejsce Wenus dlatego, że spośród wszystkich ośmiu planet odpowiadają jej najprostsze wzory.
W celu zakodowania algorytmu VFP dla dowolnej z pozostałych planet wystarczy powtórzyć całą strukturę podprogramu WENUS zmieniając jedynie jego nazwę, parametry A (u VFP są one oznaczone m.in. numerami od 1 do 33 zgodnie ze wskaźnikami przy A użytymi w mojej implementacji) i wyrażenia arytmetyczne na obliczanie heliocentrycznej ekliptycznej długości (w załączonym programie DLH) i szerokości (SZH) oraz odległości od Słońca (RH) w oparciu o stabelaryzowane dane w VFP. Czytelnik może nieco ułatwić sobie pracę przez utworzenie w swoim programie tablic współczynników A i J i obliczanie wyrażeń typu (2) w pętlach DO. Warto jednak wziąć pod uwagę to, że większość ze współczynników J ma wartość zero, dlatego mój jawny zapis sum będzie z pewnością bardziej efektywny w liczeniu (nie zawiera on mnożeń przez 0 ani 1). W bloku COMMON umieściłem te ze współczynników A, które występują jednocześnie w podprogramie SLONCE i w podprogramach na kilka planet (w których nie ma zatem potrzeby ponownego ich oblicznia). Zmienne SE i CE oraz DL w drugim bloku /SUN/ są sinusem i kosinusem kąta nachylenia równika do ekliptyki oraz poprawką w długości ekliptycznej związaną z nutacją — wszystkie obliczone dla epoki t. Oznacza to, że otrzymane współrzędne ciał niebieskich będą odniesione do ekliptyki i równika daty (precesję uwzględniono niejawnie w wyrażeniach typu (2) dla długości ekliptycznych).
Zestaw podprogramów na współrzędne opisanych trzech ciał zamyka procedura zamiany współrzędnych heliocentrycznych planet na geocentryczne (HELGEO). Algorytm zakodowany w tym podprogramie w części obliczania współrzędnych równikowych (RA, DEC i R), chociaż zapisany inaczej, jest oczywiście równoważny wzorom (8) i (9) w VFP; reszta jest moim dodatkiem, podobnie jak wzory na równanie czasu (EQT) i czas gwiazdowy w Greenwich (GST), które włączyłem do procedury SLONCE.
Całość Dodatku obejmuje ponadto prosty program demonstracyjny (PRZYKLAD), który ilustruje sposób korzystania z opisanych procedur. Oto przykład działania tego programu dla momentu czasu, dla którego w VFP podano również testowe obliczenia:
DATA: 1969 6 28, UT: 0.00
REKT. DEKL. PROM. PARAL./ODL. DLE SZE
[h] [stop.] ['] [rad]
SLONCE: 6.4449 23.3014 15.7574 8.6500" 1.67774
KSIEZYC: 16.4983 -26.7378 16.5635 60.7880' 4.36190 -0.08468
WENUS: 3.2754 15.0979 0.7864ja 0.88591 -0.05152
DJ = 2440400.50000, GST = 18.3946 h, EQT = -3.021 min.
|
Na koniec chcę jeszcze ostrzec, że niektóre kompilatory zbuntują się na zbyt złożone wyrażenia arytmetyczne w procedurze KSIEZY, i wtedy trzeba będzie je rozbić na kilka wyrażeń (sam używałem bezproblemowo kompilatora Lahey).
W czasie przygotowywania tej pracy korzystałem z pomocy finansowej w ramach problemu RPBR Nr RR.I.11/2.
LITERATURA
Borkowski K.M., 1987, ,,Dni juliańskie i daty
kalendarzowe", Post. Astr. 35, 275–279.
Todorovic-Juchniewicz B., 1985, Post. Astr. 33, 59 i 153.
Van Flandern T.C., Pulkkinen K.F., 1979, Astrophys. J. Supp. Ser.
41, 391.
DODATEK
PROGRAM PRZYKLAD
REAL*8 DJ,DG
JD(L,M,N) = (L+(M-8)/6)*1461/4+(153*MOD(M+9,12)+2)/5+N+1721119
* -(L+100+(M-8)/6)/100*3/4
DG=57.295779513083D0
1 WRITE(*,'(/A)') ' Rok, miesiac, dzien, UT[h.mm] ? '
READ(*,*) IY,M,ID,UT
DJ=JD(IY,M,ID) - .5D0 + (AINT(UT) + AMOD(UT,1.)/.6)/24D0
CALL SLONCE(DJ,RAS,DS,DLS, RS,EQT,GST)
CALL KSIEZY(DJ,RAK,DK,DLK,SK,RK)
CALL WENUS ( RAW,DW,DLW,SW,RW)
WRITE(*,2) IY,M,ID,UT,RAS*DG/15.,DS*DG,16.01967/RS,8.794/RS,DLS
*, RAK*DG/15.,DK*DG,ASIN(.272493/RK)*DG*60.,ASIN(1./RK)*DG*60.,DLK
*, SK,RAW*DG/15.,DW*DG,RW,DLW,SW,DJ,GST*DG/15.,EQT*DG*4.
2 FORMAT(1H /15X,'DATA:',I5,2I3,', UT:',F6.2/15X,'REKT. DEKL. '
*,' PROM. PARAL./ODL. DLE SZE'/16X,'[h] [stop.] ['']'
*,20X,'[rad]'/3X,'SLONCE: ',4F9.4,'"',F10.5/3X,'KSIEZYC: ',4F9.4,
*''' ',2F9.5/3X,'WENUS:',3X,2F9.4,9X,F9.4,'ja',2F9.5/6X,
*'DJ =',F14.5,', GST =',F8.4,' h, EQT =',F7.3,' min.'/)
GO TO 1
END |
SUBROUTINE SLONCE(DJ,RA,DEC,DLE,R,EQT,GST)
C
C PROCEDURA OBLICZA WSPOLRZEDNE [W RAD] ROWNIKOWE (RA I DEC), DLUGOSC
C EKLIPTYCZNA (DLE), ODLEGLOSC SLONCA OD ZIEMI (R) [W JEDN. ASTR.] ORAZ
C [W RAD] ROWNANIE CZASU (EQT) I CZAS GWIAZDOWY GREENWICH (GST) NA MO-
C MENT OKRESLONY PRZEZ DJ, TJ. DZIEN JULIANSKI.
C COMMON BLOCK /SUN/ I /PLAN/ ZAWIERAJA DANE DO OBLICZEN POZYCJI PLANET.
C
COMMON /SUN/ XI0,ETA0,SE,CE,DL /PLAN/ T,A1,A7,A8,A13,A16,A19
DOUBLE PRECISION DJ,D1,D2,T
F(D1,D2) = 6.283185307179D0*DMOD(D1 + D2*T,1D0)
T = DJ-2451545d0
A1 = F(.606434D0,.03660110129D0)
A5 = F(.347343D0,-1.4709391D-4)
A7 = F(.779072D0,.273790931D-2)
A8 = F(.993126D0,.27377785D-2)
A13 = F(.140023D0,4.45036173D-3)
A16 = F(.053856D0,.145561327D-2)
A19 = F(.056531D0,.23080893D-3)
DL = -8.34E-5*SIN(A5)
DLE = A7+DL+((6893.-4.6543D-4*SNGL(T))*SIN(A8)+72.*SIN(A8+A8)-
*7.*COS(A8-A19)+6.*SIN(A1-A7)+5.*SIN(4.*A8-8.*A16+3.*A19)-
*5.*COS(2.*(A8-A13))-4.*SIN(A8-A13)+4.*COS(4.*(A8-2.*A16)+3.*A19)+
*3.*(SIN(2.*(A8-A13))-SIN(A19)-SIN(2.*(A8-A19))))/206264.8
IF(DLE.LT..0) DLE = DLE+6.2831853
R = 1.00014-.01675*COS(A8)-.00014*COS(A8+A8)
SE = SIN(.40905013 - 6.214E-9*SNGL(T) + 4.36E-5*COS(A5))
CE = SQRT(1. - SE*SE)
XI0 = R*COS(DLE)
ETA0= R*SIN(DLE)
RA = ATAN2(CE*ETA0,XI0)
DEC = ASIN(SE*ETA0/R)
EQT = AMOD(A7 - RA,6.283185)
IF(ABS(EQT).GT.3.) EQT = EQT - SIGN(6.2831853,EQT)
IF(RA.LT..0) RA = RA + 6.2831853
GST = DMOD(T*6.283185307179D0+A7+DL,6.283185307179D0)
IF(GST.LT.0.) GST = GST + 6.2831853
RETURN
END |
SUBROUTINE KSIEZY(DJ,RA,DEC,DLE,SZE,R)
C
C PROCEDURA OBLICZA WSPOLRZEDNE [W RAD] EKLIPTYCZNE (DLE I SZE) I
C ROWNIKOWE (RA I DEC) ORAZ ODLEGLOSC KSIEZYCA OD ZIEMI [W PROMIE-
C NIACH ZIEMI] NA MOMENT OKRESLONY PRZEZ DZIEN JULIANSKI (DJ).
C
DOUBLE PRECISION DJ,D1,D2,T
F(D1,D2) = 6.28318530718D0*DMOD(D1 + D2*T,1D0)
T = DJ-2451545D0
A1 = F(.606434D0,.03660110129D0)
A2 = F(.374897D0,.03629164709D0)
A3 = F(.259091D0,.0367481952D0)
A7 = F(.779072D0,.273790931D-2)
A4 = A1-A7
A8 = F(.993126D0,.27377785D-2)
A12= F(.505498D0,4.45046867D-3)
DLE= A1-8.34E-5*SIN(A1-A3)+(22640.*SIN(A2)-4586.*SIN(A2-A4-A4)+
* 2370.*SIN(A4+A4)+769.*SIN(A2+A2)-668.*SIN(A8)-125.*SIN(A4)-
* 412.*SIN(A3+A3)-212.*SIN(2.*(A2-A4))-206.*SIN(A2-A4-A4+A8)+
* 192.*SIN(A2+A4+A4)+165.*SIN(A4+A4-A8)+148.*SIN(A2-A8)-
* 110.*SIN(A2+A8)-55.*SIN(2.*(A3-A4))-45.*SIN(A2+A3+A3)+40.*
* SIN(A2-A3-A3)-38.*SIN(A2-4.*A4)+36.*SIN(3.*A2)-31.*SIN(2.*(A2-
* A4-A4))+28.*SIN(A2-A4-A4-A8)-24.*SIN(A4+A4+A8)+19.*SIN(A2-A4)+
* 18.*SIN(A4+A8)+15.*SIN(A2+A4+A4-A8)+14.*(SIN(2.*(A2+A4))+
* SIN(4.*A4))-13.*SIN(3.*A2-A4-A4)-11.*SIN(A2+16.*A7-18.*A12)+
* 10.*SIN(A2+A2-A8)+9.*(SIN(A2-2.*(A3+A4))+COS(A2+16.*A7-18.*A12)-
* SIN(2.*(A2-A4)+A8))-8.*SIN(A2+A4)+8.*(SIN(2.*(A4-A8))-SIN(2.*A2-
* (-A8)))-7.*(SIN(A8+A8)+SIN(A2+2.*(A8-A4))-SIN(A1-A3))-
* 6.*(SIN(A2+2.*(A4-A3))+SIN(2.*(A3+A4)))-4.*(SIN(A2-4.*A4+A8)+
* SIN(2.*(A2+A3)))+3.*(SIN(A2-3.*A4)-SIN(A2+2.*A4+A8)-
* SIN(2.*A2-4.*A4+A8)+SIN(A2-A8-A8)+SIN(A2-2.*(A8+A4)))+
* 2.*(SIN(A2+A2-A4)+SIN(4.*A4-A8)+SIN(4.*A2)+SIN(A2+4.*A4)-
* SIN(2.*(A3-A4)+A8)-SIN(2.*(A2-A4)-A8))+5.6569*
* (SNGL(T)/36525.+1.)*SIN(A2+16.*A7-18.*A12+.785398))/206264.8
IF(DLE.LT..0) DLE = DLE+6.2831853
SZE = (18461.*SIN(A3)+1010.*SIN(A2+A3)+1000.*SIN(A2-A3)-
* 624.*SIN(A3-A4-A4)-199.*SIN(A2-A3-A4-A4)-167.*SIN(A2+A3-A4-A4)+
* 117.*SIN(A3+A4+A4)+62.*SIN(A2+A2+A3)+33.*SIN(A2-A3+A4+A4)+
* 32.*SIN(A2+A2-A3)-30.*SIN(A3-A4-A4+A8)-16.*SIN(A2+A2+A3-A4-A4)+
* 15.*SIN(A2+A3+2.*A4)+12.*SIN(A3-A4-A4-A8)-8.*SIN(A1)-
* 9.*SIN(A2-A3-A4-A4+A8)+8.*SIN(A3+A4+A4-A8)+7.*(SIN(A2+A3-A8)-
* SIN(A2+A3-2.*A4+A8)-SIN(A2+A3-4.*A4))-6.*(SIN(A3+A8)+SIN(3.*A3)-
* SIN(A2-A3-A8))+5.*(SIN(A3-A4)+SIN(A3-A8)-SIN(A3+A4)-
* SIN(A2+A3+A8)-SIN(A2-A3+A8))+4.*(SIN(3.*A2+A3)-SIN(A3-4.*A4))+
* 3.*(SIN(A2-3.*A3)-SIN(A2-A3-4.*A4))+2.*(SIN(3.*A2-A3)+
* SIN(A2+A2-A3-A4-A4)+SIN(A2-A3+A4+A4-A8)-SIN(A2+A2-A3-4.*A4)+
* SIN(A2+A2-A3+A4+A4)-SIN(3.*A3-A4-A4)))/206264.8
R = 60.36298-3.27746*COS(A2)-.57994*COS(A2-A4-A4)-
* .46357*COS(A4+A4)-.08904*COS(A2+A2)+.03865*COS(2.*(A2-A4))-
* .03237*COS(A4+A4-A8)-.02688*COS(A2+A4+A4)-.02358*COS(A2-A4-A4+
* A8)+.01247*COS(A2-A3-A3)+1E-5*(704.*COS(A8)+529.*COS(A4+A4+A8)-
*2030.*COS(A2-A8)+1719*COS(A4)+1671.*COS(A2+A8)-524.*COS(A2-4.*A4)
* +398.*COS(A2-A4-A4-A8)-366.*COS(3.*A2)-295.*COS(A2+A2-4.*A4)-
* 263.*COS(A4+A8)+249.*COS(3.*A2-A4-A4)-221.*COS(A2+A4+A4-A8)+
* 185.*COS(2.*(A3-A4))-161.*COS(2.*(A4-A8))+
* 147.*COS(A2-2.*(A3-A4))-142.*COS(4.*A4)+139.*COS(2.*(A2-A4)+A8)-
* 118.*COS(A2-4.*A4+A8)-116.*COS(2.*(A2+A4))-110.*COS(A2+A2-A8))
SE = SIN(.40905013 - 6.214E-9*SNGL(T) + 4.36E-5*COS(A1-A3))
CE = SQRT(1. - SE*SE)
SD = SIN(DLE)
DEC= ASIN(SD*SE*COS(SZE) + SIN(SZE)*CE)
RA = AMOD(ATAN2(SD*CE-TAN(SZE)*SE,COS(DLE))+6.2831853,6.2831853)
RETURN
END |
SUBROUTINE WENUS(RA,DEC,DLE,SZE,R)
C
C PROCEDURA OBLICZA WSPOLRZEDNE ROWNIKOWE (RA I DEC) I EKLIPTYCZNE
C (DLE I SZE) [RAD] ORAZ ODLEGLOSC OD ZIEMI [JEDN. ASTR.] WENUS NA
C MOMENT ODLEGLY O T DNI OD POLUDNIA DNIA 1 STYCZNIA 2000 ROKU.
C WYMAGA ONA UPRZEDNIEGO WYWOLANIA PODPROGRAMU SLONCE.
C
COMMON /PLAN/ T,A1,A7,A8,A13,A16,A19
DOUBLE PRECISION T,D1,D2
F(D1,D2) = 6.28318530718D0*DMOD(D1 + T*D2,1D0)
A12 = F(.505498D0,4.45046867D-3)
A14 = F(.292498D0,4.45040017D-3)
DLH = A12+(SIN(A13)*(2794.-SNGL(T)/1826.25)-181.*SIN(A14+A14)
*+12.*SIN(A13+A13) - 10.*COS(2.*(A8-A13)) + 7.*COS(3.*(A8-A13)))
* /206264.8
SZH = (12215.*SIN(A14) + 166.*SIN(A13)*COS(A14))/206264.8
RH = .72335 - .00493*COS(A13)
CALL HELGEO(DLH,SZH,RH,RA,DEC,DLE,SZE,R)
RETURN
END |
SUBROUTINE HELGEO(DLH,SZH,RH,RA,DEC,DLE,SZE,R) C C PODPROGRAM TEN PRZELICZA HELIOCENTRYCZNE WSPOLRZEDNE PLANET (DLH, C SZH I RH) NA WSPOLRZEDNE GEOCENTRYCZNE (RA, DEC, DLE, SZE I R) C COMMON /SUN/ XI0,ETA0,SE,CE,DL C = RH*COS(SZH) X = C*COS(DLH + DL) + XI0 C = C*SIN(DLH + DL) + ETA0 B = RH*SIN(SZH) Y = C*CE - B*SE Z = C*SE + B*CE R = SQRT(X*X + Y*Y + Z*Z) RA = AMOD(ATAN2(Y,X) + 6.2831853,6.2831853) DEC = ASIN(Z/R) DLE = AMOD(ATAN2(Y*CE + Z*SE,X) + 6.2831853,6.2831853) SZE = ASIN((Z*CE - Y*SE)/R) RETURN END |
| W sierpniu 2002 r. opisany program i wszystkie podprogramy ponownie skompilowalem starym Laheyowym F77L (v. 2.2) oraz nowszym LF90 (v. 4.5) nie napotykając żadnych problemów i uzyskując te same wyniki, które widnieją w tabelce oryginalnego artykułu. Poniżej dołączam dla kompletu podprogramy (w takiej postaci, jaką miały wówczas) na obliczanie pozycji pozostałych planet, które nie były zamieszczone w przedstawionym artykule, ze względu na objętość. |
SUBROUTINE MERKUR(RA,DEC,DLE,SZE,R)
COMMON /PLAN/ T,A1,A7,A8,A13,A16,A19
DOUBLE PRECISION T,D1,D2
F(D1,D2) = 6.28318530718D0*DMOD(D1 + T*D2,1D0)
A9 = F(.700695D0,.011367714D0)
A10 = F(.485541D0,.01136759566D0)
A11 = F(.566441D0,.01136762384D0)
DLH = A9+((84386.+2.2E-4*SNGL(T))*SIN(A10)+10733.*SIN(A10+
* A10)+1892.*SIN(3.*A10)-646.*SIN(A11+A11)+381.*SIN(4.*A10)-306.*
* SIN(A10-A11-A11)-274.*SIN(A10+A11+A11)-92.*SIN(2.*(A10+A11))+83.*
* SIN(5.*A10)-28.*SIN(3.*A10+2.*A11)+25.*SIN(2.*(A10-A11))+
* 19.*SIN(6.*A10)-9.*SIN(2.*(A10+A10+A11))+7.*COS(2*A10-5.*A13))/
* 206264.8
SZH = (24134.*SIN(A11)+5180.*SIN(A10-A11)+4910.*SIN(A10+A11)+
* 1124.*SIN(A10+A10+A11)+271.*SIN(3.*A10+A11)
* +132.*SIN(A10+A10-A11)+67.*SIN(4.*A10+A11)+18.*SIN(3.*A10-A11)+
* 17.*SIN(5.*A10+A11)-10.*SIN(3.*A11)-9.*SIN(A10-3.*A11))/
* 206264.8
RH = .39528-.07834*COS(A10)-.00795*COS(2.*A10)-
* .00121*COS(3.*A10)-22.E-5*COS(4.*A10)
CALL HELGEO(DLH,SZH,RH,RA,DEC,DLE,SZE,R)
RETURN
END
SUBROUTINE MARS(RA,DEC,DLE,SZE,R)
COMMON /PLAN/ T,A1,A7,A8,A13,A16,A19
DOUBLE PRECISION T,D1,D2
F(D1,D2) = 6.28318530718D0*DMOD(D1 + T*D2,1D0)
A15 = F(.987353D0,.145575328D-2)
A17 = F(.849694D0,.145569465D-2)
DLH = A15+((38451.+(37.+8.*COS(A16))*(1.+SNGL(T)/36525.))*
* SIN(A16)+2238.*SIN(A16+A16)+181.*SIN(3.*A16)-52.*SIN(2.*A17)-
* 22.*COS(A16-2.*A19)-19.*SIN(A16-A19)+17.*(COS(A16-A19)+SIN(4.*
* A16))-16.*COS(2.*(A16-A19))+13.*COS(A8-A16-A16)-10.*(SIN(A16-
* 2.*A17)+SIN(A16+2.*A17))+7.*(COS(A8-A16)-COS(2.*A8-4.*.75*A16))-
* 5.*(SIN(A13-3.*A16)+SIN(A8-A16)+SIN(A8-2.*A16))+4.*(-COS((2.*A8-
* 4.*A16))+COS(A19))+3.*(COS(A13-3.*A16)+SIN(2.*(A16-A19))))
* /206264.8
SZH = (6603.*SIN(A17)+622.*SIN(A16-A17)+615.*SIN(A16+A17)+
* 64.*SIN(2.*A16+A17))/206264.8
RH = 1.53031-.1417*COS(A16)-.0066*COS(2.*A16)-.00047*COS(3.*A16)
CALL HELGEO(DLH,SZH,RH,RA,DEC,DLE,SZE,R)
RETURN
END
SUBROUTINE JOWISZ(RA,DEC,DLE,SZE,R)
COMMON /PLAN/ T,A1,A7,A8,A13,A16,A19
DOUBLE PRECISION T,D1,D2
F(D1,D2) = 6.28318530718D0*DMOD(D1 + T*D2,1D0)
W = T/36525.+1.
A25 = F(.400589D0,.3269438D-4)
A22 = F(.882987D0,.9294371D-4)
A18 = A19+.2078289
A0 = 2.*A19-5.*A22
S = SIN(A19)
C = COS(A19)
C0 = COS(A0)
S0 = SIN(A0)
S1 = SIN(A19-A22)
C1 = COS(A19-A22)
S2 = SIN(A22)
C2 = COS(A22)
DLH = A18+(W*(5023.-74*C+68.*S-43.*S0-19.*C0-4.*(C*C-S*S)+6.*S*C
* -2.*S0*C)+19934.*S+2511.+1093.*C0+1202.*S*C-479*S0-370.*S1*C1+
* C*(137.*S0-262.*S1*C1)+S*(137.*C0+131.*(C1*C1-S1*S1))+79.*C1-
* 76.*(C1*C1-S1*S1)+66.*COS(2.*A19-3.*A22)+116.*C0*C-10.*S0*S-37.*C
* +49.*SIN(2.*A19-3.*A22)+25.*(SIN(A18+A18)+SIN(3.*A19))-
* 23.*SIN(A0-A19)+17.*(COS(A0+A22)+COS(3.*(A19-A22)))-14.*S1-
* 13.*SIN(A0+A19+A22)+9.*(C2-COS(A18+A18)-S2-SIN(3.*A19-A22-A22)+
* SIN(A0+A19+A19)+SIN(2.*(A19-3.*A22)+3.*A25))-8.*(C0*C0-S0*S0)+
* 7.*(COS(A0+A19+A22)-COS(A19-3.*A22)-2.*S0*C0-SIN(A19-3.*A22))+
* 6.*(COS(A0+A19+A19)-SIN(3.*(A19-A22)))+5.*(C2*C2-S2*S2)-
* 4.*(4.*S1*C1*(C1*C1-S1*S1)-C*SIN(A18+A18)+COS(3.*A22)-COS(2.*A19-
* A22)+COS(3.*A19-2.*A22))+3.*(COS(5.*A22)+COS(5.*(A19-2.*A22))+
* S2))/206264.8
SZH = (-4692.*C+259.*S+454.*S*S+W*(30.*S+21.*C)+29.*S*C0+
* C*(3.*(S0-4.+8.*S*(S+S+1.))+C0+C0)-12.*S*S0)/206264.8
RH = 5.20883-.25122*C-.00604*(C*C-S*S)+.0026*COS(2.*(A19-A22))+
* 1.E-5*(-170.*COS(A0+A19)-106.*SIN(2.*(A19-A22))-W*(91.*S+84.*C)+
* 69.*SIN(A0+A22+A22)-67.*SIN(A0-A19)+63.*S1-51.*COS(A0+A22+A22)-
* 46.*S+S0*(66.*C-29.*S)+C0*(66.*S-29.*C)+27.*COS(A19-2.*A22)-
* 22.*C*(1.-4.*S*S)-21.*S0)
CALL HELGEO(DLH,SZH,RH,RA,DEC,DLE,SZE,R)
RETURN
END
SUBROUTINE SATURN(RA,DEC,DLE,SZE,R)
COMMON /PLAN/ T,A1,A7,A8,A13,A16,A19
DOUBLE PRECISION T,D1,D2
F(D1,D2) = 6.28318530718D0*DMOD(D1 + T*D2,1D0)
W = T/36525.+1.
A25 = F(.400589D0,.3269438D-4)
A22 = F(.882987D0,.9294371D-4)
A21 = A22+1.5727316
A0 = 2.*A19-5.*A22
C0 = COS(A0)
S0 = SIN(A0)
S = SIN(A22)
C = COS(A22)
DLH = A21+(W*(5014.-229.*C-142.*S+101.*S0+60.*C0+22.*SIN(A0+A22)
* -22.*S*C-17.*(C*C-S*S)-6.*(2.*S0*(S0+C0)-1.-SIN(A0-A22))+
* 4.*COS(A19-A22-A22)+3.*(COS(2.*A19-4.*A22)-SIN(A21+A21)-COS(A21+
* A21)+COS(A19+A19-6.*A22))+2.*SIN(A19-A22-A22))+
* 23045.*S-2689.*C0+2507.-826.*COS(A0+A22)+1604.*S*C+1177.*S0+
* 425.*SIN(A19-A22-A22)-153.*COS(A0-A22)-114.*C-70.*COS(A21+A21)+
* 67.*SIN(A21+A21)+66.*SIN(A0-A22)+41.*SIN(A19-3.*A22)+39.*SIN(3.*
* A22)+31.*(SIN(A19-A22)+ SIN(2.*(A19-A22)))-29.*COS(A19+A19-3.*
* A22)-28.*SIN(A0-A22+3.*A25)+28.*COS(A19-3.*A22)-
* 22.*SIN(A22-3.*A25)+20.*(SIN(2.*A19-3.*A22)+COS(2.*A0))
* +19.*(COS(2.*A22-3.*A25)+2.*S0*C0)-16.*COS(A22-3.*A25)-
* 12.*(SIN(A0+A22)-COS(A19)+SIN(2.*(A22-A25)))-11.*COS(A0-A22-A22)+
* 10.*(SIN(A22+A22-3.*A25)+COS(2.*(A19-A22)))+9.*SIN(4.*A19-9.*A22)
* +8.*(COS(A21*2.-A22)-COS(A21*2.+A22)-SIN(A22-2.*A25)-SIN(2.*A21-
* A22)+COS(A22+A25))+7.*(SIN(2.*A21+A22)-COS(A19-2.*A22)-COS(A22*2.
* ))+5.*(SIN(A0-A22-A22)+SIN(3.*A19-4.*A22)-SIN(3.*A19-7.*A22)-
*COS(3.*(A19-A22))-COS(2.*(A22-A25)))+4.*(SIN(3.*(A19-A22))+SIN(A0+
* A19))+3.*(COS(2.*A19-6.*A22+3.*A25)+COS(3.*A19-7.*A22)+COS(4.*A19
* -9.*A22)+SIN(3.*(A19-A22-A22))+SIN(2.*A19-A22)+SIN(A19-4.*A22))
* +2.*(COS(3.*(A22-A25))+SIN(4.*A22)-COS(3.*A19-4.*A22)-COS(2.*
* A19-A22)-SIN(2.*A19-7.*A22+3.*A25)+COS(A19-4.*A22)+COS(2.*A0-A22)
* -SIN(A22-A25)))/206264.8
SZH = (8297.*S-3346.*C+924.*S*C-189.*(C*C-S*S)+185.+
* W*(79.*C+18.*S-10.-8.*S*C+3.*SIN(A0+A22))-71.*COS(A0+A22)+
* 46.*SIN(A0-A22)-45.*COS(A0-A22)+29.*SIN(3.*A22)-
* 20.*COS(A0+A22+A22)-14.*COS(A0)-11.*COS(3.*A22)+9.*SIN(A19-3.*
* A22)+8.*SIN(A19-A22)-6.*SIN(A0+A22+A22)+5.*(SIN(A0-A22-A22)-
* COS(A0-A22-A22))+4.*S0+3.*SIN(A19-A22-A22)*(1.+S+S)+
* 2.*(SIN(4.*A22)-COS(2.*(A19-A22))))/206264.8
RH = 9.55774+W*(-.00524*S+.00328*C-.00028)+
* 1.E-5*(-53252.*C-1878.*SIN(A0+A22)-1482.*(C*C-S*S)+817.*
* SIN(A19-A22)-539.*COS(A19-A22-A22)+349.*S0+347.*SIN(A0-A22)-
* 225.*S+149.*COS(A0-A22)-126.*COS(2.*(A19-A22))+104.*COS(A19-A22)+
* 101.*C0+98.*COS(A19-3.*A22)-73.*COS(A0+A22+A22)-62.*COS(3.*A22)+
* 42.*SIN(2.*A22-3.*A25)+41.*SIN(2.*(A19-A22))-40.*(SIN(A19-3.*
* A22)-COS(A0+A22))-23.*SIN(A19)+20.*SIN(A0-A22-A22))
CALL HELGEO(DLH,SZH,RH,RA,DEC,DLE,SZE,R)
RETURN
END
SUBROUTINE URAN(RA,DEC,DLE,SZE,R)
COMMON /PLAN/ T,A1,A7,A8,A13,A16,A19
DOUBLE PRECISION T,D1,D2
F(D1,D2) = 6.28318530718D0*DMOD(D1 + T*D2,1D0)
A22 = F(.882987D0,.9294371D-4)
A25 = F(.400589D0,.3269438D-4)
A26 = F(.664614D0,.3265562D-4)
A28 = F(.725368D0,.1672092D-4)
W = T/36525.+1.
S = SIN(A25)
C = COS(A25)
C2 = COS(A22-A25-A25)
S2 = SIN(A22-A25-A25)
S3 = SIN(A22-3.*A25)
C3 = COS(A22-3.*A25)
S1 = SIN(A22-A25)
DLH = A25+2.950458+(19397.*S+570.*2.*S*C+143.*S2+102.*S3+76.*C3-
* 49.*SIN(A19-A25)+29.*(SIN(2.*A19-6.*A22+3.*A25)+COS(2.*(A25-A28))
* )-28.*COS(A25-A28)+23.*SIN(3.*A25)-21.*COS(A19-A25)+
* 20.*(SIN(A25-A28)+C2)-19.*COS(A22-A25)+17.*SIN(2.*A25-3.*A28)+
* 14.*SIN(3.*(A25-A28))+13.*S1-12.*C+10.*SIN(2.*(A25-A28))-
* 9.*(SIN(A26+A26)-COS(A25+A25-3.*A28))+8.4853*SIN(2.*(A19-3.*A22+
* A25+.3926991))+5.*SIN(A22-4.*A25)-4.*(SIN(3.*A25-4.*A28)-
* COS(3.*(A25-A28)))-3.*COS(A28)-2.*SIN(A28)+
* W*(110.*S-536*C-30.*(C*C-S*S)+8.*C2+7.*(C3-S3+2.*S*C)+
* W*(32.-12.*C-9.*S)))/206264.8
SZH = ((2775.-C)*SIN(A26)+261.*S*COS(A26))/206264.8
RH = 19.21216-.90154*C-.02121*(C*C-S*S)-.00585*C2 -
* .00451*COS(A19-A25)+.00336*S1+.00198*SIN(A19-A25)+
* .00118*C3+.00107*S2-.00081*COS(3.*(A25-A28))-
* W*(S*.02488+C*.00508+.00206*S*C)
CALL HELGEO(DLH,SZH,RH,RA,DEC,DLE,SZE,R)
RETURN
END
SUBROUTINE NEPTUN(RA,DEC,DLE,SZE,R)
COMMON /PLAN/ T,A1,A7,A8,A13,A16,A19
DOUBLE PRECISION T,D1,D2
F(D1,D2) = 6.28318530718D0*DMOD(D1 + T*D2,1D0)
A22 = F(.882987D0,.9294371D-4)
A25 = F(.400589D0,.3269438D-4)
A28 = F(.725368D0,.1672092D-4)
A29 = F(.480856D0,.1663715D-4)
W = T/36525.+1.
S = SIN(A28)
C = COS(A28)
S9 = SIN(A19-A28)
C9 = COS(A19-A28)
DLH = (3523.*S-W*(43.*C+4.*S*(1.-W))-50.*SIN(A29+A29)+
* 29.*S9+38.*S*C-18.*C9+18.38478*SIN(A22-A28+.7853982)-
* 9.*(SIN(A25+A25-3.*A28)-COS(2.*(A25-A28)))-5.*COS(A25+A25-3.*
* A28)+4.*COS(A25-A28-A28))/206264.8+A28+.7636834
SZH = ((6404.-33.*W)*SIN(A29)+110.*S*COS(A29))/206264.8
RH = 30.07175-W*.00314*S-.25701*C-.00787*COS(A25-A28-A28+4.3735793
* )+.00409*C9+.0025*S9-.00194*SIN(A22-A28)+.00185*COS(A22-A28)
CALL HELGEO(DLH,SZH,RH,RA,DEC,DLE,SZE,R)
RETURN
END
SUBROUTINE PLUTON(RA,DEC,DLE,SZE,R)
COMMON /PLAN/ T,A1,A7,A8,A13,A16,A19
DOUBLE PRECISION T,D1,D2
F(D1,D2) = 6.28318530718D0*DMOD(D1 + T*D2,1D0)
W = T/36525.+1.
A31 = F(.663854D0,1.115482D-5)
A32 = F(.04102D0,1.104864D-5)
A33 = A32+1.987591
S = SIN(A32)
C = COS(A32)
S1 = SIN(A33)
C1 = COS(A33)
DLH = A31+(S*(106892.+S*(S1*C1*(3900.+2280.*C)-S*(C*6712+12516.))
* +C*33312.+W*(W*227.+200.)+S1*S1*(7574.+C*(2156.+2280.*C)))-
* S1*C1*(9136.- 90.*C))*4.8481368E-6
SZH = (57726.*S1+29359.*S*C1-1155.*C*S1+3870.*SIN(A32+A32+A33)+
* 1138.*SIN(3.*A32+A33)+472.*SIN(A32+A32-A33)+353.*SIN(4.*A32+A33)-
* 255.*S*C1*(1.-4.*S1*S1)+(33.*C-119.)*S1*(4.*C1*C1-1.))/206264.8
RH = 41.86345-C*(8.90288+C*(1.93438+C*(.90596+.39968*C)))
CALL HELGEO(DLH,SZH,RH,RA,DEC,DLE,SZE,R)
RETURN
END |