BlackHolePerturbationToolkit / KerrGeodesics

Code to compute quantities related to bound timelike Kerr geodesics
MIT License
27 stars 13 forks source link

Negative values of a #22

Closed barrywardell closed 1 year ago

barrywardell commented 2 years ago

Most functions currently allow negative values for the spin parameter, a, but do not necessarily produce correct results in that case as several places assume a is positive and retrograde orbits are specified by setting the inclination angle to be pi. This is particularly problematic when interacting with the Teukolsky package as code similar to the following silently produces incorrect results:

orbit = KerrGeoOrbit[-0.5, 10.0, 0, -1];
TeukolskyPointParticleMode[-2, 2, 2, 0, 0, orbit]
nielsw2 commented 2 years ago

We should decide if in the Toolkit we always take 0 <= a <= 1 and let x decide between prograde and retrograde, or if we want to allow both positive and negative a and also let x choose between prograde and negative. I suggest we initially pick the former and the consider if the it important to let a vary smoothly between prograde and retrograde (for fixed x).

barrywardell commented 2 years ago

We should also keep in mind packages other than KerrGeodesics. In particular, the Teukolsky package doesn't necessarily know anything about geodesics and prograde vs retrograde orbits, but there is still the possibility of having -1<a<0.

nielsw2 commented 2 years ago

At the ICERM workshop we have decided that 0 <= a <=1 and the sign of x will differentiate between prograde and retrograde orbits. All functions should not allow a < 0 and if the user inputs that then an error should be returned that explains the convention.

MvdMeent commented 2 years ago

I strongly disagree with that choice, as it is user unfriendly. The (a,p,e,x) parametrization of orbital configuration space naturally behaves as a double cover. Functions should be implemented such that they correctly reflect that.

znasipak commented 2 years ago

I agree with @MvdMeent. For people studying equatorial Kerr orbits, allowing $a$ to be positive and negative provides a smooth way of transitioning from prograde to retrograde orbits.

barrywardell commented 2 years ago

The problem is that the code currently silently produces wrong results for negative a. For example,

Plot[{KerrGeoAngularMomentum[a, 10, 0.1, 1], Sign[a] KerrGeoAngularMomentum[Abs[a], 10, 0.1, Sign[a]]}, {a, -0.9, 0.9}]

produces

KerrGeodesicsProblem

@znasipak is going to work on trying to make everything correct for $a<0$, but if that ends up being very difficult I propose we initially disallow negative $a$.

MvdMeent commented 2 years ago

The code currently also regularly gives wrong results for x=-1.

E.g.

ListPointPlot3D@
 Flatten[Table[{a, x, 
    KerrGeoConstantsOfMotion[a, 10.0`32, 0, x][
     "\[ScriptCapitalL]"]}, {a, -1, 1, 1/10}, {x, -1, 1, 1/10}], 1]

produces: image

MvdMeent commented 2 years ago

I think the correct thing to do, is to make negatives a's work when we can, and print a warning message for functions that have not been tested for negative a.

MvdMeent commented 2 years ago

Some of the bad behaviour can be fixed with

KerrGeoEnergy[a_?Negative,p_,e_,x_]:= KerrGeoEnergy[-a,p,e,-x]

znasipak commented 2 years ago

I'm currently working on an implementation that will allow for negative $a$ and return expected amplitudes and fluxes when combined with the Teukolsky package

nielsw2 commented 1 year ago

Just looking to see if this can be resolved before Capra. It looks like it probably cannot, unless Zach you implemented some fixes? (I don't see a pull request).

There also seem to be other issues with retrograde orbits that we've not discussed, e.g., the following should give the same result but does not

KerrGeoFrequencies[0.99, 10, 0, -1]
KerrGeoFrequencies[-0.99, 10, 0, 1]

Also, shouldn't the azimuthal orbital frequency be negative for retrograde orbits?

znasipak commented 1 year ago

It looks like this is now fixed by Maarten's recent pull request #48

nielsw2 commented 1 year ago

As Zach noted, this has been fixed with pull request #48 which has now been merged into the master branch.