GLProgramming.com

home :: about :: development guides :: irc :: forums :: search :: paste :: links :: contribute :: code dump

-> Click here to learn how to get live help <-


New Paste :: Recent Pastes:: No Line Numbers


C++ precession of perihelion of mercury by taby2
1
 
// precession of perihelion of mercury

#include <iostream>
using std::cout;
using std::endl;

#include <cmath>

namespace
{
    const double c2 = pow(299792458.0, 2.0);
    const double pi = 3.1415926535897932384626433832795;
    const double G = 6.6742e-11;

    // mass of Sun, in kilograms
    const double M = 1.988435e30;

    // semi-major axis of the process known as Mercury's orbit, in metres
    const double a = 5.7909176e10;

    // eccentricity of Mercury's orbit
    const double e2 = pow(0.20563069, 2.0);

    // time of Mercury's orbit, in Earth days and Earth seconds
    const double t_d = 87.96934;
    const double t_s = t_d*24.0*60.0*60.0;
};

double revolutions_to_arcseconds_per_century(void);
double method1(void);
double method2(void);

int main(void)
{
    // in arcseconds of a degree per century
    cout << method1() << endl;
    cout << method2() << endl;

    return 0;
}

double revolutions_to_arcseconds_per_century(double angle)
{
    // 2pi radians per orbit
    angle *= 2.0*pi;

    // 180/pi * 60 seconds per degree * 60 arcseconds per second
    angle *= 206264.80624709635515647335733078;

    // orbits per 100 Earth years
    angle *= 365.0*100.0 / t_d;

    return angle;
}

double method1(void)
{
    // precession rate per orbit
    double angle = (3.0*G*M) / (c2*a*(1.0 - e2));

    return revolutions_to_arcseconds_per_century(angle);
}

double method2(void)
{
    // precession rate per orbit
    double angle = (12.0*pi*pi*a*a) / (t_s*t_s*c2*(1 - e2));

    return revolutions_to_arcseconds_per_century(angle);
}