Incorporate angles related functions initially in LinearResponse

define angles_gradient w.r.t. anomaly and export it + radius from anomaly function

Project.toml

@@ -1,6 +1,6 @@
 name = "OrbitalElements"
 uuid = "a3b07092-bde3-4843-b84f-c597d614ec7b"
-version = "2.0.0dev-2"
+version = "2.0.0dev-3"
 
 [deps]
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"

src/OrbitalElements.jl

@@ -48,7 +48,8 @@ export ae_to_frequencies_jacobian
 export αβ_from_frequencies, frequencies_from_αβ
 export αβ_from_frequencies_derivatives, frequencies_from_αβ_derivatives
 export αβ_to_frequencies_jacobian, frequencies_to_αβ_jacobian
-# angles gradient
+# Angle to radius and anomaly
+export radius_from_anomaly
 export angles_gradient
 
 # Resonant mappings structures

src/mappings/analytic/isochrone.jl

@@ -87,14 +87,14 @@ the analytic expression for Theta for the isochrone profile are given in
 Fouvry, Hamilton, Rozier, Pichon eq. (G10). It has the major advantage of 
 always being well-posed for -1<u<1
 """
-function Θ(
-    u::Float64,
+function _Θ(
+    w::Float64,
     a::Float64,
     e::Float64,
     model::AnalyticIsochrone,
     params::OrbitalParameters=OrbitalParameters()
 )
-    r = radius_from_anomaly(u, a, e, model, params)
+    r = radius_from_anomaly(w, a, e, model, params)
     rp, ra = rpra_from_ae(a, e)
     # dimensionless pericentre, apocentre, and radius
     x, xp, xa = (r, rp, ra) ./ radial_scale(model)
@@ -109,7 +109,7 @@ function Θ(
             * (sp + sa)
             / (x + xp) 
             / (x + xa) 
-            / (4 - u^2)
+            / (4 - w^2)
         )
         / (sqrt(2) * frequency_scale(model))
     )

src/mappings/analytic/plummer.jl

@@ -95,7 +95,7 @@ always being well-posed for -1<u<1
 
 @QUESTION: is `w` hard-coded as Hénon's anomaly here ?
 """
