Keplerbahn berechnen



  • Ich habe a, b (große/kleine Halbachse) einer Flugbahn (Ellipse) und die Umlaufzeit T gegeben.
    Nun will ich die Position eines Planeten zur Zeit t berechnen.
    Das mache ich wie folgt:

    #define GRAD(x) ((x) / (2 * PI) * 360)
    
    typedef struct _LDPOINT {
      double x;
      double y;
    } LDPOINT, *PLDPOINT;
    
    double radius(double a, double e, double v)
    {
        return a * (1 - pow(e, 2)) / (1 + e * cos(v));
    }
    
    double wahre_anomalie(double e, double E)
    {
        if (E >= 0 && E <= PI)
        {
            return acos((cos(E) - e) / (1 - e * cos(E)));
        }
        else if (E >= PI && E <= (2 * PI))
        {
            return 2 * PI - acos((cos(E) - e) / (1 - e * cos(E)));
        }
        else
        {
            return -1.;
        }
    }
    
    double exzentrische_anomalie(double M, double e, double precision)
    {
        double E, Enew;
        E = M;
    
        for(;;)
        {
            Enew = M + e * sin(E);
            if (fabs(Enew - E) < precision)
            {
                break;
            }
            E = Enew;
        }
        return Enew;
    }
    
    double mittlere_anomalie(double t, double T)
    {
        return 2 * PI * t / T;
    }
    
    double numerische_exzentritaet(double a, double b)
    {
        return sqrt(pow(a, 2) - pow(b, 2)) / a;
    }
    
    int main(void)
    {
        double a = 100.;
        double b = 75.;
        double T = 200.;
        double t;
        LDPOINT p;
    
        for (t = 0.; t <= T; t+=5)
        {
            double e = numerische_exzentritaet(a, b);
            double M = mittlere_anomalie(t, T);
            double E = exzentrische_anomalie(M, e, 0.00000001);
            double v = wahre_anomalie(e, E);
            double r = radius(a, e, v);
            p.x = r * cos(v);
            p.y = r * sin(v);
    
            printf("t: %.4Lf E: %.4Lf v: %.4Lf r: %.4Lf x: %.4Lf y: %.4Lf\n", t, GRAD(E), GRAD(v), r, p.x, p.y);
        }
    
        return 0;
    }
    

    Dabei erhalte ich folgende Werte:

    t: 0.0000 E: 0.0000 v: 0.0000 r: 33.8562 x: 33.8562 y: 0.0000
    t: 5.0000 E: 25.0405 v: 52.3879 r: 40.0731 x: 24.4571 y: 31.7444
    t: 10.0000 E: 44.6187 v: 84.5394 r: 52.9191 x: 5.0359 y: 52.6789
    t: 15.0000 E: 59.7309 v: 103.6559 r: 66.6595 x: -15.7376 y: 64.7751
    t: 20.0000 E: 72.0537 v: 116.3435 r: 79.6194 x: -35.3313 y: 71.3509
    t: 25.0000 E: 82.5803 v: 125.5923 r: 91.4584 x: -53.2301 y: 74.3720
    t: 30.0000 E: 91.8773 v: 132.8025 r: 102.1668 x: -69.4196 y: 74.9597
    t: 35.0000 E: 100.2883 v: 138.7019 r: 111.8133 x: -84.0039 y: 73.7941
    t: 40.0000 E: 108.0355 v: 143.7059 r: 120.4785 x: -97.1044 y: 71.3149
    t: 45.0000 E: 115.2708 v: 148.0698 r: 128.2366 x: -108.8335 y: 67.8225
    t: 50.0000 E: 122.1029 v: 151.9605 r: 135.1515 x: -119.2879 y: 63.5321
    t: 55.0000 E: 128.6126 v: 155.4924 r: 141.2771 x: -128.5489 y: 58.6038
    t: 60.0000 E: 134.8621 v: 158.7475 r: 146.6580 x: -136.6840 y: 53.1605
    t: 65.0000 E: 140.9007 v: 161.7865 r: 151.3312 x: -143.7492 y: 47.2999
    t: 70.0000 E: 146.7687 v: 164.6560 r: 155.3269 x: -149.7903 y: 41.1016
    t: 75.0000 E: 152.4995 v: 167.3930 r: 158.6700 x: -154.8444 y: 34.6318
    t: 80.0000 E: 158.1219 v: 170.0273 r: 161.3800 x: -158.9417 y: 27.9475
    t: 85.0000 E: 163.6612 v: 172.5844 r: 163.4726 x: -162.1053 y: 21.0987
    t: 90.0000 E: 169.1402 v: 175.0860 r: 164.9592 x: -164.3529 y: 14.1305
    t: 95.0000 E: 174.5798 v: 177.5518 r: 165.8480 x: -165.6967 y: 7.0845
    t: 100.0000 E: 180.0000 v: 180.0000 r: 166.1438 x: -166.1438 y: 0.0000
    t: 105.0000 E: 185.4202 v: 182.4482 r: 165.8480 x: -165.6967 y: -7.0845
    t: 110.0000 E: 190.8598 v: 184.9140 r: 164.9592 x: -164.3529 y: -14.1305
    t: 115.0000 E: 196.3388 v: 187.4156 r: 163.4726 x: -162.1053 y: -21.0987
    t: 120.0000 E: 201.8781 v: 189.9727 r: 161.3800 x: -158.9417 y: -27.9475
    t: 125.0000 E: 207.5005 v: 192.6070 r: 158.6700 x: -154.8444 y: -34.6318
    t: 130.0000 E: 213.2313 v: 195.3440 r: 155.3269 x: -149.7903 y: -41.1016
    t: 135.0000 E: 219.0993 v: 198.2135 r: 151.3312 x: -143.7492 y: -47.2999
    t: 140.0000 E: 225.1379 v: 201.2525 r: 146.6580 x: -136.6840 y: -53.1605
    t: 145.0000 E: 231.3874 v: 204.5076 r: 141.2771 x: -128.5489 y: -58.6038
    t: 150.0000 E: 237.8971 v: 208.0395 r: 135.1515 x: -119.2879 y: -63.5321
    t: 155.0000 E: 244.7292 v: 211.9302 r: 128.2366 x: -108.8335 y: -67.8225
    t: 160.0000 E: 251.9645 v: 216.2941 r: 120.4785 x: -97.1044 y: -71.3149
    t: 165.0000 E: 259.7117 v: 221.2981 r: 111.8133 x: -84.0039 y: -73.7941
    t: 170.0000 E: 268.1227 v: 227.1975 r: 102.1668 x: -69.4196 y: -74.9597
    t: 175.0000 E: 277.4197 v: 234.4077 r: 91.4584 x: -53.2301 y: -74.3720
    t: 180.0000 E: 287.9463 v: 243.6565 r: 79.6194 x: -35.3313 y: -71.3509
    t: 185.0000 E: 300.2691 v: 256.3441 r: 66.6595 x: -15.7376 y: -64.7751
    t: 190.0000 E: 315.3813 v: 275.4606 r: 52.9191 x: 5.0359 y: -52.6789
    t: 195.0000 E: 334.9595 v: 307.6121 r: 40.0731 x: 24.4571 y: -31.7444
    t: 200.0000 E: 360.0000 v: 360.0000 r: 33.8562 x: 33.8562 y: -0.0000
    

    http://www.rainerstumpe.de/HTML/kepler_ber.html?/HTML/body_kepler_ber.html
    Auf dieser Internetseite gibt es eine Tabelle zu einer Flugbahn mit den selben Parametern. Dort sind die Werte zwar ähnlich, doch für Rechenungenauigkeiten sind die Unterschiede zu groß. Berechne ich die Werte nicht richtig oder was stimmt nicht?

    Infos zur Keplerbahn:
    http://de.wikipedia.org/wiki/Kepler-Gleichung
    http://de.wikibooks.org/wiki/Astronomische_Berechnungen_für_Amateure/_Himmelsmechanik/_Keplergleichung

    MfG Caster



  • a) #define PI 3.14159 eventuell zu #define PI atan(1)4 ändern
    b) pow(e, 2) macht man zu "e
    e", eventuell genauer, auf jeden Fall schneller



  • Listing schrieb:

    a) #define PI 3.14159 eventuell zu #define PI atan(1)4 ändern
    b) pow(e, 2) macht man zu "e
    e", eventuell genauer, auf jeden Fall schneller

    a) #define PI weglassen und M_PI aus <cmath> benutzen



  • Die pow()s hab ich mal geändert und PI ist eh schon als M_PI define'd. Aber das ändert nichts an meinen Werten, da der Fehler wie gesagt warscheinlich NICHT auf solche Ungenauigkeiten zurückzuführen ist.



  • Dort sind die Werte zwar ähnlich, doch für Rechenungenauigkeiten sind die Unterschiede zu groß. Berechne ich die Werte nicht richtig oder was stimmt nicht?

    Ich rate einfach mal: haben die vielleicht die Relativität berücksichtigt?


Anmelden zum Antworten