# Kapitel 3_b: Planetenbewegungen (Keplerbahnen)
# © 2015 Friedrich U. Mathiak,
# mathiak@mechanik-info.de
# 
> interface(displayprecision=5): restart: with(plots): with(ScientificConstants): Digits:=20:
# Wir beschaffen uns zuerst die Definition der universellen Gravitationskonstanten , die in Maple als physikalische Konstante mit der international üblichen Bezeichnung G abgespeichert ist. 
> GetConstant(G);
 Newtonian_constant_of_gravitation, symbol = G, 



                                                           3  

                   -11                      -13           m   

   value = 6.673 10   , uncertainty = 1.0 10   , units = -----

                                                             2

                                                         kg s 
;
> Gamma:=evalf(Constant(G));
                                  -11

                          6.673 10   
;
# Weiterhin entnehmen wir der Maple-Datenbasis die Massen von Sonne und Erde
> GetConstant(mass_of_Sun); GetConstant(mass_of_Earth);
                                                    30  

       mass_of_Sun, symbol = M[Sun], value = 1.99 10  , 



                           27            

         uncertainty = 5 10  , units = kg
                                                      24  

     mass_of_Earth, symbol = M[Earth], value = 5.97 10  , 



                         21            

       uncertainty = 5 10  , units = kg
;
> Mso:=evalf(Constant(mass_of_Sun));  Mer:=evalf(Constant(mass_of_Earth));
                                  30

                           1.99 10  
                                  24

                           5.97 10  
;
# Die Masse des Mondes  befindet sich  nicht in der Maple-Datenbasis wissenschaftlicher Konstanten. Wir erweitern deshalb die vorhandene Datenbasis um diesen Wert:
> AddConstant(mass_of_Moon,symbol= M[Moon],value = 7.349*10^(22),units = kg);
> GetConstant(mass_of_Moon);  Mmo:=evalf(Constant(mass_of_Moon));
mass_of_Moon, symbol = M[Moon], 



                                  22                           

  value = 7.3490000000000000000 10  , uncertainty = undefined, 



  units = kg
                                           22

                   7.3490000000000000000 10  
;
# Für alle Planeten unseres Sonnensystems ist der Wert  ¼ = “ÅMso  eine Konstante:
> AddConstant(mu_of_Sun,symbol= mu[Sun],value = Gamma*Mso,units = m^3/s^2);
> GetConstant(mu_of_Sun);  mu:=evalf(Constant(mu_of_Sun));
                                                      20  

      mu_of_Sun, symbol = mu[Sun], value = 1.327927 10  , 



                                          3

                                         m 

        uncertainty = undefined, units = --

                                          2

                                         s 
                                    20

                         1.327927 10  
;
# Zur Normierung von Längen und Zeiten werden folgende Bezugsgrößen eingeführt:
# 1.) Die  Astronomische Einheit (Einheitenzeichen AE auch AU in [m] ) entspricht etwa der mittleren Länge a der großen Halbachse der Erdbahn um die Sonne:
> AU:=convert(1, 'units', 'astronomical_unit', 'meters'); #Astronomische Einheit in [m]
                          149597870691
;
# 2.) Das siderische Jahr Tsid  ist die Zeit in [s] für einen vollständigen Umlauf von 360° der Erde um die Sonne  
> Tsid:=evalf(2*Pi*sqrt(AU^3/mu));  #Siderisches Jahr
                                           7

                   3.1548657678288133682 10 
;
# Zur Berechnung der Bahnparameter der Planeten stellen wir eine Prozedur zur Verfügung. Dazu werden die Parameter  ¼, a und e benötigt, aus denen dann  b, A, p, Ã und T berechnet werden können. Die Berechnung der Länge L = 4 a EllipticE(e)  der Bahnellipse erfolgt mittels des vollständigen elliptischen Integrals 2. Gattung  und  kann bekanntlich durch elementare Funktionen nicht ausgedrückt werden. Für dieses Integral gilt folgende Reihenentwicklung:
> series(EllipticE(e),e,3);
                     1      1     2    / 4\

                     - Pi - - Pi e  + O\e /

                     2      8              
