PR for first comments from DTS

[tomkimpson]

This PR deals with the JOSS review comments .

• Update package name to be less general
• Reorganise code to enable the user the easily select the background metric
• Update metric.jl to be less messy
• Clear up use of e.g. RelativisticDynamics.delta vs delta along with useful_functions.jl
• Use GR convention to distinguish Levi-Civita tensor vs symbol
• Update definition of Rtensor for less wasteful allocations
• Update/remove schwarzchild_covariant_riemann function which is just a specialised version of riemann
• Update the documentation in orbit.jl to refer to the correct parameters file
• Rewrite constants so that only the $GM_{\odot}$ product is used, not $G$, $M$ separatley. Update definition of $c$ to be exact.
• General updates to clarity and organisation of code
• Update documentation/article to discuss SSCs more. Ensure any comments about SSCs are clarified to be either about SSCs generally, or for a particular SSC.
• Discuss the choice to formulate as an ODE integration problem rather than action-angle.

[codecov-commenter]

Codecov Report

Merging #41 (c6c443e) into main (8385e65) will not change coverage.
The diff coverage is 100.00%.

📣 This organization is not using Codecov’s GitHub App Integration. We recommend you install it so Codecov can continue to function properly for your repositories. Learn more

@@            Coverage Diff            @@
##              main       #41   +/-   ##
=========================================
  Coverage   100.00%   100.00%           
=========================================
  Files            7         6    -1     
  Lines          419       330   -89     
=========================================
- Hits           419       330   -89     

Impacted Files	Coverage Δ
src/initial_conditions.jl	100.00% <100.00%> (ø)
src/metric.jl	100.00% <100.00%> (ø)
src/orbit.jl	100.00% <100.00%> (ø)
src/timestepping.jl	100.00% <100.00%> (ø)
src/universal_constants.jl	100.00% <100.00%> (ø)
src/useful_functions.jl	100.00% <100.00%> (ø)

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

[tomkimpson]

@duetosymmetry Thanks again for all your input. I have addressed your comments as follows:

• Update package name to be less general

It would be good to agree an alternative name before I implement the change. The subtitle of the article is "Relativistic Spin-Orbital Dynamics in Julia" which I think makes the application of the package more clear. Perhaps a rename to just RelativisticSpinOrbitalDynamics.jl? @dfm is there any issue on the JOSS side if I update the github repo name?

---

• Reorganise code to enable the user the easily select the background metric
• Update metric.jl to be less messy

The module metric.jl has been re-structured as follows.

Firstly, the metric_gAB functions have all been removed - these were only there previously to allow AD with respect to the metric during unit tests. Previously we were using Zygote.jl which requires pure functions (hence the metric_gAB) and avoidance of any array mutation We have now switched over to Enzyme.jl (which it looks like is now becoming the preferred package) and can operate directly on mutable arrays.

Secondly, all functions which are specific to the (Kerr) metric are now contained within metric.jl. For instance, the mapping from Keplerian orbital parameters to E,L,Q constants was originally in initial_conditions.jl - it has now been moved to metric.jl since it is a construction specific to Kerr. This makes it much easier to swap in/out alternative metrics in the future.

---

• Clear up use of e.g. RelativisticDynamics.delta vs delta along with useful_functions.jl

We have cleaned this up to just use e.g. delta() throughout the src/ code.

The original idea of the useful_functions.jl module was to provide a module which housed functions that were either shared between other modules or whose purpose was slightly separate from the module from which they were called. For instance, the calculation of the Levi-civita tensor seems distinct from the explicit expressions for the ODEs, and one can easily imagine this function being used elsewhere. I have often seen utility function type modules like this in many codes.

However, this module has now been significantly paired back, with some functions moved only to the test/ area (e.g. kretschman scalar) and some (e.g. delta) moved into metric.jl since they are specific to Kerr. Some functions do remain in this module, for example calculating the Levi-civita tensor or the covariant spin tensor - since these functions serve a particular purpose but are more general and somewhat distinct from the modules that call them.

As I say, having a utility functions module is quite common and it provides an area to add any future general functions that may be shared between modules, but I don't feel particularly strongly about this point; I am happy to remove and allocate the functions to the particular modules they are called from if preferred?

---

• Use GR convention to distinguish Levi-Civita tensor vs symbol

This has now been updated in useful_functions.jl. We have two functions, levi_civita_symbol() and levi_civita_tensor()

---

• Update definition of Rtensor for less wasteful allocations

