add compatibility with DistributionFunctions.jl; some renaming

Project.toml

@@ -13,6 +13,7 @@ DistributionFunctions = "c869f47d-2815-40f9-b874-25ebe83f43af"
 
 [compat]
 julia = "1.6.7"
+OrbitalElements = "2"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

src/GFunc.jl

@@ -5,6 +5,29 @@
 #
 ########################################################################
 
+
+"""
+    ndFdJ(EL::Tuple{Float64,Float64},ΩΩ::Tuple{Float64,Float64},res::Resonance, distributionfunction::DistributionFunction)
+
+Distribution function derivative for `distributionfunction` for a given `E`,`L`.
+
+    TODO: make frequencies an optional input
+"""
+function ndFdJ(EL::Tuple{Float64,Float64},ΩΩ::Tuple{Float64,Float64},resonance::Resonance, df::EnergyAngularMomentumDistributionFunction)
+
+        DFDEval = DFDE(EL,df)
+        DFDLval = DFDL(EL,df)
+        Ω1,Ω2   = ΩΩ
+        n1,n2   = resonance.number[1],resonance.number[2]
+        ndotΩ   = n1*Ω1 + n2*Ω2
+
+        return DFDEval * ndotΩ + DFDLval * n2
+
+    end
+
+
+
+
 """
     MakeGu(ndFdJ,n1,n2,Wdata,tabu[,params])
 

src/ModeComputation/Determinant.jl

@@ -10,18 +10,38 @@ VERBOSE flag rules
 
 
 """
-    detXi(IMat,tabM[,ξ])
+    permitivity_determinant(identity_matrix, tabM; ξ=1.0) -> ComplexF64
 