;
# Für kreisnahe Planetenbahnen (e << 1) erhalten wir damit die Näherung: EllipticE(e) H ((1/2)*Pi)[1-((1/2)*e)^2];
> Proc_Calc_15:=proc(a,e,mu)
> #----------------------------------------------------------------------------------------------
> #Eingabe: 
> #          a: Große Halbachse der Planetenbahn in [m]
> #          e: Nummerische Exzentrizität des Planeten
> #         mu: Gamma*Mso in [m^3/s^2]
> #Ausgabe:  
> #          b: Kleine Halbachse der Planetenbahn in [m]
> #          p: Halbparameter der Ellipse in [m]
> #          A: Flächeninhalt der Ellipse in [m^2]
> #          L: Ellipsenumfang in [m]
> #      sigma: Länge des Vektors der Bahnebenennormalen in [m^2/s] 
> #          T: Umlaufzeit in [s]
> #     lambda: Länge des Laplace-Vektors in [m^3/s^2]
> #  rmax,rmin: Apozentrumdistanz (rmax), Perizentrumdistanz (rmin) in [m]
> #  vmax,vmin: Geschwindigkeiten (Extremwerte) im Perizentrum(vmax) und Apozentrum(vmin) in [m/s]
> #-----------------------------------------------------------------------------------------------
> local sdp;
> global b,p,A,L,T,TN,sigma,Tsid,lambda,rmin,rmax,vmin,vmax;
> description "Berechnung die Keplerbahn eines Planeten im Kraftfeld der Sonne";
> b:=a*sqrt(1-e^2);
> A:=evalf(Pi*a*b);
> p:=a*(1-e^2);
> L:=evalf(4*a*EllipticE(e));
> sigma:=sqrt(p*mu);
> T:=2*A/sigma;
> TN:=T/Tsid;
> sdp:=sigma/p;
> lambda:=mu*e;
> rmax:=p/(1-e);
> rmin:=p/(1+e);
> vmax:=sdp*(1+e);
> vmin:=sdp*(1-e);
> #Ausgabe der Ergebnisse
> print(`a = `,a);
> print(`b = `,b);
> print(`e = `,e);
> print(`A = `,A);
> print(`L = `,L);
> print(`T = `,T);
> print(`Umlaufdauer in sid. Jahren:  Tsid = `,TN);
> print(`p = `,p);
> print(`sigma  = `, sigma);
> print(`lambda = `,lambda);
> print(`rmax   = `,  rmax);
> print(`rmin   = `,  rmin);
> print(`vmax   = `,  vmax);
> print(`vmin   = `,  vmin);
> print(`Mittlere Bahngeschwindigkeit = `, L/T);
> end proc:
# Parameter der Bahnellipsen  für die Bewegung der Planeten im Kraftfeld der Sonne:
> #Proc_Calc_15(0.3870*AU,0.2056,mu); #Merkur
;
> #Proc_Calc_15(0.7233*AU,0.0068,mu); #Venus
;
> Proc_Calc_15(1.0000*AU,0.0167,mu); #Erde
                                            11

                  a = , 1.495978706910000 10  
                                              11

                b = , 1.4957700856126022303 10  
                          e = , 0.0167
                                              22

                A = , 7.0297540090143397861 10  
                                              11

                L = , 9.3988560394008213000 10  
                                              7

                T = , 3.1548657678288133684 10 
  Umlaufdauer in sid. Jahren:  Tsid = , 1.0000000000000000001
                                              11

                p = , 1.4955614934084298701 10  
                                                15

             sigma  = , 4.4564520498456797062 10  
                                           18

                   lambda = , 2.21763809 10  
                                                11

             rmax   = , 1.5209615513153970000 10  
                                                11

             rmin   = , 1.4709958625046030000 10  
                vmax   = , 30295.476441775067865
                vmin   = , 29300.228174680263825
     Mittlere Bahngeschwindigkeit = , 29791.619457296713736
;
# Geschwindigkeitsverlauf des Mittelpunktes der Erde bei ihrem Umlauf um die Sonne (e = 0.0167)
> vE:=sigma/p*sqrt(1+2*0.0167*cos(phi)+0.0167^2);
                                                       (1/2)

   29797.852308227665845 (1.00027889 + 0.0334 cos(phi))     
;
> plot(vE/vmax,phi=0..2*Pi,axes=boxed,gridlines=true,title= "\nGeschwindigkeitsverlauf des Erdmittelpunktes auf seiner Bahn um die Sonne\n",titlefont=["ARIAL", 15],labeldirections=["horizontal","vertical"],axesfont = ["HELVETICA","ROMAN",8],labels=["Wahre Anomalie phi", "v/vmax"],labelfont = ["HELVETICA", 12]);