This is a nice suggestion, and not one I had thought of before. However since the code is currently nowhere near being limited in speed or memory I personally don't feel it to be necessary at this time - would you agree? Whilst I agree that explicit allocations are wasteful, they are also conceptually straightforward and directly comparable with the standard form of the MPD equations. Such a decomposition is definitely, however, something that I will keep in mind for the future, especially w.r.t alternative number formats where scalings and allocations may become important.

---

• Update/remove schwarzchild_covariant_riemann function which is just a specialised version of riemann

This was only present for testing that the Riemann tensor for Kerr correctly reduces to Schwarzchild. This has been removed from the main src/ code and brought into the unit tests.

---

• Update the documentation in orbit.jl to refer to the correct parameters file.

This has been updated to the correct src/system_parameters.jl

---

• Rewrite constants so that only the $GM_{\odot}$ product is used, not $G$, $M$ separatley. Update definition of $c$ to be exact.

This was a nice idea and the constants and definitions have now been updated accordingly.

---

• General updates to clarity and organisation of code

Some updates to code organisation have already been discussed. In summary, the code has now been restructured as follows:

• metric.jl now contains all metric-specific quantities
• model.jl has been removed - the model struct is now defined in orbit.jl.
• useful_functions.jl has been paired back

---

• Update documentation/article to discuss SSCs more. Ensure any comments about SSCs are clarified to be either about SSCs generally, or for a particular SSC.

In the summary section of the article we now mention that the choice of SSC is equivalent to choosing the centre of mass and that we specifically take the TD condition for this work.

I am cautious about spending too much time in the article itself discussing choices and nuances of SSCs since it is quite involved and distracts from the main message. I have therefore relegated the discussion of this point to the docs page and pointed readers to my go-to reference of Costa & Natário, 2015. I can of course expand more in the article itself if preferred.

---

• Discuss the choice to formulate as an ODE integration problem rather than action-angle.

Formulation as an AA is interesting - thanks @duetosymmetry for bringing this up as an idea!

Ideally, looking long term I would like this package to be background independent. As you point out AA is only for a select few specific spacetimes. However, I think it could be worthwhile in the future exploring how to setup the code to use AA for those spacetimes where it is available, and so take advantage of long term precision.

On the other hand, with AA one is often primarily interested in exploring the fundamental frequencies of the system, rather than determining the evolution of position, spin, momentum etc. For pulsars, one is often interested in these quantities; for instance the orientation of the spin axis will influence the received pulse, the momentum will influence the Doppler shift, etc.

I could also be misunderstanding something, but it seems like action angle coordinates for the evolution with respect to a distant observer (rather than Carter-Mino time) are currently only well established for Schwarzchild (https://ui.adsabs.harvard.edu/abs/2022arXiv220311952W/abstract) and have only recently (Jan 2023!) been applied to Kerr (https://ui.adsabs.harvard.edu/abs/2023arXiv230108150K/abstract)?

A sentence on the choice not to use AA has been added to the end of the second paragraph of the Statement of Need.

---

Let me know if you have any further comments or queries. One additional thing I need to do is complete the swap over from Zygote.jl to Enzyme.jl in the testing and notebooks, but this is not quite ready yet.

Cheers!

[duetosymmetry]

Thanks, @tomkimpson, for your herculean effort! I am kind of under water this week (and possibly for the next 2, with a proposal deadline looming...) so I can only make trivial comments right now. I agree that Rtensor allocation is not important here (as compared to an NR code where one would be computing this at every grid point). Regarding the AA formulation for (nearly) geodesic orbits in Kerr: The 2023 paper you referred to above is for the inspiral, but geodesic order is much farther ahead. The fact that AA can be done for Kerr geodesics was known since Carter (1968), but strictly going to AA variables is not necessary to get all the benefits. The fact that the Carter-Mino time parameter makes it possible to solve for radial and polar motions independently became widely known after Mino (2003) but you still get many benefits without it. The analytical solutions for geodesic order were written down by Fujita and Hikida (2009), and repackaged into a very convenient form in 1906.05090. The fact that the program continues to work to O(s^1) was shown by Vojtech Witzany, also in (2019). (This and most EMRI work expands in spin, as I mentioned before, since if you want to work to O(s^2), you should also be including the quadrupole moment, etc... so it's not really worth it to work to "all orders" in spin, since it's inconsistent anyway.)

[duetosymmetry]

Forgot to comment... @farr do you have any time to take a look while I'm tied up with other work?

[tomkimpson]

All comments addressed. Merging and closing this PR - will continue discussion in PR for reviewer 2