Skip to content
Snippets Groups Projects
kepler.c 976 B
Newer Older
Bruce Cowan's avatar
Bruce Cowan committed
#include <stdio.h>
#include <stdlib.h>
#include <tgmath.h>

#define ACC 1e-9

static double
solve_keplers_equation (double m, double ecc)
{
    static int i = 1;
    double e = m; // Initial guess

    while (1)
    {
        printf ("---Iteration #%d---\n", i);

        double delta = e - ecc * sin (e) - m;
        printf ("E = %f, delta = %g\n", e, delta);
    
        if (fabs (delta) <= ACC)
            return e;

        e -= delta / (1 - ecc * cos (e));
        i++;
    }
}

int
main (int argc, char **argv)
{
    double m;
    double ecc;

    if (argc != 3)
    {
        printf ("Input M (degrees) ");
        scanf ("%lf", &m);
        printf ("Input e ");
        scanf ("%lf", &ecc);
    }
    else
    {
        sscanf (argv[1], "%lf", &m);
        sscanf (argv[2], "%lf", &ecc);
    }

    m *= (M_PI / 180.0);    // Convert m into radians
    double e = solve_keplers_equation (m, ecc);

    printf ("\nResult is E = %.17g rad\n", e);
    
    return 0;
}