On 15/06/15 00:46, Tomasz Buchert wrote: > [...] > Could it be a compiler bug that miscompiled the loop for > 0.13.3-2?
I take it back, now I think it is not compiler issue (see below). I extracted the parameters that make the loop infinite and made a minimal program that shows the problem (attached). As you can see, it converges quite rapidly, but finally |E-Ep| stays bigger than 1e-10 and E oscilates between 2 values. An engineer's way to solve it would be to let the loop run for, say, 100 iterations, and if it didn't converge, just show a warning. Cheers, Tomasz
#include <stdio.h>
#include <cmath>
#define EPSILON 1e-10
template <typename T> int sign(T val)
{
return (T(0) < val) - (val < T(0));
}
static void InitEll(const double q, const double n, const double e, const double dt, double &rCosNu, double &rSinNu)
{
const double a = q/(1.0-e); // semimajor axis
double M = fmod(n*dt,2*M_PI); // Mean Anomaly
if (M < 0.0) M += 2.0*M_PI;
// double E = M;
// for (;;) { // Newton(?) Solve Kepler's equation (similar to Meeus second method, Astro.Alg. 1998 p.199)
// const double Ep = E;
// E -= (M-E+e*sin(E))/(e*cos(E)-1);
// if (fabs(E-Ep) < EPSILON) break;
// }
// GZ: Comet orbits are quite often near-parabolic, where this may still only converge slowly.
// Better always use Laguerre-Conway. See Heafner, Ch. 5.3
double E=M+0.85*e*sign(sin(M));
for (;;)
{
const double Ep=E;
const double f2=e*sin(E);
const double f=E-f2-M;
const double f1=1.0-e*cos(E);
E+= (-5.0*f)/(f1+sign(f1)*std::sqrt(fabs(16.0*f1*f1-20.0*f*f2)));
printf("%.15lf %.15lf\n", E, fabs(E-Ep));
if (fabs(E-Ep) < EPSILON) break;
}
// Note: q=a*(1-e)
const double h1 = q*std::sqrt((1.0+e)/(1.0-e)); // elsewhere: a sqrt(1-e²) ... q / (1-e) sqrt( (1+e)(1-e)) = q sqrt((1+e)/(1-e))
rCosNu = a*(cos(E)-e);
rSinNu = h1*sin(E);
}
int main() {
double q = 6.2417020000000001;
double n = 2.0430262950080074e-11;
double e = 0.99999300000000002;
double dt = -540.44803517358378;
double c = -4.7238161603831923;
double s = -2.8365257674233852;
InitEll(q, n, e, dt, c, s);
return 0;
}
signature.asc
Description: Digital signature
