/*
 * Note: often just:
 * (1+e2)*(sin(E) + 0.25*e*(1+e2)*sin(2*E) + e2*sin(3*E)/3.0);
 * or:
 * (1+e2)*sin(E)*((1+e2) + 0.5*e*(1+e2)*cos(E) - 4*e2*sin^2(E)/3.0);
 * suffices.
 */
double TrueAnomalySeries(double E, double e)
{
    double e2 = 0.25*e*e, beta, t, v;
    //beta = 2*(1.0 - sqrt(1-e*e))/(e*e);
    //beta = (1.0+e2*(1.0+2*e2*(1.0+2.5*e2)));
    beta = (1.0+e2*(1.0+2*e2*(1.0+e2*(2.5+7.0*e2)));
    t = e*beta;
    //v = beta*(sin(E) + 0.25*t*(sin(2*E) + t*(sin(3*E)/3.0 + 0.125*t*sin(4*E))));
    //v = beta*(sin(E) + 0.25*t*(sin(2*E) + t*(sin(3*E)/3.0 + 0.25*t*       (0.5*sin(4*E) + 0.2*t       *sin(5*E)))));
    v = beta*(sin(E) + 0.25*t*(sin(2*E) + t*(sin(3*E)/3.0 + 0.25*e*(1+e2)*(0.5*sin(4*E) + 0.2*e*(1+e2)*sin(5*E)))));
    return e*v;
}