Merge pull request #16 from JuliaStellarDynamics/v0.10.1

v1.0.0dev-0
JuliaStellarDynamics · Mar 3, 2024 · c9f8f05 · c9f8f05
2 parents 99e221d + 312be33
commit c9f8f05
Show file tree

Hide file tree

Showing 14 changed files with 303 additions and 181 deletions.
diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
@@ -18,7 +18,7 @@ jobs:
         arch:
           - x64
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v1
         with:
           version: ${{ matrix.version }}
@@ -36,6 +36,6 @@ jobs:
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v1
+      - uses: codecov/codecov-action@v4
         with:
           file: lcov.info
diff --git a/Project.toml b/Project.toml
@@ -1,16 +1,11 @@
 name = "FiniteHilbertTransform"
 uuid = "68cdf181-988b-4042-8999-0b9541c1c3c2"
 authors = ["Mike Petersen <petersen@iap.fr>", "Mathieu Roule <roule@iap.fr>"]
-version = "0.10.0"
+version = "1.0.0dev-0"
 
 [deps]
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
-HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
-IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a"
-JLLWrappers = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
-Quadmath = "be4d8f0f-7fa4-5f49-b795-2f01399ab2dd"
 
 [compat]
 julia = "1.6.7"

diff --git a/README.md b/README.md
@@ -2,8 +2,11 @@
 # FiniteHilbertTransform.jl
 
 [![Test Builds](https://github.com/JuliaStellarDynamics/FiniteHilbertTransform.jl/actions/workflows/CI.yml/badge.svg)](https://github.com/JuliaStellarDynamics/FiniteHilbertTransform.jl/actions/workflows/CI.yml)
-[![codecov](https://codecov.io/github/JuliaStellarDynamics/FiniteHilbertTransform.jl/graph/badge.svg?token=LO51PVUTU1)](https://codecov.io/github/JuliaStellarDynamics/FiniteHilbertTransform.jl)
+[![codecov](https://codecov.io/github/JuliaStellarDynamics/FiniteHilbertTransform.jl/coverage.svg)](https://codecov.io/github/JuliaStellarDynamics/FiniteHilbertTransform.jl)
 [![image](https://github.com/JuliaStellarDynamics/FiniteHilbertTransform.jl/actions/workflows/documentation.yml/badge.svg?branch=documentation)](https://juliastellardynamics.github.io/FiniteHilbertTransform.jl/)
+[![image](https://img.shields.io/badge/julia-stable-blue)](https://github.com/JuliaStellarDynamics/OrbitalElements.jl/actions/workflows/devCI.yml)
+[![image](http://img.shields.io/badge/license-MIT-brightgreen.svg)](https://github.com/JuliaStellarDynamics/OrbitalElements.jl/blob/v2.0/LICENSE)
+[![image](http://img.shields.io/badge/DOI-10.48550/arXiv.2311.10630-blue.svg)](http://dx.doi.org/10.48550/arXiv.2311.10630)
 
 **FiniteHilbertTransform.jl** is a Julia package designed to compute the finite version of the Hilbert transformations. This toolbox is inspired by Tricomi's work on the finite Hilbert transform from 1957. In the context of gravitational dynamics, the finite Hilbert transform may be used as a scheme for analytic continuation to the lower half of the complex plane. See [Fouvry & Prunet (2022)](https://ui.adsabs.harvard.edu/abs/2022MNRAS.509.2443F/abstract), or [Petersen et al. (2024)](https://ui.adsabs.harvard.edu/abs/2023arXiv231110630P/abstract) for details.
 
@@ -16,9 +19,10 @@ To invoke Julia in the Terminal, you need to make sure that the `julia` command-
 See [here](https://julialang.org/downloads/platform/#optional_add_julia_to_path) for detailed instructions.
 
 Once Julia installed, obtain the `FiniteHilbertTransform.jl` library[^1][^2] and compile it by running:
-    ```
-    julia -e 'using Pkg; Pkg.add(url="https://github.com/JuliaStellarDynamics/FiniteHilbertTransform.jl.git")'
-    ```
+
+```
+julia -e 'using Pkg; Pkg.add(url="https://github.com/JuliaStellarDynamics/FiniteHilbertTransform.jl.git")'
+```
 
 ---
 ## Quickstart

diff --git a/examples/plasmademo.png b/examples/plasmademo.png
diff --git a/examples/run_plasma.jl b/examples/run_plasma.jl
@@ -72,53 +72,6 @@ function CGFuncPlasma(uNodes::Vector{Float64}, qSELF::Float64, xmax::Float64)
 end
 
 
-"""
-Compute the values of I-Xi(omg) for a given complex frequency.
-
-# Arguments
-- `omg::Complex{Float64}`: Complex frequency.
-- `taba::Vector{Float64}`: Vector of coefficients a_k(u).
-- `xmax::Float64`: Maximum value of x.
-- `FHT::FiniteHilbertTransform.AbstractFHT`: AbstractFHT structure.
-
-# Returns
-- `IminusXi::Complex{Float64}`: Value of I-Xi(omg).
-
-# Example
-```julia
-GetLegendreIminusXiPlasma(1.0 + 1.0im, taba, 10.0, FHT)
-```
-"""
-function GetLegendreIminusXiPlasma(omg::Complex{Float64},
-                                   taba::Vector{Float64},
-                                   xmax::Float64,
-                                   FHT::FiniteHilbertTransform.AbstractFHT)
-
-
-    # Rescale the COMPLEX frequency
-    K_u = size(taba,1)
-
-    varpi = omg/xmax
-
-    # compute the Hilbert-transformed Legendre functions
-    FiniteHilbertTransform.GettabD!(varpi,FHT)
-
-    xi = 0.0 + 0.0*im # Initialise xi
-
-    # loop over the Legendre functions
-    for k=1:(K_u)
-
-        # add the contribution
-        xi += taba[k]*FHT.tabDLeg[k]
-    end
-
-    # compute 1.0 - xi
-    IminusXi = 1.0 - xi
-
-    return IminusXi # Output
-end
-
-
 
 """
 Compute coefficients a_k(u) for all values of u by looping over Legendre weights, Legendre polynomials, and G(u) values.
@@ -136,7 +89,7 @@ Compute coefficients a_k(u) for all values of u by looping over Legendre weights
 ComputeALegendre(FHT, tabG)
 ```
 """
-function ComputeALegendre(FHT::FiniteHilbertTransform.LegendreFHT,tabG::Vector{Float64})
+function ComputeA(FHT::FiniteHilbertTransform.AbstractFHT,tabG::Vector{Float64})
 
     taba, warnflag = FiniteHilbertTransform.GetaXi(FHT,tabG)
 
@@ -164,7 +117,7 @@ ComputeIminusXi([1.0+1.0im, 2.0+2.0im], [0.1, 0.2, 0.3], 10.0, [struct_tabLeg_1,
 function ComputeIminusXi(tabomega::Vector{Complex{Float64}},
                          taba::Vector{Float64},
                          xmax::Float64,
-                         struct_tabLeg::Vector{FiniteHilbertTransform.LegendreFHT})
+                         struct_tabFHT)
 
     # get constants
     K_u    = size(taba,1)
@@ -180,7 +133,7 @@ function ComputeIminusXi(tabomega::Vector{Complex{Float64}},
         thr = Threads.threadid()
 
         # compute I-Xi(omg) using the parallel containers
-        val = GetLegendreIminusXiPlasma(tabomega[iomega],taba,xmax,struct_tabLeg[thr])
+        val = FiniteHilbertTransform.GetIminusXi(tabomega[iomega]/xmax,taba,struct_tabFHT[thr])
 
         # fill in tabIminusXi
         tabIminusXi[iomega] = val
@@ -218,7 +171,24 @@ function setup_legendre_integration(Ku::Int64, qself::Float64, xmax::Float64, PA
     tabG = CGFuncPlasma(FHT.tabu, qself, xmax)
 
     # Compute the coefficients for integration
-    taba, warnflag = ComputeALegendre(FHT, tabG)
+    taba, warnflag = ComputeA(FHT, tabG)
+
+    # Set up the table for integration
+    FHTlist = [deepcopy(FHT) for k = 1:Threads.nthreads()]
+
+    return taba, FHTlist
+end
+
+
+function setup_chebyshev_integration(Ku::Int64, qself::Float64, xmax::Float64, PARALLEL::Bool=false)
+    # Filling in the arrays used in the Chebyshev quadrature
+    FHT = FiniteHilbertTransform.ChebyshevFHT(Ku)
+
+    # Compute the function G(u)
+    tabG = CGFuncPlasma(FHT.tabu, qself, xmax)
+
+    # Compute the coefficients for integration
+    taba, warnflag = ComputeA(FHT, tabG)
 
     # Set up the table for integration
     FHTlist = [deepcopy(FHT) for k = 1:Threads.nthreads()]
@@ -228,6 +198,8 @@ end
 
 
 
+
+
 """
     get_tabomega(tabOmega::Vector{Float64}, tabEta::Vector{Float64})
 
@@ -319,7 +291,7 @@ function parse_commandline()
         "--Etamax"
             help     = "Maximum imaginary frequency"
             arg_type = Float64
-            default  = -0.01
+            default  = 3.0
         "--Cmode"
             help     = "Continuation mode for damped calculation (legendre/chebyshev,rational)"
             arg_type = String
@@ -364,6 +336,7 @@ function main()
     qself    = parsed_args["qSELF"]
     xmax     = parsed_args["xmax"]
     K_u      = parsed_args["K_u"]
+    CMODE    = parsed_args["Cmode"]
 
 
     if (PARALLEL)
@@ -379,29 +352,24 @@ function main()
     # (flat) array of omega values to check
     tabomega = get_tabomega(tabOmega,tabEta)
 
-    taba,struct_tabLeg = setup_legendre_integration(K_u,qself,xmax,PARALLEL)
-    @time tabIminusXi = ComputeIminusXi(tabomega,taba,xmax,struct_tabLeg)
-
-
-    # Prefix of the directory where the files are dumped
-    prefixnamefile = "data/"
-
-    # Name of the file where the data is dumped
-    namefile = prefixnamefile*"data_"*parsed_args["Cmode"]*"_Plasma_Ku_"*string(K_u)*
-               "_qSELF_"*string(qself)*"_xmax_"*string(xmax)*".hf5"
+    if (CMODE=="legendre")
+        taba,struct_tabLeg = setup_legendre_integration(K_u,qself,xmax,PARALLEL)
+    elseif (CMODE=="chebyshev")
+        taba,struct_tabLeg = setup_chebyshev_integration(K_u,qself,xmax,PARALLEL)
+    else
+        throw(DomainError(CMODE,"Integration type must be legendre or chebyshev."))
+    end
 
-    # you can save the data by uncommenting this:
-    #print("Dumping the data | ")
-    #@time FiniteHilbertTransform.dump_tabIminusXi(namefile,tabomega,tabIminusXi) # Dumping the values of det[I-Xi]
+    @time tabIminusXi = ComputeIminusXi(tabomega,taba,xmax,struct_tabLeg)
 
     epsilon_real = reshape(real.(tabIminusXi),parsed_args["nEta"],parsed_args["nOmega"])
     epsilon_imag = reshape(imag.(tabIminusXi),parsed_args["nEta"],parsed_args["nOmega"])
     epsilon = abs.(epsilon_real .+ im * epsilon_imag)
 
 
     # Plot
-    contour(tabOmega,tabEta,log10.(epsilon), levels=10, color=:black, #levels=[-2.0, -1.5, -1.0, -0.5, -0.25, 0.0], 
-            xlabel="Re[ω]", ylabel="Im[ω]", xlims=(-4, 4), ylims=(-3, 0),
+    contour(tabOmega,tabEta,log10.(epsilon), levels=10, color=:black, 
+            xlabel="Re[ω]", ylabel="Im[ω]", xlims=(parsed_args["Omegamin"],parsed_args["Omegamax"]), ylims=(parsed_args["Etamin"],parsed_args["Etamax"]),
             clims=(-2, 0), aspect_ratio=:equal, legend=false)
     savefig("plasmademo.png")
 

diff --git a/src/Chebyshev/Chebyshev.jl b/src/Chebyshev/Chebyshev.jl
@@ -14,12 +14,13 @@ struct ChebyshevFHT <: AbstractFHT
 
     tabu::Array{Float64,1}     # u values (sampling points)
     tabw::Array{Float64,1}     # w values (weights at sampling points)
-    tabP::Matrix{Float64}     # P_k(u) values (Ku x Ku)
+    tabP::Matrix{Float64}     # a values
     tabc::Vector{Float64}     # prefactor at each sampling point
+    taba::Array{Float64,1}
 
     # arrays for the continuation
-    tabPLeg::Array{ComplexF64,1} # Static container for tabPLeg
-    tabQLeg::Array{ComplexF64,1} # Static container for tabQLeg
+    tabTLeg::Array{ComplexF64,1} # Static container for tabTLeg
+    tabULeg::Array{ComplexF64,1} # Static container for tabULeg
     tabDLeg::Array{ComplexF64,1} # Static container for tabDLeg
 
 end
@@ -35,7 +36,7 @@ function ChebyshevFHT(Ku::Int64;name::String="Chebyshev")
 
     tabu,tabw,tabc,tabP = tabCquad(Ku)
 
-    return ChebyshevFHT(name,Ku,tabu,tabw,tabP,tabc,
+    return ChebyshevFHT(name,Ku,tabu,tabw,tabP,tabc,zeros(Float64,Ku),
                                   zeros(ComplexF64,Ku),zeros(ComplexF64,Ku),zeros(ComplexF64,Ku))
 
 end
@@ -59,22 +60,22 @@ function GettabD!(omg::ComplexF64,
             println("FiniteHilbertTransform.Chebyshev.get_Chebyshev_Xi: Using DAMPED Chebyshev integration.")
         end
 
-        get_Xi_DAMPED(omg,FHT.taba)
+        get_Xi_DAMPED(omg,FHT)
 
     elseif (imag(omg) == 0.0)
         if verbose > 2
             println("FiniteHilbertTransform.Chebyshev.get_Chebyshev_Xi: Using NEUTRAL Chebyshev integration.")
         end
 
-        get_Xi_NEUTRAL(omg,FHT.taba)
+        get_Xi_NEUTRAL(omg,FHT)
 
     else
         # by default use unstable integration
         if verbose > 2
             println("FiniteHilbertTransform.Chebyshev.get_Chebyshev_Xi: Using UNSTABLE Chebyshev integration.")
         end
 
-        get_Xi_UNSTABLE(omg,FHT.taba)
+        get_Xi_UNSTABLE(omg,FHT)
     end
 end
 
@@ -190,15 +191,31 @@ include("Damped.jl")
 
 function GetaXi!(FHT::ChebyshevFHT,
                  tabGXi::AbstractVector{Float64},
-                 res::Vector{Float64},warnflag::Vector{Float64})
+                 res::Vector{Float64},warnflag::Int64)
 
      # Perfoming the discrete sine transform
      taba_temp = FFTW.r2r(tabGXi,FFTW.RODFT10,1)
 
      for i=1:FHT.Ku
-         res[i] = taba_temp[i] / FHT.Ku
+        FHT.taba[i] = taba_temp[i] / FHT.Ku
      end
- end
+
+     return FHT.taba,warnflag
+
+end
+
+function GetaXi!(FHT::ChebyshevFHT,
+    tabGXi::AbstractVector{Float64},
+    res::Vector{Float64},warnflag::Vector{Float64})
+
+    println("FiniteHilbertTransform.GetaXi!: deprecation warning: warnflag is now an integer.")
+
+    res,warnval = GetaXi!(FHT,tabG,res,0)
+    warnflag[1] = warnval
+
+    return res,warnflag
+
+end
 
 """
 
@@ -209,12 +226,45 @@ function GetaXi(FHT::ChebyshevFHT,
 
 
     # start with no warnings
-    warnflag = zeros(FHT.Ku)
+    warnflag = 0#zeros(FHT.Ku)
 
     res = zeros(Float64,FHT.Ku)
 
     GetaXi!(FHT,tabGXi,res,warnflag)
 
-    return res,warnflag
+    return FHT.taba,warnflag
 
 end
+
+
+function GetIminusXi(ϖ::Complex{Float64},taba::Vector{Float64},FHT::ChebyshevFHT;verbose::Int64=0)
+
+    if (imag(ϖ) < 0.0)
+        if verbose > 2
+            println("FiniteHilbertTransform.Chebyshev.get_Chebyshev_Xi: Using DAMPED Chebyshev integration.")
+        end
+
+        xi = get_Xi_DAMPED(ϖ,FHT)
+
+    elseif (imag(ϖ) == 0.0)
+        if verbose > 2
+            println("FiniteHilbertTransform.Chebyshev.get_Chebyshev_Xi: Using NEUTRAL Chebyshev integration.")
+        end
+
+        xi = get_Xi_NEUTRAL(ϖ,FHT)
+
+    else
+        # by default use unstable integration
+        if verbose > 2
+            println("FiniteHilbertTransform.Chebyshev.get_Chebyshev_Xi: Using UNSTABLE Chebyshev integration.")
+        end
+
+        xi = get_Xi_UNSTABLE(ϖ,FHT)
+    end
+
+    # compute 1.0 - xi
+    IminusXi = 1.0 - xi
+
+    return IminusXi # Output
+end
+
diff --git a/src/Chebyshev/Damped.jl b/src/Chebyshev/Damped.jl
@@ -1,9 +1,8 @@
 ##################################################
 # Response matrix in the DAMPED case
 ##################################################
-function get_Xi_DAMPED(omg::ComplexF64,
-                       taba::Vector{Float64})
-    sumT, sumU = get_sumT(omg,taba), get_sumU(omg,taba) # Computing the needed sum
+function get_Xi_DAMPED(omg::ComplexF64,struct_tabCheb::ChebyshevFHT)
+    sumT, sumU = get_sumT(omg,struct_tabCheb.taba), get_sumU(omg,struct_tabCheb.taba) # Computing the needed sum
     #####
     Xi = -sumT # Starting to compute the expression
     #####
@@ -20,6 +19,7 @@ function get_Xi_DAMPED(omg::ComplexF64,
     return Xi # Output
 end
 
+"""
 function get_Xi_array(omg::ComplexF64,
                        taba::Vector{Float64})
     sumT, sumU = get_sumT(omg,taba), get_sumU(omg,taba) # Computing the needed sum
@@ -38,3 +38,4 @@ function get_Xi_array(omg::ComplexF64,
     #####
     return Xi # Output
 end
+"""