-function Θ(
+function _Θ(
     w::Float64,
     a::Float64,
     e::Float64,
@@ -150,7 +150,7 @@ function _radial_action_from_ae(
     end
     # Generic computations
     function action_integrand(w::Float64)::Float64
-        drdw_over_vrad = Θ(w, a, e, model, params)
+        drdw_over_vrad = _Θ(w, a, e, model, params)
         vrad = radial_velocity(w, a, e, model, params)
         return drdw_over_vrad * vrad^2
     end

src/mappings/frequencies.jl

@@ -150,16 +150,16 @@ function αβ_from_ae(
         return res
     end
     # Generic computations
-    # 1/α = Ω₀ / π * \int_{-1}^1 du Θ(u); β = L / π * \int_{-1}^1 du Θ(u) / [r(u)]^2
-    # using Θ calculations to compute frequencies: leans heavily on Θ from 
+    # 1/α = Ω₀ / π * \int_{-1}^1 dw Θ(w); β = L / π * \int_{-1}^1 dw Θ(w) / [r(w)]^2
+    # using Θ calculations to compute frequencies: leans heavily on _Θ from 
     # frequencies/anomaly.jl
     # @IMPROVE: EDGE could be adaptive based on circularity and small-ness of rperi
     # currently using Simpson's 1/3 rule only (hard-coded)
     # @IMPROVE: free the integration scheme, i.e., make it a parameter
-    function invαβ_integrands(u::Float64)::Tuple{Float64,Float64}
-        # push integration forward on two different quantities: Θ(u), Θ(u)/r^2(u)
-        integrand = Θ(u, a, e, model, params)
-        r = radius_from_anomaly(u, a, e, model, params)
+    function invαβ_integrands(w::Float64)::Tuple{Float64,Float64}
+        # push integration forward on two different quantities: Θ(w), Θ(w)/r^2(w)
+        integrand = _Θ(w, a, e, model, params)
+        r = radius_from_anomaly(w, a, e, model, params)
         return integrand, integrand / r^2
     end
     accum1, accum2 = _integrate_simpson(invαβ_integrands, params.NINT)

src/mappings/frequencies/anomaly.jl

@@ -1,7 +1,7 @@
 """
 
 @IMPROVE: Henon anomaly is hard-coded. Make it a function of parameters ?
-Mainly in 3 places: `Θ_edge`, `ru` and `drdu`
+Mainly in 3 places: `_Θ_edge`, `radius_from_anomaly` and `radius_from_anomaly_derivative`
 """
 
 ########################################################################
@@ -51,11 +51,12 @@ _henond4f(x::Float64)::Float64 = 0.
 #
 ########################################################################
 """
-    radius_from_anomaly(u, a, e, model[, params])
+    radius_from_anomaly(w, a, e, model[, params])
 
 mapping from anomaly to radius. Default is Henon anomaly.
 
-@IMPROVE: right now, Hénon anomaly is hard-coded.
+@IMPROVE: right now, Hénon anomaly is hard-coded. (except for AnalyticIsochrone and 
+SemiAnalyticPlummer, for which it is redefined using specific anomaly)
 """
 function radius_from_anomaly(
     w::Float64,
@@ -68,11 +69,12 @@ function radius_from_anomaly(
 end
 
 """
-    radius_from_anomaly_derivative(u, a, e, model[, params])
+    radius_from_anomaly_derivative(w, a, e, model[, params])
 
 derivative of the mapping from anomaly to radius. Default is Henon anomaly.
 
-@IMPROVE: right now, Hénon anomaly is hard-coded.
+@IMPROVE: right now, Hénon anomaly is hard-coded. (except for AnalyticIsochrone and 
+SemiAnalyticPlummer, for which it is redefined using specific anomaly)
 """
 function radius_from_anomaly_derivative(
     w::Float64,

src/mappings/frequencies/derivatives.jl

@@ -11,7 +11,7 @@
 """
     αβ_from_ae_internal_derivatives(a, e, model, params)
 
-use the defined function Θ(u) to compute frequency ratios integrals
+use the defined function _Θ(w) to compute frequency ratios integrals
 AND DERIVATIVES
 """
 function αβ_from_ae_internal_derivatives(
@@ -20,7 +20,7 @@ function αβ_from_ae_internal_derivatives(
     model::Potential,
     params::OrbitalParameters=OrbitalParameters()
 )::Tuple{Float64,Float64,Float64,Float64,Float64,Float64}
-    tola, tole = params.TOLA, _eccentricity_tolerance(a, params.rc, params.TOLECC)
+    tola, tole = params.TOLA, _eccentricity_tolerance(a/params.rc, params.TOLECC)
     if (a < tola) || (e < tole) || (e > 1.0-tole)
         # For edge cases: Derivation outside the integral
         # Function to differentiate
@@ -35,11 +35,12 @@ function αβ_from_ae_internal_derivatives(
         return α, β, ∂α∂a, ∂β∂a, ∂α∂e, ∂β∂e
     end
     # Derivation inside the integral
-    # using Θ calculations to compute frequencies: leans heavily on Θ from Ufunc.jl
+    # using Θ calculations to compute frequencies: leans heavily on _Θ from 
+    # radial_velocity.jl
     # @IMPROVE: EDGE could be adaptive based on circularity and small-ness of rperi
 
-    # WARNING !! Strong assumption:
-    # r(u) = a(1+e * _henonf(w))
+    # @WARNING !! Strong assumption:
+    # r(w) = a(1+e * _henonf(w))
     function w6func(w::Float64)
         # push integration forward on eight different quantities:
         # 1. Θ                          → α
@@ -49,7 +50,7 @@ function αβ_from_ae_internal_derivatives(
         # 5. ∂Θ/∂a/r^2                  → ∂β∂a
         # 6. (∂Θ/∂e - 2af(w)Θ/r )/r^2   → ∂β∂e
 
-        integrand = Θ(w, a, e, model, params)
+        integrand = _Θ(w, a, e, model, params)
         ∂Θ∂a, ∂Θ∂e = _Θ_derivatives_ae(w, a, e, model, params)
 
         r = radius_from_anomaly(w, a, e, model, params)

src/mappings/frequencies/radial_velocity.jl

@@ -89,16 +89,16 @@ end
 #
 ########################################################################
 """
-    Θ(w, a, e, model[, params])
+    _Θ(w, a, e, model[, params])
 
-cured inverse radial velocity, Θ(u) = (dr/du)/v_rad, at anomaly `u` on orbit `(a,e)`.
+cured inverse radial velocity, _Θ(w) = (dr/dw)/v_rad, at anomaly `w` on orbit `(a,e)`.
 
 @IMPROVE: right now, Hénon anomaly is hard-coded.
 @IMPROVE: the expansion is trying both Taylor expansion and if fails, extrapolation with 
 particular care need at boundary switch. Fix it to use only one.
 @IMPROVE: find a better name
 """
-function Θ(
+function _Θ(
     w::Float64,
     a::Float64,
     e::Float64,
@@ -127,7 +127,7 @@ end
 """
     _Θedge(w, a, e, model, params)
 
-same as `Θ(...)` for w close to +1,-1 (curing the 0/0 limit at peri/apocentre)
+same as `_Θ(...)` for w close to +1,-1 (curing the 0/0 limit at peri/apocentre)
 
 @IMPROVE: right now, Hénon anomaly is hard-coded.
 @IMPROVE: the expansion is trying both Taylor expansion and, if fails, extrapolation with 
@@ -170,9 +170,9 @@ function _Θedge(
         w2 = wl * (1 - 2*params.EDGE)
         w3 = wl * (1 - 3*params.EDGE)
 
-        Θ1 = Θ(w1, a, e, model, params)
-        Θ2 = Θ(w2, a, e, model, params)
-        Θ3 = Θ(w3, a, e, model, params)
+        Θ1 = _Θ(w1, a, e, model, params)
+        Θ2 = _Θ(w2, a, e, model, params)
+        Θ3 = _Θ(w3, a, e, model, params)
 
         # It in fact is an extrapolation
         return _interpolate_order_2(w, w1, Θ1, w2, Θ2, w3, Θ3)
@@ -223,7 +223,7 @@ function angles_gradient(
     a::Float64,
     e::Float64,
     model::Potential,
-    params::OrbitalParameters;
+    params::OrbitalParameters=OrbitalParameters();
     L::Float64=0.0,
     Ω1::Float64=0.0,
     Ω2::Float64=0.0
@@ -239,22 +239,22 @@ function angles_gradient(
     # Current location of the radius, r=r(w)
     rval = radius_from_anomaly(w, a, e, model, params)
     # Current value of the radial frequency integrand (almost dθ/dw)
-    gval = Θ(w, a, e, model, params)
+    gval = _Θ(w, a, e, model, params)
     # Angles gradient (dθ1/dw, dθ2/dw)
     return Ω1*gval, (Ω2 - L/(rval^(2)))*gval
 end
 
 
 ########################################################################
 #
-# Θ(u) derivatives w.r.t. a and e
+# Θ(w) derivatives w.r.t. a and e
 #
 ########################################################################
 
 """
     _Θ_derivatives_ae(w, a, e, model, params)
 
-numerical differentiation of Θ w.r.t. semimajor axis and eccentricity
+numerical differentiation of `_Θ` w.r.t. semimajor axis and eccentricity
 
 @IMPROVE: right now the derivative procedure is hard-coded. Should incorparate it 
 in the OrbitalParameters structure (not only the steps)
@@ -270,7 +270,7 @@ function _Θ_derivatives_ae(
 )::Tuple{Float64,Float64}
 
     # Function to differentiate
-    fun(atmp::Float64, etmp::Float64) = Θ(w, atmp, etmp, model, params)
+    fun(atmp::Float64, etmp::Float64) = _Θ(w, atmp, etmp, model, params)
     # Perform differentiation
     _, ∂Θ∂a, ∂Θ∂e = _derivatives_ae(fun, a, e, params)