> 
;
> #Proc_Calc_15(1.5237*AU,0.0934,mu); #Parameter der Bahnellipsen für die Bewegung der Planeten Mars
;
> #Proc_Calc_15(5.2030*AU,0.0484,mu); #Parameter der Bahnellipsen für die Bewegung der Planeten Jupiter
;
> #Proc_Calc_15(9.5826*AU,0.0565,mu); #Parameter der Bahnellipsen für die Bewegung der Planeten Saturn
;
> #Proc_Calc_15(19.201*AU,0.0472,mu); #Parameter der Bahnellipsen für die Bewegung der Planeten Uranus
;
> #Proc_Calc_15(30.047*AU,0.0113,mu); #Parameter der Bahnellipsen für die Bewegung der Planeten Neptun
;
> #Proc_Calc_15(39.482*AU,0.2488,mu); #Parameter der Bahnellipsen für die Bewegung der Planeten Pluto
;
# 
# Lösung der Keplerschen Zeitgleichung
> Proc_Calc_16:=proc(a,e,M,flag::integer)
> #---------------------------------------------------
> #Eingabe:
> #          a: Große Halbachse
> #          e: Nummerische Exzentrizität des Planeten
> #          M: Mittlere Anomalie [M = 2*Pi*(t-t0)/T]
> #       flag:
> #             flag = 0: Eingabe von M im Bogenmaß
> #             flag = 1: Eingabe von M im  Gradmaß 
> #Ausgabe:
> #         EA: Exzentrische Anomalie im Bogenmaß
> #        phi: Wahre Anomalie im Bogenmaß
> #          r: Abstand Planet-Sonne in [m]
> #----------------------------------------------------
> local NN,cf,sf,MR,p,sigma,sdp;
> global EA,phi,r,mu,vr,vphi,v,x1,x2;
> description "Lösung der Keplerschen Zeitgleichung";
> if flag = 1 then
>   MR:=evalf(M*Pi/180.);
>     else
>   MR:=M;
> end if;
> EA   :=fsolve(E-e*sin(E)-MR,E=MR);
> NN   :=evalf(1-e*cos(EA));
> cf   :=(cos(EA)-e)/NN;
> sf   :=sqrt(1-e^2)*sin(EA)/NN;
> phi  :=arctan(sf,cf);
> p    :=a*(1-e^2);
> r    :=p/(1-e*cos(phi));
> x1   :=r*cos(phi);
> x2   :=r*sin(phi);
> sigma:=sqrt(p*mu);
> sdp  :=sigma/p;
> vr   :=sdp*e*sin(phi);
> vphi :=sdp*(1+e*cos(phi));
> v    :=sqrt(vr^2+vphi^2);
> #Ausgabe der Ergebnisse
> print(`Wahre Anomalie     phi   = `, phi);
> print(`Exzentrische Anomalie E  = `,  EA);
> print(`Mittlere Anomalie M  = `,  MR);
> print(`Abstand Planet-Sonne  r  = `,   r);
> print(`Geschwindigkeit vr       = `,  vr);
> print(`Geschwindigkeit vphi     = `,vphi);
> print(`Geschwindigkeit v        = `,   v);
> print(`Kart. Koordinaten [x1, x2] des Planeten = `, [x1,x2]);
> end proc:
# Die Erde erreichte das Perihel ihrer Bahn am 3. Januar 2014 um 16 h MEZ. Damit ist t0 = (2Å24+16)Å60Å60  [s]
> t0:=(2*24+16)*3600.;
                                    5

                          2.30400 10 
;
# Im selben Jahr ist der Frühlingsanfang am 20. März 17:57 MEZ. Damit ist
> t:=((31+28+19)*24+17)*3600 + 57*60;
                            6803820
;
# Für die Mittlere Anomalie errechnen wir
> M:=evalf(2*Pi*(t-t0)/Tsid); evalf(convert(%,degrees));
                     1.3091528769017830382
                 75.008934583881962125 degrees
;
# Wir lösen die Keplersche Zeitgleichung für diesen Wert und erhalten
> Proc_Calc_16(1.0*AU,0.0167,M,0);
       Wahre Anomalie     phi   = , 1.3415855760862478102
       Exzentrische Anomalie E  = , 1.3253523703018652431
         Mittlere Anomalie M  = , 1.3091528769017830382
                                                         11

    Abstand Planet-Sonne  r  = , 1.5012578513137739580 10  
       Geschwindigkeit vr       = , 484.60928331944879908
       Geschwindigkeit vphi     = , 29910.916983566220655
       Geschwindigkeit v        = , 29914.842485884320131
   Kart. Koordinaten [x1, x2] des Planeten = , 



     [                        10                          11]

     [3.4109927576910705933 10  , 1.4619939877444914071 10  ]
;
> rPlanet:=unapply(a*(1-e^2)/(1 + e*cos(theta)),a,e,theta);
                      /     2\   

                    a \1 - e /   

(a, e, theta) -> ----------------

                 1 + e cos(theta)
;
> rErde:=unapply(rPlanet(1.0000,0.0167,theta),theta);
            0.999721110000    

theta -> ---------------------

         1 + 0.0167 cos(theta)
;
> rMars:=unapply(rPlanet(1.5237,0.0934,theta),theta);
            1.510407911628    

theta -> ---------------------

         1 + 0.0934 cos(theta)