-determinant of the susceptibility matrix I - ξ*M(ω) 
-for known M(ω) and active fraction ξ (default 1).
+Calculate the determinant of the permitivity matrix `I - ξ*M(ω)`.
+
+# Arguments
+- `identity_matrix::AbstractMatrix{ComplexF64}`: The identity matrix.
+- `tabM::AbstractMatrix{ComplexF64}`: The matrix `M(ω)` representing the frequency-dependent response.
+- `ξ::Float64=1.0`: The active fraction, defaulting to 1.
+
+# Returns
+- `ComplexF64`: The determinant of the susceptibility matrix `I - ξ*M(ω)`.
+
+# Details
+This function computes the determinant of the matrix `I - ξ*M(ω)`, assuming the resulting matrix is symmetric. The function uses Julia's `Symmetric` type to inform the determinant computation that the matrix is symmetric, which can optimize the computation. The result is the real portion of the determinant.
+
+# Examples
+```julia
+identity_matrix = [1 + 0im 0 + 0im; 0 + 0im 1 + 0im]
+tabM = [0.5 + 0.1im 0.2 + 0.2im; 0.2 + 0.2im 0.3 + 0.1im]
+ξ = 0.9
+
+det_value = permitivity_determinant(identity_matrix, tabM, ξ=ξ)
+println(det_value)
+# Output: ComplexF64 value representing the determinant of I - ξ*M(ω)
 """
-function detXi(IMat::AbstractMatrix{ComplexF64},
+function permitivity_determinant(identity_matrix::AbstractMatrix{ComplexF64},
                tabM::AbstractMatrix{ComplexF64};
                ξ::Float64=1.0)::ComplexF64
 
     # Computing the determinant of (I-M).
     # ATTENTION, we tell julia that the matrix is symmetric
-    val = det(Symmetric(IMat-ξ*tabM))
+    val = det(Symmetric(identity_matrix-ξ*tabM))
 
     # only save the real portion
     return val # Output
@@ -38,12 +58,12 @@ function RunDeterminant(ωlist::Array{ComplexF64},
 
     # Preparinng computations of the response matrices
     tabMlist, tabaMcoef, tabωminωmax, FHTlist = PrepareM(Threads.nthreads(),FHT,params)
-    IMatlist = makeIMat(params.nradial,Threads.nthreads())
+    identity_matrixlist = make_identity_matrix(params.nradial,Threads.nthreads())
 
     # how many omega values are we computing?
     nω = length(ωlist)
     # allocate containers for determinant and min eigenvalue
-    tabdetXi = zeros(ComplexF64,nω)
+    tabpermitivity_determinant = zeros(ComplexF64,nω)
 
     (params.VERBOSE > 0) && println("LinearResponse.Determinant.RunDeterminant: computing $nω frequency values.")
 
@@ -58,7 +78,7 @@ function RunDeterminant(ωlist::Array{ComplexF64},
             tabM!(ωlist[i],tabMlist[k],tabaMcoef,tabωminωmax,FHTlist[k],params)
         end
 
-        tabdetXi[i] = detXi(IMatlist[k],tabMlist[k],ξ=ξ)
+        tabpermitivity_determinant[i] = permitivity_determinant(identity_matrixlist[k],tabMlist[k],ξ=ξ)
     end
 
     h5open(DetFilename(params,ξ=ξ), "w") do file
@@ -68,11 +88,11 @@ function RunDeterminant(ωlist::Array{ComplexF64},
         write(file,"omega",real(ωlist))
         write(file,"eta",imag(ωlist))
         # Results
-        write(file,"det",tabdetXi)
+        write(file,"det",tabpermitivity_determinant)
         # Parameters
         WriteParameters(file,params)
     end
-    return tabdetXi
+    return tabpermitivity_determinant
 end
 
 
@@ -91,7 +111,7 @@ function FindZeroCrossing(Ωguess::Float64,ηguess::Float64,
 
     # Preparinng computations of the response matrices
     MMat, tabaMcoef, tabωminωmax = PrepareM(params)
-    IMat = makeIMat(params.nradial)
+    identity_matrix = make_identity_matrix(params.nradial)
 
     omgval = Ωguess + im*ηguess
     domega = 1.e-4
@@ -109,11 +129,11 @@ function FindZeroCrossing(Ωguess::Float64,ηguess::Float64,
 
         tabM!(omgval,MMat,tabaMcoef,tabωminωmax,FHT,params)
 
-        centralvalue = detXi(IMat,MMat,ξ=ξ)
+        centralvalue = permitivity_determinant(identity_matrix,MMat,ξ=ξ)
 
         tabM!(omgvaloff,MMat,tabaMcoef,tabωminωmax,FHT,params)
 
-        offsetvalue = detXi(IMat,MMat,ξ=ξ)
+        offsetvalue = permitivity_determinant(identity_matrix,MMat,ξ=ξ)
 
         # ignore the imaginary part
         derivative = real(offsetvalue - centralvalue)/domega

src/ModeComputation/FindPole.jl

@@ -2,7 +2,7 @@
 
 
 function GoStep(omgval::ComplexF64,
-                IMat::Array{ComplexF64,2},MMat::Array{ComplexF64,2},
+                identity_matrix::Array{ComplexF64,2},MMat::Array{ComplexF64,2},
                 FHT::FiniteHilbertTransform.AbstractFHT,
                 tabaMcoef::Array{Float64,4},
                 tabωminωmax::Matrix{Float64},
@@ -14,43 +14,43 @@ function GoStep(omgval::ComplexF64,
     tabM!(omgval,MMat,tabaMcoef,tabωminωmax,FHT,params)
 
     # compute the determinant of I-M
-    detXival = detXi(IMat,MMat,ξ=ξ)
+    permitivity_determinant_val = permitivity_determinant(identity_matrix,MMat,ξ=ξ)
     # if this value is bad, should immediately break
 
     # now compute a Jacobian via finite differences
     riomgval = omgval + DSIZE
     upomgval = omgval + DSIZE*1im
 
     tabM!(riomgval,MMat,tabaMcoef,tabωminωmax,FHT,params)
-    detXivalri = detXi(IMat,MMat,ξ=ξ)
+    permitivity_determinant_valri = permitivity_determinant(identity_matrix,MMat,ξ=ξ)
 
     tabM!(upomgval,MMat,tabaMcoef,tabωminωmax,FHT,params)
-    detXivalup = detXi(IMat,MMat,ξ=ξ)
+    permitivity_determinant_valup = permitivity_determinant(identity_matrix,MMat,ξ=ξ)
 
-    dXirir = real(detXivalri-detXival)/DSIZE
-    dXirii = imag(detXivalri-detXival)/DSIZE
-    dXiupr = real(detXivalup-detXival)/DSIZE
-    dXiupi = imag(detXivalup-detXival)/DSIZE
+    dXirir = real(permitivity_determinant_valri-permitivity_determinant_val)/DSIZE
+    dXirii = imag(permitivity_determinant_valri-permitivity_determinant_val)/DSIZE
+    dXiupr = real(permitivity_determinant_valup-permitivity_determinant_val)/DSIZE
+    dXiupi = imag(permitivity_determinant_valup-permitivity_determinant_val)/DSIZE
 
     # jacobian = [dXirir dXiupr ; dXirii dXiupi]
-    # curpos = [real(detXival);imag(detXival)]
+    # curpos = [real(permitivity_determinant_val);imag(permitivity_determinant_val)]
     # increment = jacobian \ (-curpos)
 
-    increment1, increment2 = OrbitalElements._inverse_2d(dXirir,dXiupr,dXirii,dXiupi,-real(detXival),-imag(detXival))
+    increment1, increment2 = OrbitalElements._inverse_2d(dXirir,dXiupr,dXirii,dXiupi,-real(permitivity_determinant_val),-imag(permitivity_determinant_val))
 
     # omgval += increment[1] + increment[2]*1im
     omgval += increment1 + increment2*1im
     if params.VERBOSE > 0
-        println("detXi=$detXival")
+        println("permitivity_determinant=$permitivity_determinant_val")
         println("omega=$omgval")
     end
 
-    return omgval,detXival
+    return omgval,permitivity_determinant_val
 end
 
 
 function FindDeterminantZero(startingomg::ComplexF64,
-                             IMat::Array{ComplexF64,2},MMat::Array{ComplexF64,2},
+                             identity_matrix::Array{ComplexF64,2},MMat::Array{ComplexF64,2},
                              FHT::FiniteHilbertTransform.AbstractFHT,
                              tabaMcoef::Array{Float64,4},
                              tabωminωmax::Matrix{Float64},
@@ -61,15 +61,15 @@ function FindDeterminantZero(startingomg::ComplexF64,
                              NSEARCH::Int64=100)
 
     # initial values
-    detXival = 1.0
+    permitivity_determinant_val = 1.0
     omgval = startingomg
     stepnum = 0
 
-    while abs(detXival) > TOL
+    while abs(permitivity_determinant_val) > TOL
 
-        omgval,detXival = GoStep(omgval,IMat,MMat,FHT,tabaMcoef,tabωminωmax,params,ξ=ξ,DSIZE=DSIZE)
+        omgval,permitivity_determinant_val = GoStep(omgval,identity_matrix,MMat,FHT,tabaMcoef,tabωminωmax,params,ξ=ξ,DSIZE=DSIZE)
 
-        if abs(detXival) > 1.0
+        if abs(permitivity_determinant_val) > 1.0
             break
         end
 
@@ -82,8 +82,8 @@ function FindDeterminantZero(startingomg::ComplexF64,
     end
 
     # what happens when this goes wrong?
-    if abs(detXival) < TOL
-        return omgval,detXival
+    if abs(permitivity_determinant_val) < TOL
+        return omgval,permitivity_determinant_val
     end
 
     return NaN,NaN
@@ -99,9 +99,9 @@ function FindPole(startingomg::ComplexF64,
 
     # Preparinng computations of the response matrices
     MMat, tabaMcoef, tabωminωmax = PrepareM(params)
-    IMat = makeIMat(params.nradial)
+    identity_matrix = make_identity_matrix(params.nradial)
 
-    bestomg,detval = FindDeterminantZero(startingomg,IMat,MMat,FHT,tabaMcoef,tabωminωmax,params,ξ=ξ,TOL=TOL,DSIZE=DSIZE)
+    bestomg,detval = FindDeterminantZero(startingomg,identity_matrix,MMat,FHT,tabaMcoef,tabωminωmax,params,ξ=ξ,TOL=TOL,DSIZE=DSIZE)
 
     (params.VERBOSE >= 0) && println("Best O for n1max=$(params.n1max),nradial=$(params.nradial) == $bestomg")
 

src/ModeComputation/Matrices.jl

@@ -9,53 +9,105 @@ VERBOSE flag rules
 """
 
 """
-@TO DESCRIBE
+    RunMatrices(ωlist::Array{ComplexF64}, FHT::FiniteHilbertTransform.AbstractFHT, params::LinearParameters)
+
+Compute the response matrices for a given list of frequencies and parameters.
+
+# Arguments
+- `ωlist::Array{ComplexF64}`: An array of complex frequencies.
+- `FHT::FiniteHilbertTransform.AbstractFHT`: An instance of a finite Hilbert transform.
+- `params::LinearParameters`: A structure containing the linear parameters required for the computation.
+
+# Returns
+- `tabRMreal::Array{Float64, 3}`: A 3D array containing the real parts of the response matrices for each frequency.
+- `tabRMimag::Array{Float64, 3}`: A 3D array containing the imaginary parts of the response matrices for each frequency.
+
+# Description
+This function prepares the computation of response matrices for a given list of complex frequencies (`ωlist`). It initializes necessary data structures, iterates over each frequency to compute the response matrix, and stores the real and imaginary parts of these matrices. The results are then saved to an HDF5 file.
+
+# Details
+- The function uses multithreading to parallelize the computation across available threads.
+- The computed matrices are saved in an HDF5 file, with both real and imaginary parts stored separately.
+- The function provides optional timing information for the second frequency computation if verbosity (`params.VERBOSE`) is enabled.
+
+# Examples
+```julia
+# Example usage
+ωlist = [1.0 + 0.1im, 2.0 + 0.2im, 3.0 + 0.3im]
+FHT = FiniteHilbertTransform()
+params = LinearParameters(nradial=3, VERBOSE=1)
+tabRMreal, tabRMimag = RunMatrices(ωlist, FHT, params)
+```
 """
 function RunMatrices(ωlist::Array{ComplexF64},
-                     FHT::FiniteHilbertTransform.AbstractFHT,
-                     params::LinearParameters)
+    FHT::FiniteHilbertTransform.AbstractFHT,
+    params::LinearParameters)
+
+    # Validate inputs
+    if isempty(ωlist)
+        throw(ArgumentError("ωlist cannot be empty"))
+    end
+    if params.nradial <= 0
+        throw(ArgumentError("params.nradial must be positive"))
+    end
 
-    # Preparinng computations of the response matrices
-    tabMlist, tabaMcoef, tabωminωmax, FHTlist = PrepareM(Threads.nthreads(),FHT,params)
+    # Preparing computations of the response matrices
+    try
+        tabMlist, tabaMcoef, tabωminωmax, FHTlist = PrepareM(Threads.nthreads(), FHT, params)
+    catch e
+        println("Error preparing matrices: ", e)
+        return nothing
+    end
 
-    # how many omega values are we computing?
+    # Number of frequency values
     nω = length(ωlist)
-    # allocate containers for the matrices
+    # Allocate containers for the matrices
     nradial = params.nradial
-    tabRMreal = zeros(Float64,nradial,nradial,nω)
-    tabRMimag = zeros(Float64,nradial,nradial,nω)
+    tabRMreal = zeros(Float64, nradial, nradial, nω)
+    tabRMimag = zeros(Float64, nradial, nradial, nω)
 
     (params.VERBOSE > 0) && println("LinearResponse.Xi.RunMatrices: computing $nω frequency values.")
 
-    # loop through all frequencies
+    # Loop through all frequencies using multithreading with dynamic scheduling
     Threads.@threads for i = 1:nω
-
         k = Threads.threadid()
 
-        if (i==2) && (params.VERBOSE>0) # skip the first in case there is compile time built in
-            @time tabM!(ωlist[i],tabMlist[k],tabaMcoef,tabωminωmax,FHTlist[k],params)
-        else
-            tabM!(ωlist[i],tabMlist[k],tabaMcoef,tabωminωmax,FHTlist[k],params)
-        end
+        try
+            # Time the second frequency computation if verbosity is enabled
+            if (i == 2) && (params.VERBOSE > 0)
+                @time tabM!(ωlist[i], tabMlist[k], tabaMcoef, tabωminωmax, FHTlist[k], params)
+            else
+                tabM!(ωlist[i], tabMlist[k], tabaMcoef, tabωminωmax, FHTlist[k], params)
+            end
 
-        for q = 1:nradial
-            for p = 1:nradial
-                tabRMreal[p,q,i] = real(tabMlist[k][p,q])
-                tabRMimag[p,q,i] = imag(tabMlist[k][p,q])
+            # Store real and imaginary parts of the matrix
+            for q = 1:nradial
+                for p = 1:nradial
+                    tabRMreal[p, q, i] = real(tabMlist[k][p, q])
+                    tabRMimag[p, q, i] = imag(tabMlist[k][p, q])
+                end
             end
+        catch e
+            # what kinds of diagnostics would be helpful here for failure modes?
+            (params.VERBOSE > 0) && println("Error in thread $k while processing frequency $i: ", e)
         end
     end
 
-    h5open(MatFilename(params), "w") do file
-        # ω grid
-        write(file,"omega",real(ωlist))
-        write(file,"eta",imag(ωlist))
-        # Matrices
-        write(file,"MatricesR",tabRMreal)
-        write(file,"MatricesI",tabRMimag)
-        # Parameters
-        WriteParameters(file,params)
+    # Save results to an HDF5 file
+    try
+        h5open(MatFilename(params), "w") do file
+            # ω grid
+            write(file, "omega", real(ωlist))
+            write(file, "eta", imag(ωlist))
+            # Matrices
+            write(file, "MatricesR", tabRMreal)
+            write(file, "MatricesI", tabRMimag)
+            # Parameters
+            WriteParameters(file, params)
+        end
+    catch e
+        println("Error writing to HDF5 file: ", e)
     end
 
     return tabRMreal, tabRMimag
-end
+end