check chebyshev plasma calculation: good!

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

examples/run_plasma.jl

@@ -89,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)
 
@@ -117,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)
@@ -133,7 +133,7 @@ function ComputeIminusXi(tabomega::Vector{Complex{Float64}},
         thr = Threads.threadid()
 
         # compute I-Xi(omg) using the parallel containers
-        val = FiniteHilbertTransform.GetIminusXi(tabomega[iomega]/xmax,taba,struct_tabLeg[thr])
+        val = FiniteHilbertTransform.GetIminusXi(tabomega[iomega]/xmax,taba,struct_tabFHT[thr])
 
         # fill in tabIminusXi
         tabIminusXi[iomega] = val
@@ -171,7 +171,7 @@ 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()]
@@ -180,6 +180,25 @@ function setup_legendre_integration(Ku::Int64, qself::Float64, xmax::Float64, PA
 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()]
+
+    return taba, FHTlist
+end
+
+
+
+
 
 """
     get_tabomega(tabOmega::Vector{Float64}, tabEta::Vector{Float64})
@@ -272,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
@@ -317,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)
@@ -332,7 +352,14 @@ 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)
+    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
+
     @time tabIminusXi = ComputeIminusXi(tabomega,taba,xmax,struct_tabLeg)
 
     epsilon_real = reshape(real.(tabIminusXi),parsed_args["nEta"],parsed_args["nOmega"])
@@ -341,8 +368,8 @@ function main()
 
 
     # 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")
 

src/FiniteHilbertTransform.jl

@@ -21,7 +21,6 @@ module FiniteHilbertTransform
 
 
 # make an abstract FiniteHilbertTransform type
-# @WARNING: Should be defined before any basis definition
 abstract type AbstractFHT  end
 
 # bring in the generic integration tools

test/runtests.jl

@@ -66,6 +66,21 @@ FiniteHilbertTransform.GettabD!(ϖ,FHT)
     @test imag(FHT.tabPLeg[1]) == 0.0
 end
 
+@testset "Legendre Damped Frequency: Heaviside checks" begin
+    ϖ = 0.99 - 0.02im
+    FiniteHilbertTransform.GettabD!(ϖ,FHT)
+    @test real(FHT.tabDLeg[1]) ≈ -4.48863636973214 atol=tol
+    ϖ = 1.0 - 0.02im
+    FiniteHilbertTransform.GettabD!(ϖ,FHT)
+    @test real(FHT.tabDLeg[1]) ≈ -4.605220183488258 atol=tol
+    ϖ = 1.01 - 0.02im
+    FiniteHilbertTransform.GettabD!(ϖ,FHT)
+    @test real(FHT.tabDLeg[1]) ≈ -4.498635453116724 atol=tol
+    ϖ = -1.01 - 0.02im
+    FiniteHilbertTransform.GettabD!(ϖ,FHT)
+    @test real(FHT.tabDLeg[1]) ≈ 4.498635453116724 atol=tol
+end
+
 # check the quadrature approximation
 # by checking with all ones, we know the analytic value
 @testset "Legendre Quadrature" begin