gradient handling signature changes

Author: michael-petersen

.github/workflows/ci.yml

@@ -19,7 +19,7 @@ jobs:
               Pkg.add("BenchmarkTools")
               Pkg.add("Test")
               Pkg.add(url="https://github.com/JuliaStellarDynamics/OrbitalElements.jl.git")
-              Pkg.add(url="https://github.com/JuliaStellarDynamics/DistributionFunctions.jl.git")
+              Pkg.add(url="https://github.com/JuliaStellarDynamics/DistributionFunctions.jl.git")'
             ls
             julia --project=df -e '
               using Pkg

Project.toml

@@ -1,6 +1,6 @@
 name = "DistributionFunctions"
 uuid = "c869f47d-2815-40f9-b874-25ebe83f43af"
-version = "0.0.6dev"
+version = "0.0.7dev"
 
 [deps]
 OrbitalElements = "a3b07092-bde3-4843-b84f-c597d614ec7b"

src/Analytic/Isochrone/isochrone.jl

@@ -2,12 +2,12 @@
 # Isochrone distribution functions (analytic)
 #####################################
 abstract type IsochroneDistributionFunction <: DistributionFunction end
-struct IsotropicIsochrone <: IsochroneDistributionFunction
-    potential::IsochronePotential # Potential model
+struct IsotropicIsochrone{modelT<:IsochronePotential} <: IsochroneDistributionFunction
+    potential::modelT # Potential model
 end
-struct OsipkovMerrittIsochrone <: IsochroneDistributionFunction
+struct OsipkovMerrittIsochrone{modelT<:IsochronePotential} <: IsochroneDistributionFunction
     ra::Float64                # Anisotropy radius
-    potential::IsochronePotential # Potential model
+    potential::modelT # Potential model
 end
 
 """

src/Analytic/Isochrone/isotropic.jl

@@ -10,10 +10,10 @@ end
 
 
 
-"""isochroneisotropicDF(E[,bc,M,G])
+"""Distribution(E[,bc,M,G])
 the isotropic DF
 """
-function F(EL::Tuple{Float64,Float64}, df::IsotropicIsochrone)
+function Distribution(EL::Tuple{Float64,Float64}, df::IsotropicIsochrone)
     E,L = EL
 
     scaleEnergy = - df.potential.G * df.potential.M / df.potential.bc
@@ -29,10 +29,10 @@ function F(EL::Tuple{Float64,Float64}, df::IsotropicIsochrone)
 end
 
 
-"""isochroneisotropicdDFdE(E[,bc,M,G])
+"""DFDE(EL::Tuple{Float64,Float64}, df::IsotropicIsochrone)
 the isotropic DF derivative w.r.t. E
 """
-function dFdE(EL::Tuple{Float64,Float64}, df::IsotropicIsochrone)
+function DFDE(EL::Tuple{Float64,Float64}, df::IsotropicIsochrone)
 
     E,L = EL
 
@@ -50,14 +50,12 @@ function dFdE(EL::Tuple{Float64,Float64}, df::IsotropicIsochrone)
 end
 
 
-function ndFdJ(EL::Tuple{Float64,Float64},ΩΩ::Tuple{Float64,Float64},resonance::Resonance, df::IsotropicIsochrone)
-    Ω1,Ω2 = ΩΩ
-    ndotΩ = resonance.number[1]*Ω1 + resonance.number[2]*Ω2
-
-    dFdEval = dFdE(EL,df)
+"""DFDL(EL::Tuple{Float64,Float64}, df::IsotropicIsochrone)
+the isotropic DF derivative w.r.t. L: always zero!
+"""
+function DFDL(EL::Tuple{Float64,Float64}, df::IsotropicIsochrone)::Float64
 
-    # Value of n.dF/dJ
-    result = ndotΩ*dFdEval
+    return 0.0
 
-    return result
 end
+

src/Analytic/Isochrone/osipkovmerritt.jl

@@ -51,7 +51,7 @@ end
 Saha distribution function
 ra is the anisotropy radius
 """
-function F(EL::Tuple{Float64,Float64}, df::OsipkovMerrittIsochrone)
+function Distribution(EL::Tuple{Float64,Float64}, df::OsipkovMerrittIsochrone)::Float64
 
     Q       = osipkovmerritt_Q(EL, df)
     scaleDF = dfscale(df)
