Cancellation in the calculation of the orbital frequencies

[nielsw2]

The equations for the Mino time t and phi frequencies (Eq. 21 of https://arxiv.org/pdf/0906.1420.pdf) have the term $r3-r+$ in the denominator and there are points in the parameter space where this is zero. As an example, consider:

a0 = 0.9;
p0 = 7.035330;
e0 = 0.251521;

rp = M + Sqrt[M^2 - a^2];

orbit = KerrGeoOrbit[a0, p0, e0, 1];

orbit["RadialRoots"] - Block[{M = 1, a = a0}, rp]

which gives

{7.96361, 4.18553, 3.41949*10^-14, -1.43589}

This does not seem to cause an issue in Mathematica (though no doubt careful fine tuning would lead to a loss of precision) but it is causing a problems with these equations are CForm'ed for use in FEW. There must be a way to re-write the equations to factor out the problematic term.

[MvdMeent]

For the record this happens exactly when L = 2 r_p E/a.

[lorenzsp]

I isolated the problematic terms and checked if the numerator rounds off to zero every time rp=r3. Fortunately, this is the case for the 10000 evolutions I tested. The C++ code was rewritten in the following way:

// This term is zero when r3 - rp == 0.0
    double prob1 = (2 * M * En * rp - a * L) * (EllipticK(pow(kr, 2)) - (r2 - r3)/(r2 - rp) * EllipticPi(hp, pow(kr, 2)));
    if (abs(prob1)!=0.0){
        prob1 = prob1/(r3 - rp);
    }
    double CapitalUpsilonPhi = (CapitalUpsilonTheta)/(sqrt(Epsilon0zp)) + (2 * a * CapitalUpsilonr)/(M_PI * (rp - rm) * sqrt((1 - pow(En, 2)) * (r1 - r3) * (r2 - r4))) * ( prob1 - (2 * M * En * rm - a * L)/(r3 - rm) * (EllipticK(pow(kr, 2)) - (r2 - r3)/(r2 - rm) * EllipticPi(hm, pow(kr,2))) );

    // This term is zero when r3 - rp == 0.0
    double prob2 = ((4 * pow(M, 2) * En - a * L) * rp - 2 * M * pow(a, 2) * En) * (EllipticK(pow(kr, 2)) - (r2 - r3)/(r2 - rp) * EllipticPi(hp, pow(kr, 2)));
    if (abs(prob2)!=0.0){
        prob2 = prob2/(r3 - rp);
    }
    double CapitalGamma = 4 * pow(M, 2) * En+ (2 * CapitalUpsilonr)/(M_PI * sqrt((1 - pow(En, 2)) * (r1 - r3) * (r2 - r4))) * (En/2 * ((r3 * (r1 + r2 + r3) - r1 * r2) * EllipticK(pow(kr, 2)) + (r2 - r3) * (r1 + r2 + r3 + r4) * EllipticPi(hr,pow(kr, 2)) + (r1 - r3) * (r2 - r4) * EllipticE(pow(kr, 2))) + 2 * M * En * (r3 * EllipticK(pow(kr, 2)) + (r2 - r3) * EllipticPi(hr,pow(kr, 2))) + (2* M)/(rp - rm) * (
        prob2
    - ((4 * pow(M, 2) * En - a * L) * rm - 2 * M * pow(a, 2) * En)/(r3 - rm) * (EllipticK(pow(kr, 2)) - (r2 - r3)/(r2 - rm) * EllipticPi(hm,pow(kr, 2)))
    )
    );

[barrywardell]

Is this related to issue #50?