0
KEPLERSE BANEN ONCONVENTIONEEL Mathematische verhandeling voor een alternatieve afleiding van de perkenwet door: ing. R.C. Ott
21 december 1998
1 INHOUDSOPGAVE Blz. INLEIDING EN BESCHOUWING
2
1. DE MATHEMATISCHE ELLIPSVORM
3
2. PERKENWET IN CONVENTIONELE NOTATIE ONCONVENTIONEEL BEREKEND
4
2.1 METRISCHE AFLEIDING ALS FUNCTIE VAN DE ANOMALIE
4
2.2 ALFLEIDING NAAR HET TIJDSDOMEIN:
7
2.3 DE KWALITEIT EN NAUWKEURIGHEID VAN DE BESCHREVEN BAANAFLEIDING
9
2.3.1 COMPUTER NAUWKEURIGHEID.
11
2.3.2 PRECISIE VAN DE WARE ANOMALIE
12
2.4 DE BRUIKBAARHEID VAN DE METHODE:
12
2.5 VOORBEELD BEREKENING
13
3. ELLIPSBAAN CURIOSA
15
3.1 DE GEMIDDELDE VOERSTRAAL
15
3.2 DE TIJDSGEMIDDELDE VOERSTRAAL
15
3.3 DE HOEKGEMIDDELDE VOERSTRAAL
17
BIJLAGEN BIJLAGE I PROGRAMMA ROUTINES IN PASCAL
18
BIJLAGE II DE TIJDSGEMIDDELDE VOERSTRAAL
21
REFERENTIES
23
2 INLEIDING EN BESCHOUWING In dit rapport geef ik een verhandeling van een alternatieve interpretatie van de 2e wet van Kepler. Deze mathematische afleidingen van de Keplerse banen zijn niet conventioneel, maar geven naar mijn inzicht een bruikbare methode. Hoe dan ook de methode is in ieder geval illustratief. Ik heb voor de juistheid van de diverse (tussen) afleidingen gebruik gemaakt van een groot aantal zelf geschreven test routines in PASCAL (Delphi voor Windows) en het tekenprogramma AUTOCAD r13 is hierbij een belangrijk hulpmiddel geweest om inzichten te verkrijgen. Alsmede hebben almanakgegevens, gegevens op internet en een aantal conventionele afleidingen gediend ter vergelijk. Ik heb getracht een zo beknopt mogelijk duidelijk beeld te geven van de eigenschappen en de kwaliteiten van de afgeleide methodiek alsook de wijze waarop de afleidingen tot stand zijn gekomen. Ook kom ik tot een aantal opmerkelijke mathematische verschijnselen waarvan ik op dit moment het nut ervan nog niet precies weet. Ik heb enige van die opmerkelijkheden als (voor mij) curiosa in hoofdstuk 3 beschreven. Ook is in dit verslag de broncode opgenomen van een aantal PASCAL routines voor het beschrijven van de Keplerse banen volgens de afgeleide methode. Dit verslag heeft voor mij gedeeltelijke een open eind; nieuwe bevindingen en informatie kunnen nieuwe inzichten geven en een aanvullingen op dit verslag. Reinier C. Ott December 1998
Ing. R.C. Ott Uiverweide 33 6708 LA Wageningen 0317-410492 e-mail
[email protected] URL: http://www.dutch.nl/rcott/astronom.htm
3 1. DE MATHEMATISCHE ELLIPSVORM Een ellips is gedefinieerd als een vlakke kromme, bestaande uit een verzameling punten waarvoor de som der afstanden tot 2 vaste punten, de zo geheten brandpunten van de ellips, constant is. (en wel gelijk is aan de lange as van de ellips) (ref.: Vademecum van de Wiskunde , Prisma) Voor het beschrijven van elliptische (Keplerse) banen is het nodig om de voerstraal te berekenen als functie van de doorlopen hoek vanuit een der brandpunten. Uit de definitie van de ellips kan een ellips worden geconstrueerd met behulp van een niet rekbaar touwtje waarvan het begin en einde aan elkaar zijn geknoopt en waarvan de beide brandpunten en een te tekenen punt van de ellips met dit touwtje steeds een strak gespannen driehoek vormen. De onderstaande figuur toont hoe de voerstraal R vanuit een der brandpunten kan worden geconstrueerd als functie van de hoek ϕ.
Figuur 1.1 Voor het doorlopen van de hoek ϕ om centrum 1 geldt: • De lengte van het denkbeeldige touwtje ( L+R+S ) blijft steeds constant. En daar ook de brandpuntsafstanden ongewijzigd blijven zal dus tevens R+S constant blijven. C = R+S ( C is de constante) • L = Rmax - Rmin R+S = L + 2.Rmin =>
C = R+S = Rmax + Rmin
• S2 = R2 + L2 - 2 R L COSϕ met S = C - R is af te leiden:
(Cosinusregel)
[1.1]
4 2. PERKENWET IN CONVENTIONELE NOTATIE ONCONVENTIONEEL BEREKEND. 2.1: METRISCHE AFLEIDING ALS FUNCTIE VAN DE ANOMALIE In de klassieke astronomie wordt de radiusvector in de elliptische vergelijking op een andere wijze gerepresenteerd dan het voorgaande hoofdstuk en wel zodanig dat de ware anomalie te noemen “ν” direct wordt beschreven vanuit het perifocus van de ellips. Zonder deze vergelijking voor de voerstraal af te leiden luidt deze: P R(ν) =
[2.1.1] 1 + e⋅COS ν
De parameter P wordt de semi-parameter genoemd, dit is de elliptische vormfactor en wordt berekend uit de excentriciteit (e) en de perifocus (of perihelium) van de ellips (q) volgens: P = q⋅(1 + e)
[2.1.2]
Of desgewenst uit het peri- en apfocus (Q) : P = 2⋅q⋅Q / (q+Q)
[2.1.3]
Voor het beschrijven van de ware anomalie geldt volgens Kepler steeds dat in gelijke tijden gelijke perken (= oppervlakken vanuit het te beschouwen brandpunt) worden beschreven. Er zal moeten worden getracht om een algemene vergelijking voor deze perken op te stellen alsook voor functie van ν. In de onderstaande figuur bevindt de Zon zich in het bewuste brandpunt van de elliptische baan die door de planeet wordt beschreven :
Figuur 2.1.1 Voor de gearceerde oppervlakken in de bovenstaande figuur geldt nu: A1 / t1 = A2 / t2
; dAν /dt = k = Atot/P
Let op: verwar P niet met P !
Atot P
= Totale oppervlak van de ellips = Omloopperiode van de planeet
5 Wanneer een (oneindig) klein segmentje op een willekeurige plaats in de ellips wordt beschouwd mag worden gesteld dat het oppervlak ervan ongeveer voldoet aan de halve basis maal hoogte. In het geval dat de hoogte van het driehoekje voldoende klein is mag ook de hoogte gelijk worden gesteld aan het stukje van de doorlopen boog R⋅dν , ν moet natuurlijk wel in radialen worden uitgedrukt.
Wanneer de algemene vergelijking [2.1.1] van Rν met het bovenstaande wordt betrokken geldt:
Voor het gehele gearceerde oppervlak tussen de grenzen a en b geldt na integratie de volgende vergelijking voor het (perk)oppervlak als functie van de ware anomalie ν.
[2.1.4] Een willekeurig perkoppervlak Aa-b waarvan de gens a niet samenvalt met het perifocus wordt berekend uit het verschil oppervlak: A a-b = A q-b - A q-a Hierin is A q-b , A q-a het perkoppervlak tussen de grens b respectievelijk grens a en het perifocus van de ellips.
6 Wanneer de grenzen va n de anomalie niet specifiek zijn aangegeven geldt vervolgens:
[2.1.5]
Na partiële integratie wordt de kwadraatterm gereduceerd.
Ook voor de gereduceerde integraal blijkt (gelukkig) een oplossing:
Hiermee wordt de totale vergelijking Aν door substitutie:
[2.1.6] Opmerking: De functie Aν is slechts geldig voor 0 < ν < π : Uit de symmetrie van de ellips over de lange as kan de anomalie ν over 2π worden berekend. De constante factor aan de voorzijde van de accolades is vereenvoudigd tot het halve product van het peri- en aphelium. Hierbij deze vereenvoudigingswijze: P2 e2 - 1
P2 =
2qQ/(q+Q)⋅ 2qQ/(q+Q) =
(e-1)⋅(e+1)
= q⋅Q
[2.1.7]
2q/(q+Q) ⋅ 2Q/(q+Q)
Opmerkelijk is hierbij dat het product q⋅Q als een oppervlakte schaalfactor voor de ellips mag worden beschouwd, de parameter e zegt uitsluitend iets over de mate van de ellipticiteit. Dit betekend, zonder echter een sluitend bewijs, dat ook het totale oppervlak van de ellips moet voldoen aan : Atot = π⋅a⋅b [2.1.8] ( Dit is de gebruikelijke manier om het oppervlak van een totale ellips te berekenen) Hierin is de betekenis van: a is de halve lange as met als mathematisch verband: a = q /(1-e) of a = Q/(1+e) b is de halve korte as en is te beschrijven als: b = a⋅
1 - e2
[2.1.9] [2.1.10]
Ter illustratie volgt hieronder een bewijs dat de totale ellips welke berekend wordt door vergelijking [2.1.6] voldoet aan de algemene bewering van verg.[2.1.8] Met het gegeven voor a mag b ook worden geschreven als :
7
Q
1 - e2
Q (1 - e)
b=
of desgewenst als: b = 1+e
[2.1.11] 1 - e2
Wanneer nu het oppervlak van de ellips wordt berekend voor ν → π volgens vergelijking [2.1.6] geldt voor het totale oppervlak Atot = 2 ⋅ | A(ν→ π) | Let op: Het domein van ν is voor de vergelijking gedefinieerd tussen 0 - π, π zelf niet meegerekend! Hierna wordt vergelijking [2.1.6] teruggebracht tot:
π⋅q⋅Q Atot =
[2.1.12]
1 - e2
Wanneer q en Q worden gesubstitueerd door respectievelijk vergelijking [2.1.9] en [2.1.11] b⋅ dus met q = a⋅( 1 - e)
en
1 - e2
Q= 1- e
blijkt dus direct te gelden dat A(ν=2⋅π) = π ⋅a ⋅b
2.2 ALFLEIDING NAAR HET TIJDSDOMEIN: Tot nu toe is steeds gesproken over een metrische betrekking; de plaats van een hemellichaam is beschreven in relatie met het doorlopen oppervlak (perk) van haar (elliptische) baan. Voor verdere berekeningen van de 2e wet van Kepler wil men graag weten wat de ware anomalie ν van een hemellichaam is als functie van de tijd. Deze tijd wordt gerekend vanaf het moment dat de perihelium passage plaats vindt. Het probleem zal dus zijn om ν uit vergelijking [2.1.6] terug te rekenen en daarvoor zal het nodig zijn om middels een geschikte iteratie methode de juiste waarde voor ν te verkrijgen. Zoals reeds vermeld luidt de 2e wet van Kepler : dAν /dt = k = Atot/P
Atot P
= Totale oppervlak van de ellips = Omloopperiode van de planeet
[2.2.1]
De aldus berekende ellips-sector (perk) is een functie van ν maar eveneens van de tijd. Alleen bij een cirkelvormige baan is de ware anomalie evenredig in tijd en wordt in dat geval uitgedrukt als middelbare anomalie M . Dus slechts bij een cirkelvormige baan geldt : ν = M. Echter M is wel een goed uitgangspunt voor een goede begin schatting voor de iteratie methode, teneinde de gewenste waarde met voldoende nauwkeurigheid voor ν te verkrijgen. Uitgaande voor de cirkelvormige baan geldt dan een evenredig verband tussen het doorlopen oppervlak A(t) van de baan in de tijd. A(t) = ½ R2 . M(t) M(t)
= 2π . t / P
R
= (constante) voerstraal hemellichaam - Zon
( M in Rad )
[2.2.2] [2.2.3]
In het geval van de werkelijke elliptische baan is ν(t) niet evenredig meer met de tijd , echter het enige dat volgens Kepler ongewijzigd blijft voor de vertaling van een cirkelvormige naar een elliptische baan is het oppervlak van de doorlopen perk. A(t) = k ⋅ t
8 Voor het verder doorrekenen, uitgaande van de cirkelvormige baan wordt voor de eerste iteratiestap gesteld : ν(t) = M(t) ! Nu zal moeten worden getracht opnieuw ν(t) te berekenen , geldend voor de perk van de elliptische baan. Voor dit terugrekenen naar het juiste argument wordt de iteratie methode van Newton-Raphson toegepast omdat deze snel convergeert en de ingrediënten bezit die reeds bekend zijn.
νi+1
=
νi
f (νi) -
[2.2.4]
f ’(νi)
Normaliter wordt deze methode toegepast om te berekenen waar de functie f(ν) de nuldoorgang(en) heeft. Echter nu moet worden bepaald wat het argument is voor de waarde dat de functie voldoet aan het gestelde oppervlak. Hiervoor wordt f (νi) geschreven als een variabele E die de foutterm beschrijft E i = A(νi) - A(t)
A(νi) wordt berekend uit de afgeleide vergelijking [2.1.6] A(t) wordt eenvoudig uit het evenredige verband met de tijd bepaald en behoud dus tijdens het iteratieproces dezelfde waarde.
Er dient dus net zo lang geïtereerd te worden dat E i een voldoende kleine waarde heeft, ofwel dat het gezochte oppervlak met voldoende nauwkeurigheid gelijk is aan A(t) Verder is nodig de afgeleide functie van A(ν). En omdat A(ν) in beginsel als primitieve is geschreven (zie algemene bewering [2.1.4] van de ellips ) is dA/dν eenvoudig te berekenen. Dit laatste is het grote voordeel om de Newton-Raphson methode toe te passen. De iteratieformule wordt hiermee:
ν i+1
=
νi
A(νi) - A(t)
-
[2.2.5] dAi / dνi
Hierin de verandering van Aν als functie van ν:
[2.2.6]
Na enkele iteratieve berekeningen wordt al een zeer hoge nauwkeurigheid verkregen. De inverse functie van ν : Omgekeerd kan de perk doorlooptijd t exact worden berekend wanneer het argument van ν hiervoor gegeven wordt. Dit is dan de tijd die sinds de perihelium passage is verstreken tot het moment waar zich de planeet (of een ander hemellichaam) zich in haar baan bevindt.
9 2.3 DE KWALITEIT EN NAUWKEURIGHEID VAN DE BESCHREVEN BAANAFLEIDING Voor de bepaling van de ware anomalie ν is dus een iteratie proces nodig. Belangrijk hierbij is de mate waarin dit proces convergeert, met andere woorden het aantal herhalingsberekeningen dat noodzakelijk is om de gewenste nauwkeurigheid van ν te bereiken. Veelal het grootste probleem bij een numerieke benadering is de keuze van de beginwaarde om het iteratieve proces te starten. Zoals vermeld wordt in het voorgaande de 1e stap van ν berekend uit de lineaire verhouding in de tijd gerelateerd aan de totale omlooptijd. Deze lineaire verhouding is slechts van toepassing indien de baan van het hemellichaam zuiver cirkelvormig is en dus wordt uitgedrukt als middelbare anomalie M. Daar een hemellichaam zich op ieder willekeurige plaats in haar baan kan bevinden zal dit gehele bereik getoetst moeten worden. Echter omdat de werkelijke ellipsbaan symmetrisch is om de absiden-lijn kan worden volstaan voor een toetsing van M=0° tot 180°. Evenzo zal moeten worden onderzocht in welke mate van ellipsvormigheid nog sprake kan zijn om ν met voldoende snelheid en nauwkeurigheid te berekenen. De resultaten van het onderzoek zijn in de onderstaande figuur weergegeven. De figuur laat een 2 dimensionale fractal zien waarbij de kleur en de intensiteit het aantal iteraties weergeeft om een nauwkeurigheid van ∆ν = 10-9 Rad te bereiken. Voor de berekening is gebruik gemaakt van de PASCAL routines Bijlage I.
Figuur 2.3.1 De toenemende ellipticiteit wordt van beneden naar boven weergegeven. De onderzijde van de fractal stelt het aantal iteraties voor dat nodig is voor het berekenen van een cirkelvormige baan, terwijl de bovenzijde ervan wordt gekenmerkt door extreme ellipsvormen.
10 In figuur 2.3.1 naar rechts is de toenemende middelbare anomalie voor over het bereik van 0 tot 180° weergegeven. De betekenis van de kleuren is als volgt, waarbij deze kleur een weergave is van het aantal iteraties dat nodig is om ν te berekenen: Blauwe kleur Groen Geel Rood
: : : :
0 tot en met 10 iteraties 10 t/m 25 iteraties 25 t/m 100 iteraties Meer dan 100 iteraties
Bovendien geldt dat naarmate de kleur lichter van tint is het aantal iteraties toeneemt. Het aardige van een fractal is dat er kan worden ‘ingezoomd’ , de onderstaande figuur toont een uitvergroting van het witte omkaderde gebied voor de grenzen ( e = 0.58 tot 0.68 en M = 50° tot 90°)
Figuur 2.3.2 Opvallend is dat bij uitvergroten, dus met een toenemend aantal (tussen) waarden, opnieuw een aantal ‘rode’ gebiedjes ontstaan. Er ontstaat dus kennelijk een duidelijke scheiding tussen de ‘stabiele’ blauwe gebieden en de meer chaotische gedetailleerde plaatsen. Uit de figuren is af te lezen dat de excentriciteit van de ellipsbaan, tot een waarde van e =0.6, voor elke middelbare anomalie minder dan 10 iteraties nodig zijn om de gewenste eindnauwkeurigheid voor ν te bereiken. Dit alles kan nog worden bevestigd in de onderstaande tabel. In deze tabel wordt de excentriciteit van de ellips uitgezet tegen het maximaal aantal iteraties dat optreedt tijdens de gehele omloop van de het hemellichaam. Excentriciteit e 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45
max. aantal Iteraties 1 4 4 4 5 5 5 5 6 6
Excentriciteit e 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95
max. aantal Iteraties 7 7 8 46 64 75 105 >1000 242 352 Tabel 2.3.1
Voor het doorlopen van M is een stapgrootte ∆M van 0.025° gehanteerd.
11 Uit tabel 2.3.1 is direct de discrepantie af te lezen in het aantal iteraties dat optreedt wanneer een excentriciteit van e=0.6 wordt overschreden. Direct hieruit kan de vraag worden gesteld of hierin een verschuiving optreedt wanneer de nauwkeurigheidseis van ∆ν = 10-9 wordt gewijzigd. Tevens zal moeten worden getoetst of de computer zelf in staat is ν met de gewenste nauwkeurigheid te berekenen. 2.3.1 COMPUTER NAUWKEURIGHEID. De berekeningen zijn uitgevoerd in eerste instantie op basis van de drijvende punt variabelen “REAL” en in tweede instantie vervangen door het PASCAL type “DOUBLE” en “EXTENDED”. In PASCAL zijn deze gegevenstypen als volgt gedeclareerd: Type
Waardebereik (absoluut)
REAL
2.9x10
-39
- 1.7x10
+38
5.0x10-324 - 1.7x10+308
DOUBLE EXTENDED
3.4x10
-4932
- 1.1x10
+4932
Aantal bytes 6 8 10 Tabel 2.3.1.1
Wanneer alle drijvende punt variabelen in de PASCAL routines (bijlage 1.1) worden gewijzigd in achtereenvolgens de bovenstaande typen dan ontstaat het onderstaande resultaat.
Figuur 2.3.1.1 REAL
DOUBLE
EXTENDED
Op een of andere wijze levert het DOUBLE en EXTENDED type een marginale afronding ten ongunste van het iteratieproces. De verschillen blijven klein en zijn alleen in het kleurrijke chaotische gebied waar te nemen. Anderzijds moet ook worden afgevraagd hoe de goniometrische functies worden berekend, of beter gezegd, door de computer worden benaderd. Het verschil tussen de drijvende punt typen uitten zich slechts op enkele plaatsen waar e de waarde van 1 nadert en dus waar een hoge mate van chaos heerst, voor de rest zijn de fractals identiek. Voor vergelijking [2.1.6] treedt instabiliteit op naarmate e toeneemt naar 1 en ν >90° In het limiet geval dat e=1 geldt: p ⋅Q
SIN ν ⋅{
Aν = 2
-
∞ ⋅ ARCTAN 0 }
[2.3.1.1]
1 + COS ν
Mathematisch geldt natuurlijk dat de ARCTAN 0 in het interval -π/2 tot +π/2 een argument van 0 oplevert. Dus hiermee verkrijgt het tweede lid de waarde van ∞ ⋅ 0 wat elke waarde kan opleveren. Wanneer e de waarde van 1 nadert is maar de vraag in hoeverre de ARCTAN functie in staat is de zie verg. [2.1.6] ) te kunnen compenseren. Een steeds groter wordende wortelterm ( 2/ 1 - e2 voordeel is dat de term niet snel divergeert naar oneindig (pas bij zeer sterke elliptische banen (e>0.9999…) is dit het geval. Anderzijds levert de SIN /( 1 + COS ) term een probleem voor ondefinieerbaarheid op wanneer e->1 tezamen met respectievelijk ν -> 180° Tevens is getracht de goniometrische functies SIN, COS en ARCTAN te benaderen met een polynoom-reeks. Echter het probleem is de slechte mate van convergentie van de inverse boog functie arctangens. Opmerkelijk is echter steeds dat de gebieden zelf onveranderd blijven waar stabiliteit en chaos heerst. En omwille van de marginale afwijkingen die kunnen optreden kan in de berekeningen worden volstaan met het ‘normale’ REAL type.
12 2.3.2 PRECISIE VAN DE WARE ANOMALIE In het voorgaande is een eindnauwkeurigheid in ν gesteld van ∆ν = 10-9 Rad. Als test op de iteratiefractal is onderzocht in hoever het uiterlijk van deze fractal en dus daarmee de kwaliteit afhangt van de eindprecisie van ν. Zoals vermeld is het blauwe, en dus stabiele, gebied direct bruikbaar. In de benadering wordt uitgegaan van alle elliptische banen vanaf cirkelvormig tot een ellipticiteit van e =0.6 geldend voor het interval M=0 tot M=180° Als test-case is een lineair gemiddeld aantal iteraties genomen die in dat in dit beschreven gebied kan optreden. De test is uitgevoerd met vaste stapgrootte ∆M =1° en ∆e=0.01 gekozen. Het resultaat tussen de gewenste precisie van ν en het benodigd gemiddeld aantal iteraties is in de onderstaande grafiek weergegeven.
Figuur 2.3.2.1 Tevens is het opmerkelijk dat de opgegeven precisie nauwelijks van invloed is op een verschuiving van de grenslijn tussen het stabiele blauwe gebied en het kleurrijke chaotische gebied. Het maximaal aantal iteraties blijft tot aan die grenslijn (en dus ook tot e=0.6) minder of gelijk aan 10 ! In het chaotische gebied treedt enige verschuiving naar de rode kleur op wat duidt op een toename van het aantal iteraties bij toenemende precisie. 2.4 DE BRUIKBAARHEID VAN DE METHODE: Zoals vermeld, is de kwaliteit van de beschreven methodiek afhankelijk van het aantal iteraties en dus tevens van M en e De fractal (figuur 2.3.1) laat duidelijk de scheiding tussen de stabiele en de chaotische gebieden zien. Alleen in het blauwe (stabiele) gebied is de methode goed bruikbaar. Voor dit bruikbare gebied geldt: • Voor elke beginwaarde M over een het bereik 0 tot 180° is een eenduidige waarde voor ν te vinden. • Er zijn relatief weinig berekeningsstappen nodig; dus een gunstige snelheid van het berekeningsproces. De vraag is nu om vooraf een toetsing te kunnen maken om uit te bepalen of de beschreven methode te gebruiken is voor de te berekenen positie van het hemellichaam in haar baan. Er kan direct worden verondersteld dat de methode geschikt is voor alle elliptische banen met een e<=0.6 Dit is op zichzelf correct maar er meer uit te halen. Het komt erop neer dat het stabiele gebied alle (M,e) punten bevat die zich onder de scheidingslijn bevinden van het blauwe en het chaotische gebied. Een goede polynoombenadering van deze scheidingslijn waaraan de bruikbaarheid kan worden getoetst indien e voldoet aan : e < -0,02764⋅M5 + 0,21536⋅M4 - 0,62536⋅M3 + 0,98589⋅M2 - 0,90426⋅M + 0,98000 [2.4.1] M moet worden opgegeven in radialen. Er dient echter te worden opgemerkt dat vergelijking [2.4.1] slechts geldig is voor M gelegen in het interval tussen 0 en π . Voor het interval van π en 2π dient M te worden berekend als 2π - M
13 2.5 VOORBEELD BEREKENING Gevraagd: De ware anomalie van komeet Schoemaker-Levy 6 op 2.0 juli 1999 uit opgegeven baanelementen geldig voor het equinox J2000. Oplossing: Opgave baanelementen komeet Schoemaker-Levy 6 uit de sterrengids 1999 : Epoche = JD 2451320.5 Perihelium doorgang T = 1999 mei 2.37579 baanexcentriciteit e = 0.7052556 Perihelium q = 1.1345121 AE Omlooptijd P = 7.55 jaar of beter via 3e Wet Kepler uit e en q : P=7.5518524… jaar Als eerste dient de middelbare anomalie M berekend te worden. Deze wordt gevonden uit de lineaire tijdsverhouding. Voor M in radialen geldt (zie vergelijking [2.2.3]) : M = 2π ⋅ t / P De tijdsspanne t vanaf het perihelium : 2.37579 Mei - 2.0 Juli geeft een verschil van 60.62421 dagen ofwel t in Juliaanse jaren = 0.1659526... jaar Hieruit volgt M = 0.138074.. Rad. Echter e is groter dan 0.6 en daarom dient vergelijking [2.4.1] als toetsing of de beschreven methode te gebruiken is. Uit vergelijking [2.4.1] : -0,02764⋅M5 + 0,21536⋅M4 - 0,62536⋅M3 + 0,98589⋅M2 - 0,90426⋅M + 0,98000 = 0.87237 En daar e kleiner is dan 0.87237 voldoet de vergelijking aan de bewering en mag de verdere afleiding volgens de beschreven methodiek worden voortgezet. Vervolgens dient het aphelium Q van de elliptische baan te worden berekend. Uit substitutie van de vergelijkingen [2.1.9] 1+e Q = q⋅
hieruit volgt een aphelium Q = 6.5637654… AE 1-e
Uitgaande van vergelijking [2.2.1] dient de oppervlakte verhouding k te worden bepaald: k = Atot/P Hiervoor is het totaal oppervlak Atot van de elliptische baan nodig. Dit oppervlak is direct te berekenen uit verg [2.1.12]
π⋅q⋅Q Atot =
1-e
2
Atot = 32.99858… AE2
en hiermee wordt de k factor vastgesteld op k = 4.369594…. AE2 / jaar Het oppervlak A(t) van de ingesloten perk kan worden berekend uit de lineaire oppervlak verhouding k ⋅ t en bedraagt in ons geval k ⋅ 0.1659526 jaar = 0.725145… AE2 en met de eerste waarde νi=1 = M = 0.1380735325 Rad wordt tevens de startwaarde voor het iteratieproces verkregen.
14 Met behulp van de methode van Newton-Raphson (zie vergelijking [2.2.5]) kan de ware anomalie ν worden benaderd. De onderstaande tabel geeft het aantal iteratiestappen aan om ν te bepalen. Iteratie stap
A(νi)
dAi / dνi
A(νi) - A(t) (foutterm)
νI
νi+1
[Rad]
[Rad]
1
0.0890925883
0.6486550485
-0.6360528146
0.1380735325
1.1186452802
2
0.8687205272
1.0936221507
0.1435751242
1.1186452802
0.9873612493
3
0.7335586841
0.9706491486
0.0084132811
0.9873612493
0.9786935643
4
0.7251761432
0.9635758213
0.0000307402
0.9786935643
0.9786616620
5
0.7251454033
0.9635500054
0.0000000004
0.9786616620
0.9786616616 Tabel 2.5.1
De parameters in de overige kolommen worden berekend uit de afgeleide vergelijkingen [2.1.6] en [2.2.6]. In de laatste kolom blijkt na 5 iteratiestappen het gevraagde antwoord voor ν met een berekeningsnauwkeurigheid van ∆ν = 10-9 Rad. In de 4e kolom is duidelijk de mate van convergentie van het iteratieproces te zien, de foutterm wordt in een weinig iteratiestappen naar 0 gereduceerd. Hiermee is de berekening voldaan en het gezochte antwoord luidt: De ware anomalie ν voor komeet Schoemaker-Levy 6 op 2 juli 1999 bedraagt 0.9786616616.. Rad ofwel 56.07318..° Duidelijk is te zien dat door de sterke ellipsvorm van de baan de ware anomalie behoorlijk afwijkt van de middelbare anomalie.
15 3. ELLIPSBAAN CURIOSA 3.1 DE GEMIDDELDE VOERSTRAAL Voor het berekenen van de radiusvector (voerstraal van een hemellichaam) als functie van de ware anomalie geldt volgens vergelijking [2.1.1]: P R(ν) =
met : 1 + e⋅COS ν
P = 2⋅q⋅Q / (q+Q)
Over de gehele elliptische baan wordt doorgaans de halve lange as vaak als een soort gemiddelde op de voerstraal berekend : Rgem = a = ( q + Q ) / 2
[3.1.1]
Met deze beschreven manier voor het bepalen van de gemiddelde voerstraal R vanuit het perifocus wordt uitgegaan van de extremen van de ellipsbaan. 3.2 DE TIJDSGEMIDDELDE VOERSTRAAL Uit de perkenwet van Kepler komt duidelijk naar voren dat een hemellichaam haar elliptische baan niet eenparig als functie in de tijd doorloopt (zie hoofdstuk 2) Er kan dus ook worden afgevraagd wat de gemiddelde voerstraal in de tijd is tijdens het doorlopen van de volledige omloopperiode door de baan. Immers de perihelium passage wordt in een korter tijdsbestek doorlopen dan bij het aphelium het geval is. Deze te bepalen gemiddelde voerstraal zou bijvoorbeeld nuttig kunnen zijn om te berekenen wat de gecumuleerde Zon invloed van een planeet is gedurende haar omloop om de Zon. (of mogelijk kunnen nieuwe mathematische inzichten vereenvoudigingen betekenen van de Keplerse banen) En hopelijk kan hiervoor een mathematische relatie voor Rgem worden afgeleid met betrekking tot de vorm van de elliptische baan.. Hiervoor is het volgende experiment gedaan: De omloopperiode P is verdeeld over een groot aantal eindige tijd-elementjes met allen eenzelfde tijdsduur, te beschouwen als ∆t . De doorlooptijd t door de elliptische baan wordt daarmee gedefinieerd als: n
tn =
Σ ∆t i
[3.2.1]
i =1
En voor de totale omloopperiode geldt dus tn = P Na iedere tijd t wordt de ware anomalie ν volgens de beschreven methodiek en PASCAL routine ‘PERKHOEK’ (zie Bijlage I) berekend met evenals de daarmee samenhangende voerstraal R(ν) volgens verg. [2.1.1] Vervolgens wordt het gemiddelde van R over de gehele omloopperiode berekend als. n
Rgem
=
Σ
R(ν(t=i) )
i =1
n
[3.2.2]
16 De resultaten van het experiment zijn weergegeven in de onderstaande tabel. Hiertoe zijn voor q en Q een aantal “mooie” waarden gekozen en na de opdeling van de omlooptijd in een groot aantal (>1000 tijdselementjes) zijn opmerkelijke waarden gevonden. Dit is opmerkelijk want de berekeningen om via de ware anomalie de voerstraal R te berekenen zijn vrij complex. q
Q
1
1
1.0000000000
Rgem in breuknotatie 1
1
2
1.5833333333
19/12
1
3
2.2500000000
9/4
1
4
2.9500000000
59/20
1
5
3.6666666666
11/3
1
6
4.3928571429
123/28
1
7
5.1250000000
41/8
1
8
5.8611111111
211/36
1
9
6.6000000000
66/10
Rgem
=>
Tabel 3.2.1 Opmerkelijk is dat de getalswaarden voor Rgem een verband lijken te hebben met de vermenigvuldiging van q+Q Door empirisch uitsplitsen van de uitkomsten en door integratie terugherleiden kan inderdaad een afleiding voor Rgem gevonden (zie voor de afleiding Bijlage II): 3 x2 + 2 x + 3 Rgem = q ⋅ -------------------∆t→0 4⋅( 1 + x )
Hierin: x = Q/q [3.2.2]
Na onderzoek met een groot aantal andere (minder mooie) waarden voor q en Q leerde dat de functie, echter zonder wiskundig bewijs te leveren, voordoet voor q ,Q ∈ +ℜ (alle positieve rationele waarden van q en Q met uitgezonderd van 0) Hieronder een grafische voorstelling op schaal van een elliptische baan met een excentriciteit van 0.6:
Figuur 3.2.1 De tijdsgemiddelde voerstraal stelt de gestippelde cirkel voor. Op het moment van de getekende positie van de planeet geldt dat de tijdsintervallen T1 en T2 aan elkaar gelijk zijn. Dit moment valt
17 echter niet samen met een der snijpunten van de werkelijke ellipstische baan en de gestippelde cirkel. (dus daar waar de beide voerstralen aan elkaar gelijk zijn) Mogelijk kan, met gebruik makend van de tijdsgemiddelde voerstraal, toch een relatie worden gelegd tot de tijdsvereffening zonder gebruik te maken van de omvangrijke afleidingen die via de 2e wet van Kepler noodzakelijk zijn. Maar waar slechts vergelijking [2.1.1] kan dienen voor het terugrekenen naar het argument van ν. Een ander nut kan mogelijk blijken voor het berekenen van de gemiddelde baansnelheid van de planeet. Voor de baansnelheid op een willekeurig punt in die baan geldt immers: d( R(ν) ⋅ν) V(t) =
dt
met d(R(ν) ⋅ν) als doorlopen boogelementje v/d ellips
Het berekenen van de gemiddelde baansnelheid uit de tijdsgemiddeldevoerstraal zou dan eenvoudig zijn met: Vgem = Rgem / P De begripsvorming van de tijdsgemiddelde voerstraal houdt mij momenteel nog bezig, mogelijk blijkt er verder geen nut uit, maar het intrigeert mij in ieder geval om de wijze waarop deze mathematisch kan worden afgeleid. 3.3 DE HOEKGEMIDDELDE VOERSTRAAL Op eenzelfde wijze waarop de tijdsgemiddeldevoerstraal is berekend kan uit vergelijking [2.1.1] worden onderzocht wat de gemiddelde voerstraal is wanneer de gehele elliptische baan wordt doorlopen met constante stappen ∆ν. Ook nu ligt het in de verwachting weer ander antwoord te verkrijgen dan op de traditionele wijze om Rgem te berekenen uit de extremen q en Q. Voor de benadering van de hoekgemiddelde voerstraal geldt: n
Rgem
Σ
=
R(ν=i)
[3.3.1]
i =1
n Hierin is n de gelijke verdeling van ∆ν over het bereik van 2π De verdere afleiding is een stuk eenvoudiger dan in het geval van de tijdsgemiddelde voerstraal, welke beschreven is in 3.2 , omdat hiervoor slechts vergelijking [2.1.1] nodig is en er niet geitereerd behoeft te worden. Zonder verdere afleiding kan worden gevonden wanneer ∆ν voldoende klein is (en daarmee n groot): q⋅Q
Rgem = ∆ν → 0
( q ,Q ∈ +ℜ )
[3.3.2]
Het meest opvallende is dat blijkt te gelden : Rgem = ∆ν → 0
b
[3.3.3]
Hieruit volgt dus de onderstaande stelling: De hoekgemiddelde voerstraal over de gehele ellips, gerekend vanuit een der brandpunten, is gelijk aan de halve korte as van diezelfde ellips!
18 BIJLAGEN BIJLAGE I PROGRAMMA ROUTINES IN PASCAL (KEPLERSE BANEN) {ÉÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍ»} {º FUNCTION R_ellips(R1,R1,nu:real):real º} {º Deze functie berekent de voerstraal als functie van nu º} {º van een ellips met: º} {º R1 = apfocus afstand º} {º R2 = perifocus afstand º} {º nu = anomalie (vanaf het perifocus) º} {º De functie is geldig voor nu= 0- 2pi (nu in RAD) º}
{º {º {º {º {º
R
P = -----------1 + e.COS(ν)
º} º} º} º} º}
{ÈÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍͼ} FUNCTION R_ellips(R1,R2,nu:real):real; VAR e:real; P:real; BEGIN e:=2*R1/(R1+R2) - 1; P:=2*R1*R2/(R1+R2); R_ellips:=P/(1+e*COS(nu)) END; {ÉÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍ»} {º FUNCTION O_ellips(R1,R1,nu:real):real º} {º Deze functie berekent de oppervlak van een ellips met: º} {º R1 = apfocus afstand º} {º R2 = perifocus afstand º} {º nu = anomalie (vanaf het perifocus) º} {º De functie is geldig voor nu= 0- 2pi (nu in RAD) º}
{º {º {º {º {º
Opp =
b P2 ⌠ dν --- . -----------2 ⌡ (1+ e.COS(ν))2 a
º} º} º} º} º}
{º berek.integraal vlg. Handbook of Chemistry & Physics º} {ÈÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍͼ} FUNCTION O_ellips(R1,R2,nu:real):real; {----------------------------------------------------------} FUNCTION Opp(R1,R2,nu:real):real; VAR e,P,k :real; int1,int2:real; a,b :real; BEGIN e:=2*R1/(R1+R2) - 1; P:=2*R1*R2/(R1+R2); k:=SQRT(1-e*e); IF nu<>pi THEN BEGIN int2:=2*ARCTAN(k*tan(nu/2)/(1+e))/k; int1:=e*SIN(nu)/(1+e*COS(nu)); Opp:=ABS(R1*R2*(int1 - int2)/2) END ELSE BEGIN a:=(R1+R2)/2; b:=a*SQRT(1-e*e); Opp:=0.5*pi*a*b {opp. halve ellips} END;
END; {----------------------------------------------------------} BEGIN IF nu>=pi THEN O_ellips:=2*Opp(R1,R2,pi)-Opp(R1,R2,2*pi-nu) ELSE O_ellips:=Opp(R1,R2,nu) END; {----------------------------------------------------------}
19 {ÉÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍ»} {º FUNCTION perkhoek(R1,R2,tel,omloop:real):real; º} {º Deze functie berekent de ware anomalie van de baan met º} {º R1 = apfocus afstand º} {º R2 = perifocus afstand º} {º tel = lineaire tijdseenheid vanaf perifocus º} {º Omloop = totale tijd voor volledig doorlopen ellips º} {º Berekeningswijze vlg. 2e wet Kepler (perkenwet): º} {º º} {º Opp1 Opp2 º} {º ---- = ---º} {º t1 t2 º} {º º} {º perkhoek = ware anomalie in Rad vanaf het perifocus º} {º dus: perkhoek= f(tel) º} {ÈÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍͼ} FUNCTION perkhoek(R1,R2,tel,tomloop:real;VAR iter:integer):real; VAR kc :real; M :real; nu :real; Opptot,Opp:real; {----------------------------------------------------------} FUNCTION hoek(R1,R2,opp,nu:real):real; VAR fout:real; M :real; {..........................................................} FUNCTION integrant(R1,R2,nu:real):real; VAR e:real; P:real; BEGIN e:=2*R1/(R1+R2) - 1; P:=2*R1*R2/(R1+R2); integrant:=SQR(P)/(2*SQR(1+e*COS(nu))); END; {..........................................................} PROCEDURE New_Rap(R1,R2,M:real;VAR nu:real); VAR fout :real; BEGIN fout:=O_ellips(R1,R2,M)-opp; nu:=M-fout/integrant(R1,R2,M); END; {..........................................................} BEGIN WHILE (ABS(nu-M)>10E-9)AND(Iter<1000) DO BEGIN M:=nu; New_Rap(R1,R2,M,nu); Iter:=Iter+1; END; WHILE nu>=2*pi DO nu:=nu-2*pi;WHILE nu<0 DO nu:=nu+2*pi; hoek:=nu; END; {----------------------------------------------------------} BEGIN Iter:=0; WHILE tel>tomloop DO tel:=tel-tomloop; WHILE tel<0 DO tel:=tel+tomloop; M:=2*pi * tel/tomloop;
{begin schatting hoek nu}
Opptot:=O_ellips(R1,R2,2*pi); kc:=Opptot/tomloop; Opp:=kc*tel; PerkHoek:=hoek(R1,R2,Opp,M); END; {----------------------------------------------------------}
20 {ÉÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍ»} {º FUNCTION perktijd(R1,R2,nu,tomloop:real):real; º} {º Deze functie berekent de oppervlak van een ellips met: º} {º R1 = apfocus afstand º} {º R2 = perifocus afstand º} {º nu = ware anomalie (vanaf het perifocus) in Rad. º} {º tomloop= totale tijd voor volledig doorlopen ellips º} {º º} {º perktijd => is tijdseenheid als functie van nu º} {º dat is verstreken vanaf het moment van de º} {º laatste perifocusdoorgang º} {º De functie is geldig voor nu= 0- 2pi (nu in RAD) º} {ÈÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍÍͼ} FUNCTION perktijd(R1,R2,nu,tomloop:real):real; VAR kc :real; Opptot,Opp:real; BEGIN WHILE nu>2*pi DO nu:=nu-2*pi; WHILE nu<0 DO nu:=nu+2*pi; Opptot:=O_ellips(R1,R2,2*pi); kc:=Opptot/tomloop; Opp:=O_ellips(R1,R2,nu); perktijd:=Opp/kc
END; {==========================================================}
21 BIJLAGE II DE TIJDSGEMIDDELDE VOERSTRAAL BEREKENING AFLEIDING VOOR Rgem : Als basis voor de afleiding dient tabel 3.2.1 en zoals daar reeds is vermeld lijken de waarden van Rgem een verband te bezitten voor het product van q+Q. (ofwel in dit geval met de aanname 1+x) Hieronder volgt de uitsplitsingstabel van tabel 3.2.1: q
Q
x= Q/q
Rgem
y=(1+x)⋅Rgem
∆y=yi+1-yi
∆(∆y) = ∆yi+1-∆yi
1
1
1
1
2
2.75
1.5
1
2
2
19/12
4.75
4.25
1.5
1
3
3
9/4
9
5.75
1.5
1
4
4
59/20
14.75
7.25
1.5
1
5
5
11/3
22
8.75
1.5
1
6
6
123/28
30.75
10.25
1.5
1
7
7
41/8
41
11.75
1.5
1
8
8
211/36
52.75
13.25
1
9
9
66/10
66
Tabel Bijlage II De kolommen 1,2 en 4 zijn overgenomen uit tabel 3.2.1. Kolom 3 is de mate van ellipticiteit als breuk representatie (Q/q) en dus niet als e. In kolom 5 zijn de getransformeerde gemiddelde waarden te vinden, welke worden gevormd door dat de noemers van Rgem worden weggedeeld door de vermenigvuldiging van 1+x. In deze kolom : y = (1+x) ⋅ Rgem In de laatste 2 kolommen wordt steeds het verschil tussen de voorgaande y waarden uitgesplitst. De uitsplitsingsprocedure dient 2x te worden uitgevoerd om in de laatste kolom gelijke waarden te doen ontstaan. Dankzij de situatie dat in die laatste kolom alle waarden aan elkaar gelijk zijn is een mathematische afleiding mogelijk voor de bepaling van y=f(x) Wanneer ervan uit wordt gegaan dat de gezochte functie f(x) continue is en geldig is voor alle rationele waarden van q en Q kunnen de volgende integraties plaats vinden: Als eerste geldt dan uit de laatste kolom: d2y/dx2 = 1.5 Hieruit volgt dus voor de voorgaande kolom: ⌠ dy/dx = 1.5 dx ⌡
=
1.5⋅x + C1
En voor y geldt dan: y
=
⌠ ( 1.5⋅x + C1 )⋅dx ⌡
=
0.75⋅x2 + C1⋅x + C2
Wanneer voor y twee tabelwaarden worden ingevuld worden de integratie constanten C1 en C2 gevonden . Met deze waarden wordt gevonden: y = 0.75⋅x2 + 0.5⋅x + 0.75
22 En omdat y is gedefinieerd als (1+x)⋅Rgem geldt: 0.75 x2 + 0.5 x + 0.75 Rgem =
1+x
Echter nu is nog een enkel probleem dat Rgem een genormeerde voerstraal is , dit impliceert dat Rgem slechts afhankelijk is van de mate van ellipticiteit (x= Q/q) en dus nog niet van de grootte van de ellips. Gelukkig is het verband in ellipsgrootte evenredig met Rgem en dient daarom met de perifocus afstand te worden gecorrigeerd. Wanneer vervolgens teller en noemer met 4 wordt vermenigvuldigd om op gehele waarde te komen wordt de afleiding als volgt gecompleteerd:
Rgem
3 x2 + 2 x + 3 = q ⋅ -------------------4⋅( 1 + x )
23 REFERENTIES
1
Astronomical algorithms
J. Meeus
ISBN 0-943396-35-2
Willmann-Bell
2
Astronomical tables of the Sun, Moon…
J. Meeus
ISBN 0-943396-45-X
Willmann-Bell
3
Atlas van de Wiskunde deel 1
ISBN 902466957
Sesam
th
4
Handbook of Chemistry & Physics 78 .
David R. Lide
ISBN 0849304784
CRC
5
Methods of orbit determination
Dan Boulet
ISBN 0-943396-34-4
Willmann-Bell
6
Sterrengids 1998
M.Drummen e.a
ISBN 90-6638-032-2
De Koepel
7
Vademecum van de wiskunde
ISBN 90 274 2269 9
Prisma