@@ -93,15 +93,27 @@ function osipkovmerritt_dFdQ(EL::Tuple{Float64,Float64}, df::OsipkovMerrittIsoch
 
 end
 
+function DFDE(EL::Tuple{Float64,Float64}, df::OsipkovMerrittIsochrone)::Float64
 
+    Q = osipkovmerritt_Q(EL,df)
 
-function ndFdJ(EL::Tuple{Float64,Float64},ΩΩ::Tuple{Float64,Float64},resonance::Resonance, df::OsipkovMerrittIsochrone)
+    # If Q is outside of the [0,1]--range, we set the function to 0.0
+    # ATTEN TION, this is a lazy implementation -- it would have been much better to restrict the integration domain
+    if (!(0.0 <= Q <= 1.0)) # If Q is outside of the [0,1]-range, we set the function to 0
+        return 0.0 # Outside of the physically allowed orbital domain
+    end
 
+    # Value of dF/dQ
+    dFdQ = osipkovmerritt_dFdQ(EL,df)
+
+    # Value of dQ/dE
+    dQdE = osipkovmerritt_dQdE(EL,df)
+
+    return dFdQ * dQdE
+end
+
+function DFDL(EL::Tuple{Float64,Float64}, df::OsipkovMerrittIsochrone)::Float64
 
-    Ω1,Ω2 = ΩΩ
-    n1,n2 = resonance.number[1],resonance.number[2]
-    ndotΩ = n1*Ω1 + n2*Ω2
-    
     Q = osipkovmerritt_Q(EL,df)
 
     # If Q is outside of the [0,1]--range, we set the function to 0.0
@@ -113,12 +125,8 @@ function ndFdJ(EL::Tuple{Float64,Float64},ΩΩ::Tuple{Float64,Float64},resonance
     # Value of dF/dQ
     dFdQ = osipkovmerritt_dFdQ(EL,df)
 
-    # Values of dQ/dE, dQ/dL
-    dQdE, dQdL = osipkovmerritt_dQdE(EL,df), osipkovmerritt_dQdL(EL,df)
-
-    # Value of n.dF/dJ
-    result = dFdQ*(dQdE*ndotΩ + n2*dQdL)
-
-    return result
+    # Value dQ/dL
+    dQdL = osipkovmerritt_dQdL(EL,df)
 
+    return dFdQ * dQdL
 end

src/Analytic/Plummer/isotropic.jl

@@ -13,7 +13,7 @@ end
 """
 The distribution function for the isotropic Plummer model
 """
-function F(EL::Tuple{Float64,Float64}, df::IsotropicPlummer)
+function Distribution(EL::Tuple{Float64,Float64}, df::IsotropicPlummer)
     E,L = EL
 
     scaleEnergy = - df.potential.G * df.potential.M / df.potential.bc
@@ -31,7 +31,9 @@ the derivative of the isotropic distribution function
 F[EN_] = (24 Sqrt[2])/(7 Pi^3) * EN^(7/2)
 D[F[EN], {EN}]
 """
-function dFdE(E::Float64,df::IsotropicPlummer)
+function DFDE(EL::Tuple{Float64,Float64},df::IsotropicPlummer)
+
+    E,L = EL
 
     scaleEnergy = - df.potential.G * df.potential.M / df.potential.bc
     scaleDF     = dfscale(df)
@@ -42,19 +44,9 @@ function dFdE(E::Float64,df::IsotropicPlummer)
 
 end
 
+function DFDL(EL::Tuple{Float64,Float64},df::IsotropicPlummer)::Float64
 
-"""
-calculate n.dF/dE, for the isotropic plummer case.
-"""
-function ndFdJ(EL::Tuple{Float64,Float64},ΩΩ::Tuple{Float64,Float64},resonance::Resonance, df::IsotropicPlummer)
-    E,L = EL
-    Ω1,Ω2 = ΩΩ
-    ndotΩ = resonance.number[1]*Ω1 + resonance.number[2]*Ω2
-
-    dFdE = dFdE(E,df)
+    return 0.0
 
-    # Value of n.dF/dJ
-    result = ndotΩ*dFdE
+end
 
-    return result
-end

src/Analytic/Plummer/osipkovmerrit.jl

@@ -48,7 +48,7 @@ end
 """
 the anisotropic distribution function from PlummerPlus
 """
-function F(EL::Tuple{Float64,Float64}, df::OsipkovMerrittPlummer)
+function Distribution(EL::Tuple{Float64,Float64}, df::OsipkovMerrittPlummer)
 
     Q       = osipkovmerritt_Q(EL, df)
     scaleDF = dfscale(df)
@@ -62,34 +62,26 @@ function F(EL::Tuple{Float64,Float64}, df::OsipkovMerrittPlummer)
 
 end
 
+function DFDE(EL::Tuple{Float64,Float64}, df::OsipkovMerrittPlummer)::Float64
 
-"""
-the derivative of the anisotropic distribution function
-"""
-function osipkovmerritt_dFdQ(EL::Tuple{Float64,Float64}, df::OsipkovMerrittPlummer)
+    Q = osipkovmerritt_Q(EL, df)
 
-    Q       = osipkovmerritt_Q(EL, df)
-    scaleDF = dfscale(df)
+    # If Q is outside of the [0,1]--range, we set the function to 0.0
+    # ATTENTION, this is a lazy implementation -- it would have been much better to restrict the integration domain
+    if (!(0.0 <= Q <= 1.0)) # If Q is outside of the [0,1]-range, we set the function to 0
+        return 0.0 # Outside of the physically allowed orbital domain
+    end
 
-    # dimensionless anisotropy radius: ALREADY SQUARED
-    gamma = (df.potential.bc/df.ra)^2
+    # Value of dF/dQ
+    dFdQ = osipkovmerritt_dFdQ(EL,df)
 
-    prefactor = (df.potential.M/scaleDF) * ((3*sqrt(2))/(4*(pi^3)))
-    return prefactor * sqrt(Q) * ( (gamma)*(3 - 16(Q^2)) + 16(Q^2)) 
+    # Value of dQ/dE, dQ/dL
+    dQdE = osipkovmerritt_dQdE(EL,df)
 
+    return dFdQ*dQdE
 end
 
-
-"""
-calculate n.dF/dJ, for the osipkov-merritt plummer case.
-
-@IMPROVE, make a distribution function fix. perhaps only integrate in regions where the DF is positive?
-"""
-function ndFdJ(EL::Tuple{Float64,Float64},ΩΩ::Tuple{Float64,Float64},resonance::Resonance, df::OsipkovMerrittPlummer)
-
-    Ω1,Ω2 = ΩΩ
-    n1,n2 = resonance.number[1],resonance.number[2]
-    ndotΩ = n1*Ω1 + n2*Ω2
+function DFDL(EL::Tuple{Float64,Float64}, df::OsipkovMerrittPlummer)::Float64
 
     Q = osipkovmerritt_Q(EL, df)
 
@@ -102,12 +94,26 @@ function ndFdJ(EL::Tuple{Float64,Float64},ΩΩ::Tuple{Float64,Float64},resonance
     # Value of dF/dQ
     dFdQ = osipkovmerritt_dFdQ(EL,df)
 
-    # Values of dQ/dE, dQ/dL
-    dQdE, dQdL = osipkovmerritt_dQdE(EL,df), osipkovmerritt_dQdL(EL,df)
+    # Value of  dQ/dL
+    dQdL =  osipkovmerritt_dQdL(EL,df)
+
+    return dFdQ*dQdL
+end
 
-    # Value of n.dF/dJ
-    result = dFdQ*(dQdE*ndotΩ + n2*dQdL)
 
-    return result
+"""
+the derivative of the anisotropic distribution function
+"""
+function osipkovmerritt_dFdQ(EL::Tuple{Float64,Float64}, df::OsipkovMerrittPlummer)
+
+    Q       = osipkovmerritt_Q(EL, df)
+    scaleDF = dfscale(df)
+
+    # dimensionless anisotropy radius: ALREADY SQUARED
+    gamma = (df.potential.bc/df.ra)^2
+
+    prefactor = (df.potential.M/scaleDF) * ((3*sqrt(2))/(4*(pi^3)))
+    return prefactor * sqrt(Q) * ( (gamma)*(3 - 16(Q^2)) + 16(Q^2)) 
 
 end
+

src/Analytic/RazorThinDiscs/mestel.jl

@@ -5,54 +5,32 @@ MestelDisc([potential])
 
 Mestel disc distribution function.
 """
-function MestelDisc(;potential::MestelPotential=MestelPotential(),q::IntorFloat=11.44,ξDF::Float64=1.0,G::Float64=1.0)
-    return MestelDisc(potential,q,ξDF,G)
+function MestelDisc(;potential::MestelPotential=MestelPotential(),q::IntorFloat=11.44,G::Float64=1.0)
+    return MestelDisc(potential,q,G)
 end
 
 """
-    MestelDF(EL::Tuple{Float64,Float64},df::MestelDisc)
+    MestelDistribution(EL::Tuple{Float64,Float64},df::MestelDisc)
 Mestel distribution function.
 """
-function F(EL::Tuple{Float64,Float64},df::MestelDisc)::Float64
+function Distribution(EL::Tuple{Float64,Float64},df::MestelDisc)::Float64
 
-    return MestelDF(EL,df)
+    return MestelDistribution(EL,df)
 end
 
 """
     dFdE(EL::Tuple{Float64,Float64},df::MestelDisc)
 Mestel DF derivative w.r.t. E.
 """
-function dFdE(EL::Tuple{Float64,Float64},df::MestelDisc)::Float64
+function DFDE(EL::Tuple{Float64,Float64},df::MestelDisc)::Float64
     return MesteldFdE(EL,df)
 end
 
 """
     dFdL(E, L[, C, q, sigma])
 Mestel DF derivative w.r.t. E.
 """
-function dFdL(EL::Tuple{Float64,Float64},df::MestelDisc)::Float64
+function DFDL(EL::Tuple{Float64,Float64},df::MestelDisc)::Float64
     return MesteldFdL(EL,df)
 end
 
-
-"""
-    ndFdJ(EL::Tuple{Float64,Float64},ΩΩ::Tuple{Float64,Float64},resonance::Resonance, df::MestelDisc)
-Zang star DF derivative w.r.t. the actions J.
-"""
-function ndFdJ(EL::Tuple{Float64,Float64},ΩΩ::Tuple{Float64,Float64},resonance::Resonance, df::MestelDisc)::Float64
-
-    E,L = EL
-    Ω1,Ω2 = ΩΩ
-    n1,n2 = resonance.number
-    ndotΩ = n1*Ω1 + n2*Ω2
-
-    if L <= 0.
-        println("WARNING: L <= 0.")
-        return 0.
-    end
-
-    dDFdE = dFdE(EL,df)
-    dDFdL = dFdL(EL,df)
-    
-    return df.ξDF * (ndotΩ*dDFdE + n2*dDFdL)
-end

src/Analytic/RazorThinDiscs/miyamoto.jl

@@ -10,7 +10,7 @@ using HypergeometricFunctions
 
 Miyamoto distribution function for Kuzmin-Toomre disc.
 """
-function miyamoto_distribution(E::Float64, L::Float64; mM::Int64=1)
+function Distribution(E::Float64, L::Float64; mM::Int64=1)
     return (
         (2mM + 3)
         * (-E)^(2mM + 2)
@@ -24,7 +24,7 @@ end
 
 Miyamoto distribution derivative w.r.t. energy.
 """
-function miyamoto_dFdE(E::Float64, L::Float64; mM::Int64=1)
+function DFDE(E::Float64, L::Float64; mM::Int64=1)
     return (
         E^(2mM)
         * (1 + mM)
@@ -42,7 +42,7 @@ end
 
 Miyamoto distribution derivative w.r.t. angular momentum.
 """
-function miyamoto_dFdL(E::Float64, L::Float64; mM::Int64=1)
+function DFDL(E::Float64, L::Float64; mM::Int64=1)
 
     return -(
         2
@@ -56,22 +56,3 @@ function miyamoto_dFdL(E::Float64, L::Float64; mM::Int64=1)
     )
 end
 
-"""
-    miyamoto_ndFdJ(n1, n2, E, L, ndotOmega[, mM])
-
-Miyamoto distribution derivatives w.r.t. the actions.
-"""
-function miyamoto_ndFdJ(
-    n1::Int64,
-    n2::Int64,
-    E::Float64,
-    L::Float64,
-    ndotOmega::Float64;
-    mM::Int64=1
-)
-
-    dFdE = miyamoto_dFdE(E, L; mM=mM)
-    dFdL = miyamoto_dFdL(E, L; mM=mM)
-    
-    return ndotOmega*dFdE + n2*dFdL
-end