Complete overhaul

Project.toml

@@ -7,6 +7,7 @@ version = "0.10.6"
 DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 FastDifferentiation = "eb9bf01b-bf85-4b60-bf87-ee5de06c00be"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+IrrationalConstants = "92d709cd-6900-40b7-9082-c6be49f344b6"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
@@ -26,9 +27,9 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
-[extensions]
-PostNewtonianForwardDiffExt = ["ForwardDiff"]
-PostNewtonianSymbolicsExt = ["Symbolics", "SymbolicUtils"]
+# [extensions]
+# PostNewtonianForwardDiffExt = ["ForwardDiff"]
+# PostNewtonianSymbolicsExt = ["Symbolics", "SymbolicUtils"]
 
 [compat]
 Aqua = "0.8"
@@ -40,6 +41,7 @@ FileIO = "1.16.3"
 ForwardDiff = "0.10.38, 1"
 HDF5 = "0.17.2"
 InteractiveUtils = "1"
+IrrationalConstants = "0.2.4"
 LinearAlgebra = "1"
 Logging = "1"
 MacroTools = "0.5.10"

docs/Project.toml

@@ -1,5 +1,7 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
+DocumenterMermaid = "a078cd44-4d9c-4618-b545-3ab9d77f9177"
 LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589"
 PlotlyJS = "f0f68f2c-4968-5e81-91da-67840de0976a"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"

docs/make.jl

@@ -1,6 +1,7 @@
 using PostNewtonian
 using Symbolics  # To document the extension
 using Documenter