;
> plot([rErde,rMars],coords = polar,gridlines=true,axes=boxed,title = "\nPlanetenbahnen\n", legend = ["Erde","Mars"], titlefont = ["ARIAL", 15], labels = ["x1-Werte in [AU]", "x2-Werte in [AU]"], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10],linestyle = [solid, solid],axesfont = ["HELVETICA","ROMAN",8],legendstyle = [font = ["HELVETICA", 9], location = bottom],scaling=constrained,color=[blue,red]);

# Für spätere Anwendungen (Hintergrundbild) erzeugen wir eine ähnliche Grafik, allerdings ohne  Titel und Legende.
> p1:=plot([rErde,rMars],coords = polar,gridlines=true,scaling=constrained,color=[blue,red]):
> T:= 'T':
# Wir animieren die Planetenbewegung für Erde und Mars, wobei wir aufgrund der kreisnahen Bahnen die im Buch TM III vorgestellten Näherungsgleichungen für den Radius r und die nummerische Exzentrizität E verwenden. Damit sparen wir uns die Lösung der nichtlinearen Keplerschen Zeitgleichung. Wir beginnen die Zeitzählung am 1. Januar 2013 um 0:00 Uhr MEZ. 
# Die Erde erreicht das Perihel ihrer Bahn am         3. Januar 2013 etwa um 7:00 Uhr MEZ. 
# Der Mars  erreicht das Perihel seiner  Bahn am  24. Januar 2013 etwa um  9:00 Uhr MEZ. 
# Damit ergeben sich folgende Periheldurchgangszeiten 
> t0E:=(2*24+7)*3600.;     #Erde
                                    5

                          1.98000 10 
;
> t0M:=(23*24+9)*3600.;    #Mars
                                     6

                          2.019600 10 
;
# Wir führen mit t = ÄÅTsid  die dimensionslose Zeit Ä ein. Darin bezeichnet  Tsid  die Umlaufzeit der Erde um die Sonne, womit sich die folgenden dimensionslosen Periheldurchgangszeiten ergeben:
> tau0E:=t0E/Tsid;tau0M:=t0M/Tsid;
                    0.0062760197920009795768
                    0.064015401878409991683
;
# Die mittlere Exzentrizität M erscheint dann in der Form:
> M:=unapply(2.*Pi*Tsid/T*(tau-tau0),T,tau,tau0);
                                          8             

                  1.9822606238545844675 10  (tau - tau0)

(T, tau, tau0) -> --------------------------------------

                                    T                   
;
> rho:=unapply(alpha*(1-e*cos(M(T,tau,tau0))),alpha,e,T,tau,tau0);
                                        / 

                                        | 

      (alpha, e, T, tau, tau0) -> alpha |1

                                        \ 



                /                        8             \\

                |1.9822606238545844675 10  (tau - tau0)||

         - e cos|--------------------------------------||

                \                  T                   //
;
> phi:=unapply(M(T,tau,tau0) + 2.*e*sin(M(T,tau,tau0)),e,T,tau,tau0);
                                               8             

                       1.9822606238545844675 10  (tau - tau0)

  (e, T, tau, tau0) -> --------------------------------------

                                         T                   



               /                        8             \

               |1.9822606238545844675 10  (tau - tau0)|

     + 2. e sin|--------------------------------------|

               \                  T                   /
;
# Speziell für Erde und Mars erhalten wir die bezogenen Größen:
> rE:=rho(1.,0.0167,Tsid,tau,tau0E);
           1.



              - 0.0167 cos(6.2831853071795864770 tau



              - 0.039433395344668839290)
;
> rM:=rho(1.5237,0.0934,1.88083*Tsid,tau,tau0M);
       1.5237 - 0.14231358 cos(3.3406449850223499608 tau



          - 0.21385273124930066021)
;
> phiE:=phi(0.0167,1.0*Tsid,tau,tau0E);
      6.2831853071795864770 tau - 0.039433395344668839290



         + 0.0334 sin(6.2831853071795864770 tau



         - 0.039433395344668839290)
;
> phiM:=phi(0.0934,1.88083*Tsid,tau,tau0M);
       3.3406449850223499608 tau - 0.21385273124930066021



          + 0.1868 sin(3.3406449850223499608 tau



          - 0.21385273124930066021)
;
# Es folgt eine Animation der Planetenbewegung  mit der Erde (blau)  und dem Mars (rot) für eine volle Umlaufzeit ( Ä = 1) der Erde. Die Sonne steht im Ursprung des Koordinatensystems.
> animate(display,[('pointplot')([rE, phiE],coords = polar,symbol=solidcircle,symbolsize = 30, colour = blue), ('pointplot')([rM, phiM],coords = polar,symbol=solidcircle, symbolsize = 20,color = red)],tau=0..1,frames=100,background = p1);

> 
;