+using DocumenterMermaid
 
 DocMeta.setdocmeta!(
     PostNewtonian, :DocTestSetup, :(using PostNewtonian); recursive=true, warn=false
@@ -33,6 +34,7 @@ makedocs(;
             "GWFrames" => "interface/gwframes.md",
         ],
         "Internals" => [
+            "Code diagram" => "internals/code_diagram.md",
             "Code structure" => "internals/code_structure.md",
             "internals/pn_systems.md",
             "internals/fundamental_variables.md",

docs/src/interface/high-level.md

@@ -53,9 +53,9 @@ using Plots  # Requires also installing `Plots` in your project
 plotlyjs()  # hide
 default(size=(800,480), linewidth=2, leg=:top)  # hide
 
-plot(inspiral.t, inspiral[:χ⃗₁ˣ], label=raw"$\vec{\chi}_1^x$")
-plot!(inspiral.t, inspiral[:χ⃗₁ʸ], label=raw"$\vec{\chi}_1^y$")
-plot!(inspiral.t, inspiral[:χ⃗₁ᶻ], label=raw"$\vec{\chi}_1^z$")
+plot(inspiral.t, inspiral[:χ₁ˣ], label=raw"$\vec{\chi}_1^x$")
+plot!(inspiral.t, inspiral[:χ₁ʸ], label=raw"$\vec{\chi}_1^y$")
+plot!(inspiral.t, inspiral[:χ₁ᶻ], label=raw"$\vec{\chi}_1^z$")
 plot!(xlabel="Time (𝑀)", ylabel="Dimensionless spin components")
 savefig("inspiral_spins.html"); nothing  # hide
 ```
@@ -78,7 +78,7 @@ They can be accessed by their symbols, like the spins above, or by their number
 in this list.  To access the `i`th variable at time step `j`, use `sol[i, j]`.
 You can also use colons: `sol[i, :]` is a vector of the `i`th variable at all
 times, and `sol[:, j]` is a vector of all the data at time step `j`.  For
-example, `inspiral[:χ⃗₁ˣ]` could also be written as `inspiral[3, :]`.
+example, `inspiral[:χ₁ˣ]` could also be written as `inspiral[3, :]`.
 
 By default, the output of `orbital_evolution` is just the time steps to which
 the adaptive ODE integrator happened to step.  If you know that you want the

docs/src/interface/symbolics.md

@@ -23,7 +23,7 @@ symbols as usual:
 julia> using PostNewtonian, Symbolics
 
 julia> symbolic_pnsystem = PostNewtonian.SymbolicPNSystem()
-PostNewtonian.SymbolicPNSystem{Vector{Num}, 9223372036854775805//2, Num}(Num[M₁, M₂, χ⃗₁ˣ, χ⃗₁ʸ, χ⃗₁ᶻ, χ⃗₂ˣ, χ⃗₂ʸ, χ⃗₂ᶻ, Rʷ, Rˣ, Rʸ, Rᶻ, v, Φ], Λ₁, Λ₂)
+PostNewtonian.SymbolicPNSystem{Vector{Num}, 9223372036854775805//2, Num}(Num[M₁, M₂, χ₁ˣ, χ₁ʸ, χ₁ᶻ, χ₂ˣ, χ₂ʸ, χ₂ᶻ, Rʷ, Rˣ, Rʸ, Rᶻ, v, Φ], Λ₁, Λ₂)
 
 julia> v = PostNewtonian.v(symbolic_pnsystem)
 v

docs/src/internals/code_diagram.md

@@ -0,0 +1,167 @@
+## Possibilities for Code Structure
+
+### `StaticArrays.FieldVector`
+
+This is a nice idea.  Given a `struct` with a series of fields, making
+it a subtype of `FieldVector` makes it an `SVector` or `MVector`
+(depending on whether it is immutable or mutable).  The fields can be
+accessed by the names you give them, but everything else is like a
+nicely completed `StaticArray`, which is handled well by all the SciML
+packages.  So, for example, we would define `BBH` as
+
+```julia
+using StaticArrays: FieldVector
+
+struct BBH{T, PNOrder} <: FieldVector{14, T}
+    M₁::T
+    M₂::T
+    χ⃗₁::T
+    χ₁ˣ::T
+    χ₁ʸ::T
+    χ₁ᶻ::T
+    χ⃗₂::T
+    χ₂ˣ::T
+    χ₂ʸ::T
+    χ₂ᶻ::T
+    R::T
+    Rʷ::T
+    Rˣ::T
+    Rʸ::T
+    Rᶻ::T
+    v::T
+    Φ::T
+end
+```
+
+The big disadvantage is that it is not possible to contain non-state
+variables in the `struct` — for example, the tidal coupling
+parameters.  Making them into state variables would make the ODE
+system larger than it needs to be.  Making them into type parameters
+would be bad design, lengthening compile times.  I don't see any other
+way to incorporate them into the `struct`, so this option is
+discarded.
+
+### `LabelledArrays`
+
+This is a nice design, and works well with SciML.  But it doesn't
+allow any other metadata — even `PNOrder`, which would be possible
+with `FieldVector`.  So this is discarded.
+
+### Emulating `LabelledArrays`' use of `Syms`
+
+Part of the type of a `LabelledArray` is a `Syms` object, which is a
+`Tuple` of `Symbol`s.  I could basically reimplement this for
+`PNSystem`s, which could bring a great deal of flexibility.  For
+example, if a `PNSystem` has an `e` field, it must be an eccentric
+system.  So, in principle, we could dispatch on `Syms`.  This seems
+clunky to me.  Also, it would remove our ability to further
+specialize, unless we defined adequate abstract super types.  I don't
+have a solid argument against this, but I think it would be better to
+just define a `symbols` function for each subtype.
+
+### Make `symbols` as function of type, and emulate `LabelledArrays` features
+
+This will be a lot of work, but basically copying all the methods of
+`LabelledArrays` seems like the way to go.
+
+### Separate types for mutable and immutable systems
+
+I think it would be too much work to maintain two separate types for
+each of the `PNSystem`s, one mutable and one immutable.  Probably
+better to just store the container type as a type parameter, and test
+whether that container `ismutabletype` — for example, inside
+`setindex!`, we could test and raise a more informative error if the
+system `!ismutabletype`.
+
+## Diagram
+
+```mermaid
+flowchart TB
+  %% define each layer as its own box
+  subgraph Core["<code>core</code><br/>Building blocks of the code"]
+    <!-- a["XYZ"]
+    b["CYS"] -->
+  end
+  subgraph Systems["<code>pn_systems</code><br/>Types encoding various binaries"]
+  end
+  subgraph Literature["<code>literature</code><br/>Modules for each reference containing pieces of PN expressions"]
+  end
+  subgraph Expressions["<code>pn_expressions</code><br/>Functions for computing physical quantities"]
+  end
+  subgraph Interface["<code>interface</code><br/>High-level functions for users to call"]
+  end
+
+  %% draw the arrows in build‐order
+  Core --> Systems
+  Systems --> Literature
+  Literature --> Expressions
+  Expressions --> Interface
+```
+
+- PostNewtonian
+  - interface: high-level functions for users to call
+    - orbital_evolution
+    - waveform
+    - pn
+  - pn_expressions: functions for computing physical quantities
+    - flux
+    - energy
+    - angular_momentum
+    - precession
+    - dynamics
+    - waveforms
+  - literature: modules for each reference containing pieces of PN expressions
+    - common variables
+    - Ref1
+      - import fundamental variables
+      - import common variables
+      - define variables not in common variables
+      - write individual expressions as separate functions with `@pn_expression`
+  - pn_systems: types encoding various binaries
+    - PNSystem
+      - AbstractBBHSystem
+        - QuasicircularBBH
+        - QuasisphericalBBH
+        - EccentricNonspinningBBH
+        - BBH
+      - AbstractBHNSSystem
+        - BHNS
+      - AbstractNSNSSystem
+        - NSNS
+  - core: the building blocks of the code
+    - PNExpression
+    - PNReference
+    - PNExpansion
+    - PNTerm
+
+- `PNSystem`
+  - Define abstract subtypes:
+    - `AbstractBBHSystem`
+      - `QuasicircularBBH`
+      - `QuasisphericalBBH`
+      - `EccentricNonspinningBBH`
+
+- `PNExpression`
+- `PNExpansion`
+- `PNTerm`
+
+- literature
+  - common variables
+  - each reference
+    - directory and module named by bibtex key
+    - modules are imported directly under top-level `PostNewtonian`
+    - explicitly import any fundamental variables used in the expressions
+    - import variables that are defined the same as the common variables
+    - define variables that are not defined in the common variables
+    - write individuals
+
+- PN expressions
+  - Assemble functions for high-level quantities like flux, etc.
+  - Inside each function, simply call the expressions from the
+    literature, prepending with module name.
+
+`@pn_expression` will
+
+1. look for symbols in the expression, and for any that matches a function
+   name in the current module, add a `let` binding for that symbol to equal
+   the function called on the `pnsystem` argument.

docs/src/internals/utilities.md

@@ -90,11 +90,6 @@ and replace them with the appropriate float values.
 ```@docs
 PostNewtonian.γₑ
 PostNewtonian.ζ3
-PostNewtonian.ln2
-PostNewtonian.ln3
-PostNewtonian.ln5
-PostNewtonian.ln³╱₂
-PostNewtonian.ln⁵╱₂
 ```
 
 ## Truncated series

docs/src/references.bib

@@ -0,0 +1,34 @@
+@article{Blanchet2024,
+    title = {{Post-Newtonian} theory for gravitational waves},
+    volume = {27},
+    issn = {1433-8351},
+    url = {https://doi.org/10.1007/s41114-024-00050-z},
+    doi = {10.1007/s41114-024-00050-z},
+    number = {1},
+    journal = {Living Reviews in Relativity},
+    author = {Blanchet, Luc},
+    month = jul,
+    year = {2024},
+}
+
+@article{TaylorPoisson2008,
+    title = {Nonrotating black hole in a {post-Newtonian} tidal environment},
+    volume = {78},
+    issn = {1550-7998, 1550-2368},
+    url = {http://arxiv.org/abs/0806.3052},
+    doi = {10.1103/PhysRevD.78.084016},
+    number = {8},
+    journal = {Physical Review D},
+    author = {Taylor, Stephanne and Poisson, Eric},
+    month = oct,
+    year = {2008},
+}
+
+@misc{Trestini2025,
+    title = {The Schott term in the binding energy for compact binaries on circular orbits at fourth {post-Newtonian} order},
+    url = {http://arxiv.org/abs/2504.13245},
+    doi = {10.48550/arXiv.2504.13245},
+    author = {Trestini, David},
+    month = apr,
+    year = {2025},
+}

ext/PostNewtonianSymbolicsExt.jl

@@ -19,10 +19,9 @@ import PostNewtonian:
     causes_domain_error!,
     prepare_pn_order,
     order_index,
+    iscall,
+    isadd,
     𝓔′,
-    apply_to_first_add!,
-    flatten_add!,
-    flatten_mul!,
     pn_expression,
     pn_expansion,
     @pn_expansion,
@@ -44,13 +43,9 @@ import PostNewtonian:
     X₁,
     X₂,
     ln,
-    ln2,
-    ln3,
-    ln5,
     ζ3,
     γₑ,
     _efficient_vector
-#apply_to_first_add!, flatten_add!, pn_expression,
 using RuntimeGeneratedFunctions: init, @RuntimeGeneratedFunction
 
 init(@__MODULE__)
@@ -59,6 +54,46 @@ function _efficient_vector(N, ::Type{Symbolics.Num})
     return Symbolics.variables(string(gensym()), 1:N)
 end
 
+### Moved from src/core/utilities/misc.jl
+
+"""
+    flatten_binary!(expr, symbols)
+
+Flatten nested binary operations — that is, apply associativity repeatedly.
+"""
+function flatten_binary!(expr, symbols)
+    while iscall(expr, symbols) && any(x -> iscall(x, symbols), expr.args[2:end])
+        args = expr.args[2:end]
+        i₊ = findfirst(x -> iscall(x, symbols), args)
+        args′ = [first(symbols); args[1:(i₊ - 1)]; args[i₊].args[2:end]; args[(i₊ + 1):end]]
+        expr.args[:] = args′[1:length(expr.args)]
+        append!(expr.args, args′[(1 + length(expr.args)):end])
+    end
+    return expr
+end
+
+flatten_add!(expr) = flatten_binary!(expr, ((+), :+))
+flatten_mul!(expr) = flatten_binary!(expr, ((*), :*))
+
+"""
+    apply_to_first_add!(expr, func)
+
+Apply `func` to the first sub-expression found in a "prewalk"-traversal of `expr` that
+satisfies [`isadd`](@ref).  If `func` acts in place, so does this function.  In either case,
+the expression should be returned.
+"""
+function apply_to_first_add!(expr, func)
+    found_add = false
+    MacroTools.prewalk(expr) do x
+        if !found_add && isadd(x)
+            found_add = true
+            func(x)
+        else
+            x
+        end
+    end
+end
+
 ### Moved from src/utilities/macros.jl
 
 hold(x) = x
@@ -68,7 +103,7 @@ Symbolics.derivative(::typeof(hold), args::NTuple{1,Any}, ::Val{1}) = 1
 function unhold(expr)
     MacroTools.postwalk(expr) do x
         m = MacroTools.trymatch(:(f_(i_)), x)
-        m === nothing || m[:f] !== hold ? x : Symbol(m[:i])
+        m ≡ nothing || m[:f] !== hold ? x : Symbol(m[:i])
     end
 end
 
@@ -260,10 +295,10 @@ end
 causes_domain_error!(u̇, ::PNSystem{VT}) where {VT<:Vector{Symbolics.Num}} = false
 
 function SymbolicPNSystem(PNOrder=typemax(Int))
-    Symbolics.@variables M₁ M₂ χ⃗₁ˣ χ⃗₁ʸ χ⃗₁ᶻ χ⃗₂ˣ χ⃗₂ʸ χ⃗₂ᶻ Rʷ Rˣ Rʸ Rᶻ v Φ Λ₁ Λ₂
+    Symbolics.@variables M₁ M₂ χ₁ˣ χ₁ʸ χ₁ᶻ χ₂ˣ χ₂ʸ χ₂ᶻ Rʷ Rˣ Rʸ Rᶻ v Φ Λ₁ Λ₂
     ET = typeof(M₁)
     return SymbolicPNSystem{Vector{ET},prepare_pn_order(PNOrder),ET}(
-        [M₁, M₂, χ⃗₁ˣ, χ⃗₁ʸ, χ⃗₁ᶻ, χ⃗₂ˣ, χ⃗₂ʸ, χ⃗₂ᶻ, Rʷ, Rˣ, Rʸ, Rᶻ, v, Φ], Λ₁, Λ₂
+        [M₁, M₂, χ₁ˣ, χ₁ʸ, χ₁ᶻ, χ₂ˣ, χ₂ʸ, χ₂ᶻ, Rʷ, Rˣ, Rʸ, Rᶻ, v, Φ], Λ₁, Λ₂
     )
 end
 

sources/PNpedia/.gitignore

@@ -0,0 +1 @@
+PNpedia