From 0978cb460d088675a8bacd314c3f5a18ef297310 Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Sun, 9 Jun 2024 21:37:37 +0100 Subject: [PATCH] Improve docs; add tests --- .github/workflows/CI.yml | 41 +++++++++ .github/workflows/documentation.yml | 2 +- Project.toml | 8 +- docs/make.jl | 2 + docs/src/bases.md | 5 ++ docs/src/common.md | 22 +++++ docs/src/index.md | 1 + docs/src/quickstart.md | 58 +++++++++++++ src/AstroBasis.jl | 27 ++++-- src/Basisdoc.jl | 22 ++++- src/CB72.jl | 92 ++++++++++++++------ src/CB73.jl | 128 +++++++++++++++++++++------- src/EXPeof.jl | 27 +++++- src/Hernquist.jl | 54 ++++++++++-- src/Kalnajs76.jl | 52 +++++++++-- src/Utils/IO.jl | 47 +++++++++- test/runtests.jl | 9 ++ test/test_razorthin.jl | 29 +++++++ test/test_spherical.jl | 30 +++++++ 19 files changed, 562 insertions(+), 94 deletions(-) create mode 100644 .github/workflows/CI.yml create mode 100644 docs/src/common.md create mode 100644 docs/src/quickstart.md create mode 100644 test/runtests.jl create mode 100644 test/test_razorthin.jl create mode 100644 test/test_spherical.jl diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml new file mode 100644 index 0000000..7fa51b2 --- /dev/null +++ b/.github/workflows/CI.yml @@ -0,0 +1,41 @@ + +name: CI +on: + - push + +jobs: + test: + name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + version: + - '1.6' + - 'nightly' + os: + - ubuntu-latest + arch: + - x64 + steps: + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v1 + with: + version: ${{ matrix.version }} + arch: ${{ matrix.arch }} + - uses: actions/cache@v1 + env: + cache-name: cache-artifacts + with: + path: ~/.julia/artifacts + key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }} + restore-keys: | + ${{ runner.os }}-test-${{ env.cache-name }}- + ${{ runner.os }}-test- + ${{ runner.os }}- + - uses: julia-actions/julia-buildpkg@v1 + - uses: julia-actions/julia-runtest@v1 + - uses: julia-actions/julia-processcoverage@v1 + - uses: codecov/codecov-action@v4.0.1 + with: + token: ${{ secrets.CODECOV_TOKEN }} \ No newline at end of file diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index 3ec399b..67fa984 100644 --- a/.github/workflows/documentation.yml +++ b/.github/workflows/documentation.yml @@ -3,7 +3,7 @@ name: Documentation on: push: branches: - - documentation # update to match your development branch (master, main, dev, trunk, ...) + - main tags: '*' pull_request: diff --git a/Project.toml b/Project.toml index 3c0544b..f0f93aa 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "AstroBasis" uuid = "f9f069e8-36d3-4f22-8c0d-26ca7e8540d8" authors = ["Mike Petersen ", "Mathieu Roule "] -version = "0.2.0" +version = "0.3.0" [deps] HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f" @@ -10,3 +10,9 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" [compat] julia = "1.6" # Equivalent to "^1.6", will allow any 1.a.b version with a>=6 + +[extras] +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["Test"] \ No newline at end of file diff --git a/docs/make.jl b/docs/make.jl index c032b7e..6366da2 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -6,6 +6,8 @@ using Documenter, AstroBasis makedocs(sitename = "AstroBasis.jl", pages=[ "Home" => "index.md", + "Quickstart" => "quickstart.md", + "Common Interface" => "common.md", "Bases" => "bases.md" ], format = Documenter.HTML(prettyurls=false)) diff --git a/docs/src/bases.md b/docs/src/bases.md index 9d8c85d..030d7cd 100644 --- a/docs/src/bases.md +++ b/docs/src/bases.md @@ -8,6 +8,11 @@ AstroBasis.CB73Basis ``` +### Hernquist (1992) +```@docs +AstroBasis.HernquistBasis +``` + ## Razor-thin bases diff --git a/docs/src/common.md b/docs/src/common.md new file mode 100644 index 0000000..5120074 --- /dev/null +++ b/docs/src/common.md @@ -0,0 +1,22 @@ +# Functions common to all bases + +To seamlessly allow for switching between different bases, the bases are all organised around common functions. + +## Accessing potential functions +```@docs +AstroBasis.getUln +``` + +```@docs +AstroBasis.tabUl! +``` + +## Accessing density functions + +```@docs +AstroBasis.getDln +``` + +```@docs +AstroBasis.tabDl! +``` diff --git a/docs/src/index.md b/docs/src/index.md index f88959a..913f118 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -2,3 +2,4 @@ *Bi-orthogonal bases for galactic dynamics in Julia* +AstroBasis.jl implements several common basis sets (see 'Bases') in galactic dynamics, with a common interface (see 'Common Interfaces'). A quickstart is availabe, which will walk you through how to construct an image of the basis sets. diff --git a/docs/src/quickstart.md b/docs/src/quickstart.md new file mode 100644 index 0000000..2e16390 --- /dev/null +++ b/docs/src/quickstart.md @@ -0,0 +1,58 @@ +# Quickstart Example + +In this example code, we will make a figure of the radial basis elements from Clutton-Brock (1973). + +```julia +import AstroBasis +using Plots +using LaTeXStrings +``` + +Now create the basis: + +```julia +println("Creating the basis ... ") +G, rb = 1., 1. +ltest, nradial = 2, 5 +basis = AstroBasis.CB73Basis(lmax=ltest,nradial=nradial,G=G, rb=rb) +``` + +Define where to make the basis points: + +```julia +# Points (rescaled radius) +println("Compute basis values ... ") +nx = 200 +rmin, rmax = 0., 3. +tabx = collect(LinRange(rmin/basis.rb,rmax/basis.rb,nx)) +``` + +Use the common function `tabUl!` to fill the table: + +```julia +# Compute the values of the potential basis elements and store them +tabU = zeros(nradial,nx) # Storage for the basis values +for j = 1:nx + # Compute all the basis elements values at a given location r (the result is stored in basis.tabUl) + AstroBasis.tabUl!(basis,ltest,tabx[j]*basis.rb) + # Store them in tabU + for i = 1:nradial + tabU[i,j] = basis.tabUl[i] + end +end +``` + +And finally plot: + +```julia +# Plot the curves +println("Plotting ... ") +labels = reshape(["n="*string(k) for k=1:nradial],(1,nradial)) #Need to be row +plU=plot(tabx, transpose(tabU), title = "Potential basis elements: Clutton-Brock (1973)",label=labels) +xlabel!(plU, L"$r / r_{\mathrm{b}}$") +ylabel!(plU, L"$U^{\ell}_n (r)\quad \ell=$"*string(ltest)) +savefig(plU,"CluttonBrock73.png") +println("The plot has been saved in the same folder as this example script under the name 'CluttonBrock73.png'.") +``` + +![`Clutton-Brock (1973)`](../../examples/CluttonBrock73_original.png) \ No newline at end of file diff --git a/src/AstroBasis.jl b/src/AstroBasis.jl index 43960fd..104066a 100644 --- a/src/AstroBasis.jl +++ b/src/AstroBasis.jl @@ -2,31 +2,44 @@ module AstroBasis using HDF5 -# make an abstract Basis type -# @WARNING: Should be defined before any basis definition -abstract type AbstractAstroBasis end +# define an abstract Basis type; define a function to get the dimension for each basis type +abstract type AbstractAstroBasis end +dimension(basis::AbstractAstroBasis) = dimension(typeof(basis)) + + +# split Basis types by dimension; specify the dimensionality for safety; correctly dispatch dimension for subtypes +abstract type RazorThinBasis <: AbstractAstroBasis end +dimension(::Type{<:RazorThinBasis}) = 2 + +abstract type SphericalBasis <: AbstractAstroBasis end +dimension(::Type{<:SphericalBasis}) = 3 + +# make dimension available +export dimension # Documentation of the function needed in a basis include("Utils/IO.jl") +export getparameters # Documentation of the function needed in a basis include("Basisdoc.jl") +export getUln,getDln,tabUl!,tabDl! # bring in the Clutton-Brock (1973) spherical basis include("CB73.jl") -# @IMPROVE, add Bessel basis - # bring in the Hernquist & Ostriker (1992) disc basis include("Hernquist.jl") +# @IMPROVE, add Bessel basis + +# @IMPROVE, add readers for EXP empirical bases + # bring in the Clutton-Brock (1972) disc basis include("CB72.jl") # bring in the Kalnajs (1976) disc basis include("Kalnajs76.jl") -# @IMPROVE, add readers for EXP empirical bases - end # module diff --git a/src/Basisdoc.jl b/src/Basisdoc.jl index c688a0c..275f1a2 100644 --- a/src/Basisdoc.jl +++ b/src/Basisdoc.jl @@ -60,13 +60,27 @@ end """ - tabUl!(basis,l) + tabUl!(basis, l, r) -Compute potential table for a given l and r, and for 0 <= n <= nmax +Compute the basis potential elements for a given basis structure. -To avoid repeated memory allocations, this overwrites the 'basis.tabUl' table +# Arguments +- `basis`: The structure containing basis parameters and prefactors. This should be a subtype of `AbstractAstroBasis`. +- `l::Int64`: The harmonic/azimuthal index. +- `r::Float64`: The radial coordinate. + +# Description +This function computes the basis elements for the provided structure. It updates the +potential or density basis elements in-place within the given structure. + +# Implementation Details +- The function initializes the recurrence relation for the adimensional subpart of the basis elements. +- It then iterates through the radial basis numbers, updating the basis elements using a recurrence relation. +- Finally, it multiplies the computed values by the appropriate prefactors. + +# Warning +> l=0 is index 1, n=0 is index 1. -@WARNING:: l=0 is index 1, n=0 is index 1. """ function tabUl!(basis,l,r) # ... [implementation sold separately] ... diff --git a/src/CB72.jl b/src/CB72.jl index b5d0e1c..ea0eed2 100755 --- a/src/CB72.jl +++ b/src/CB72.jl @@ -9,14 +9,28 @@ radial index (as lmax for azimuthal/harmonic index l). """ """ - CB72Basis - -Radial basis elements from Clutton-Brock (1972) + CB72Basis <: RazorThinBasis + +A structure for interfacing with the radial basis elements of Clutton-Brock (1972). + +# Fields +- `name::String`: Basis name (default CB72). +- `lmax::Int64`: Maximal harmonic/azimuthal index (starts at 0). +- `nradial::Int64`: Number of radial basis elements (≥ 1). +- `G::Float64`: Gravitational constant (default 1.0). +- `rb::Float64`: Radial extension (default 1.0). +- `tabPrefU::Array{Float64,2}`: Potential prefactors array. +- `tabPrefD::Array{Float64,2}`: Density prefactors array. +- `tabUl::Array{Float64,1}`: Potential elements value array. +- `tabDl::Array{Float64,1}`: Density elements value array. + +# Description +The `CB72Basis` structure is the interface to the radial basis elements as defined +by Clutton-Brock (1972). """ -struct CB72Basis <: AbstractAstroBasis +struct CB72Basis <: RazorThinBasis name::String # Basis name (default CB72) - dimension::Int64 # Basis dimension (default 2) lmax::Int64 # Maximal harmonic/azimuthal index (starts at 0) nradial::Int64 # Number of radial basis elements (≥ 1) @@ -32,22 +46,23 @@ struct CB72Basis <: AbstractAstroBasis end + """ CB72Basis([name, dimension, lmax, nradial, G, rb]) creates a CB72Basis structure (and fill prefactors) """ -function CB72Basis(;name::String="CB72", dimension::Int64=2, +function CB72Basis(;name::String="CB72", lmax::Int64=0, nradial::Int64=0, G::Float64=1., rb::Float64=1.) - basis = CB72Basis(name,dimension, + basis = CB72Basis(name, lmax,nradial, G,rb, zeros(Float64,nradial,lmax+1),zeros(Float64,nradial,lmax+1), # Prefactors arrays zeros(Int64,nradial),zeros(Float64,nradial)) # Elements value arrays - FillPrefactors!(basis) + _FillPrefactors!(basis) return basis end @@ -57,7 +72,7 @@ end # Computation of the prefactors ################################################## """ - PrefactorsAzimuthalRecurrenceCB72(previous, l, n) + _PrefactorsAzimuthalRecurrenceCB72(previous, l, n) Gives the next adimensional prefactor (recurrence over azimuthal number l), a^n_l = a^n_{l-1} * 2 * sqrt( 1 / [(n+2l)*(n+2l-1)] ), @@ -66,14 +81,14 @@ given the previous one (a^n_{l-1}). Initializing the recurrence at a^n_0 = sqrt(2), one gets the expected prefactor a^n_l = 2^(l+1/2) * sqrt( n! / (n+2l)! ). """ -function PrefactorsAzimuthalRecurrenceCB72(previous::Float64, +function _PrefactorsAzimuthalRecurrenceCB72(previous::Float64, l::Int64, n::Int64)::Float64 return 2.0 * sqrt( 1.0 / ( (n+2.0*l)*(n+2.0*l-1.0) ) ) * previous end """ - FillPrefactors!(basis::CB72Basis) + _FillPrefactors!(basis::CB72Basis) Clutton-Brock (1972) prefactors a^n_l = 2^(l+1/2) * sqrt( n! / (n+2l)! ) computed through the recurrence relation given by the function @@ -81,7 +96,7 @@ prefactors_recurrence_l_CB72. @IMPROVE precompute the a^n_l up to a given limit and save them in a file to read? """ -function FillPrefactors!(basis::CB72Basis) +function _FillPrefactors!(basis::CB72Basis) lmax, nradial = basis.lmax, basis.nradial G, rb = basis.G, basis.rb @@ -99,8 +114,8 @@ function FillPrefactors!(basis::CB72Basis) # Recurrence for n=0:nmax # Loop over the radial basis numbers. ATTENTION, index starts at n=0 for l=1:lmax # Loop over the harmonic indices. ATTENTION, harmonic index starts at l=0 (here 0 has been initialized) - tabPrefU[n+1,l+1] = PrefactorsAzimuthalRecurrenceCB72(tabPrefU[n+1,l],l,n) - tabPrefD[n+1,l+1] = PrefactorsAzimuthalRecurrenceCB72(tabPrefD[n+1,l],l,n) + tabPrefU[n+1,l+1] = _PrefactorsAzimuthalRecurrenceCB72(tabPrefU[n+1,l],l,n) + tabPrefD[n+1,l+1] = _PrefactorsAzimuthalRecurrenceCB72(tabPrefD[n+1,l],l,n) end end @@ -111,23 +126,23 @@ end # Computation of the potential basis elements ################################################## """ - ρCB72(x) + _ρCB72(x) returns the parameter -1 <= ρ <= 1 for a given dimensionless radius x=r/rb. """ -@inline function ρCB72(x::Float64)::Float64 +@inline function _ρCB72(x::Float64)::Float64 return (x*x - 1.0)/(x*x + 1.0) # Value of rho end """ - UAzimuthalInitializationCB72(x, l) + _UAzimuthalInitializationCB72(x, l) Gives the (n=0, l) adimensional potential basis element value at a given position x = r/rb, xi^0_l(x) = prod(2k - 1, k=1 -> l) / [ (1 + x^2)^{l + 1/2} ]. @IMPROVE other function to give the initialization for all l simultaneously through recurrence? """ -function UAzimuthalInitializationCB72(x::Float64, +function _UAzimuthalInitializationCB72(x::Float64, l::Int64)::Float64 return prod(1:2:(2*l-1)) / (sqrt(1.0 + x*x) * (1.0 + x*x)^(l)) end @@ -141,7 +156,7 @@ xi^n_l(x), given the 2 previous ones, xi^{n-1}_l(x) and xi^{n-2}_l(x). Initializing the recurrence at xi^0_l given by the initialization from U_init_l_CB72, one get the right next adimensional potential basis element from Clutton-Brock (1972). """ -function URadialRecurrenceCB72(u0::Float64, u1::Float64, +function _URadialRecurrenceCB72(u0::Float64, u1::Float64, ρ::Float64, l::Int64, n::Int64)::Float64 return ( 2.0 + ((2.0*l-1.0)/n) ) * ρ * u1 - ( 1.0 + ((2.0*l-1.0)/n) ) * u0 @@ -150,7 +165,32 @@ end """ tabUl!(basis::CB72Basis, l, r[, forD]) -For Clutton-Brock (1972) basis elements. +Compute the Clutton-Brock (1972) basis elements for a given `CB72Basis` structure. + +# Arguments +- `basis::CB72Basis`: The `CB72Basis` structure containing basis parameters and prefactors. +- `l::Int64`: The harmonic/azimuthal index. +- `r::Float64`: The radial coordinate. +- `forD::Bool=false`: A boolean flag indicating whether to use density prefactors. Default is `false`. + +# Description +This function computes the basis elements for the Clutton-Brock (1972) model. It updates the +potential or density basis elements in-place within the provided `CB72Basis` structure. + +# Implementation Details +- The function initializes the recurrence relation for the adimensional subpart of the basis elements. +- It then iterates through the radial basis numbers, updating the basis elements using a recurrence relation. +- Finally, it multiplies the computed values by the appropriate prefactors. + +# Example +```julia +# Assuming basis is an instance of CB72Basis with appropriate parameters +tabUl!(basis, 2, 1.0) +``` + +# Warning +> Ensure that the tabPrefU and tabPrefD arrays in CB72Basis are correctly sized to avoid runtime errors. + """ function tabUl!(basis::CB72Basis, l::Int64,r::Float64, @@ -165,17 +205,17 @@ function tabUl!(basis::CB72Basis, x = r/rb # Dimensionless radius xl = x^(l) - ρ = ρCB72(x) # Value of the rescaled parameter rho + ρ = _ρCB72(x) # Value of the rescaled parameter rho ##### # Recurrence on adimensional subpart (xi) of the potential basis elements ##### # Initialization - u0, u1 = 0.0, UAzimuthalInitializationCB72(x,l) # n = -1, 0 + u0, u1 = 0.0, _UAzimuthalInitializationCB72(x,l) # n = -1, 0 tabl[1] = u1 # Recurrence loop for n=1:nmax # Loop over the radial basis numbers. ATTENTION, index starts at n=0 - u2 = URadialRecurrenceCB72(u0,u1,ρ,l,n) + u2 = _URadialRecurrenceCB72(u0,u1,ρ,l,n) tabl[n+1] = u2 u0, u1 = u1, u2 end @@ -203,16 +243,16 @@ function getUln(basis::CB72Basis, tabPrefU = basis.tabPrefU x = r/rb # Dimensionless radius - ρ = ρCB72(x) # Value of the rescaled parameter rho + ρ = _ρCB72(x) # Value of the rescaled parameter rho ##### # Recurrence on adimensional subpart (xi) of the potential basis elements ##### # Initialization - uprev, uact = 0.0, UAzimuthalInitializationCB72(x,l) # n = -1, 0 + uprev, uact = 0.0, _UAzimuthalInitializationCB72(x,l) # n = -1, 0 # Recurrence loop for np=1:n # Loop over the radial basis numbers. ATTENTION, index starts at n=0 - unext = URadialRecurrenceCB72(uprev,uact,ρ,l,np) + unext = _URadialRecurrenceCB72(uprev,uact,ρ,l,np) uprev, uact = uact, unext end diff --git a/src/CB73.jl b/src/CB73.jl index 727cb1e..c5d4cf2 100644 --- a/src/CB73.jl +++ b/src/CB73.jl @@ -10,14 +10,55 @@ data_path_cb73() = abspath(joinpath(@__DIR__, "tables", "data_CB73_lmax_50_nmax_ """ - CB73Basis + CB73Basis <: SphericalBasis + +A structure for accessing the radial basis elements of Clutton-Brock (1973). + +# Fields +- `name::String`: Basis name (default CB73). +- `lmax::Int64`: Maximal harmonic/azimuthal index (starts at 0). +- `nradial::Int64`: Number of radial basis elements (≥ 1). +- `G::Float64`: Gravitational constant (default 1.0). +- `rb::Float64`: Radial extension (default 1.0). +- `tabPrefU::Array{Float64,2}`: Potential prefactors array. +- `tabPrefD::Array{Float64,2}`: Density prefactors array. +- `tabUl::Array{Float64,1}`: Potential elements value array. +- `tabDl::Array{Float64,1}`: Density elements value array. + +# Description +The `CB73Basis` structure is used to access the radial basis elements as defined +by Clutton-Brock (1973). + +# Example +```julia +using Random + +# Creating an instance of CB73Basis +cb73 = CB73Basis( + "CB73", # name + 10, # lmax + 5, # nradial + 1.0, # G + 1.0, # rb + rand(10, 5), # tabPrefU + rand(10, 5), # tabPrefD + rand(5), # tabUl + rand(5) # tabDl +) + +# Accessing fields +println(cb73.name) +println(cb73.lmax) +println(cb73.G) + +# Warning +> Ensure that the prefactor arrays (tabPrefU and tabPrefD) and the element value arrays +(tabUl and tabDl) have compatible dimensions to avoid runtime errors. -Radial basis elements from Clutton-Brock (1973) """ -struct CB73Basis <: AbstractAstroBasis +struct CB73Basis <: SphericalBasis name::String # Basis name (default CB73) - dimension::Int64 # Basis dimension (default 2) lmax::Int64 # Maximal harmonic/azimuthal index (starts at 0) nradial::Int64 # Number of radial basis elements (≥ 1) @@ -37,16 +78,37 @@ end """ CB73Basis([name, dimension, lmax, nradial, G, rb, filename]) -creates a CB73Basis structure (and fill prefactors) +Create a CB73Basis structure and fill prefactors. + +# Arguments +- `name::String="CB73"`: Basis name. +- `lmax::Int64=0`: Maximal harmonic/azimuthal index (starts at 0). +- `nradial::Int64=1`: Number of radial basis elements (≥ 1). +- `G::Float64=1.0`: Gravitational constant. +- `rb::Float64=1.0`: Radial extension. +- `filename::String=data_path_cb73()`: File path to CB73 prefactors data. + +# Description +The `CB73Basis` constructor function creates a new instance of the `CB73Basis` structure. +It initializes the basis with provided parameters and fills the potential and density prefactor arrays +by reading from a file specified by `filename`. + +# Example +```julia +# Creating a CB73Basis with default parameters +cb73_basis = CB73Basis() + +# Creating a CB73Basis with custom parameters +custom_basis = CB73Basis(name="CustomBasis", lmax=10, nradial=5, G=2.0, rb=2.0) """ -function CB73Basis(;name::String="CB73", dimension::Int64=3, +function CB73Basis(;name::String="CB73", lmax::Int64=0, nradial::Int64=1, G::Float64=1., rb::Float64=1., filename::String=data_path_cb73()) - tabPrefU, tabPrefD = ReadFillCB73Prefactors(lmax,nradial,rb,G,filename) + tabPrefU, tabPrefD = _ReadFillCB73Prefactors(lmax,nradial,rb,G,filename) - basis = CB73Basis(name,dimension, + basis = CB73Basis(name, lmax,nradial, G,rb, tabPrefU,tabPrefD, # Prefactors arrays @@ -55,14 +117,14 @@ function CB73Basis(;name::String="CB73", dimension::Int64=3, return basis end -"""ReadFillCB73Prefactors(lmax,nradial[,rb,G,precomputed_filename]) +"""_ReadFillCB73Prefactors(lmax,nradial[,rb,G,precomputed_filename]) reads a table of pre-computed prefactors for Gegenbauer functions comes loaded with a pre-computed large table of prefactors (probably more than you need!) @WARNING: this may only be run once per session, as it will create constants. """ -function ReadFillCB73Prefactors(lmax::Int64,nradial::Int64, +function _ReadFillCB73Prefactors(lmax::Int64,nradial::Int64, rb::Float64=1.,G::Float64=1., filename::String=data_path_cb73()) @@ -92,20 +154,20 @@ function ReadFillCB73Prefactors(lmax::Int64,nradial::Int64, end -"""rhoCB73(x) +"""_rhoCB73(x) Function that returns the parameter -1 <= rho <= 1 for a given dimensionless radius x=r/Rbasis """ -function rhoCB73(x::Float64) +function _rhoCB73(x::Float64) return (x^(2) - 1.0)/(x^(2) + 1.0) # Value of rho end -"""ClnCB73(alpha,n,rho) +"""_ClnCB73(alpha,n,rho) Definition of the Gegenbauer polynomials These coefficients are computed through an upward recurrence @IMPROVE compute all the basis elements (n)_{1<=n<=nradial} at once """ -function ClnCB73(alpha::Float64,n::Int64,rho::Float64) +function _ClnCB73(alpha::Float64,n::Int64,rho::Float64) v0 = 1.0 # Initial value for n=0 if (n == 0) @@ -129,48 +191,48 @@ function ClnCB73(alpha::Float64,n::Int64,rho::Float64) return v # Output of the value end -"""UlnpCB73(lharmonic,np,r,prefactor_table[,rb]) +"""_UlnpCB73(lharmonic,np,r,prefactor_table[,rb]) Definition of the potential radial basis elements from Clutton-Brock (1973) Be careful: l=0 is index 1, np=0 is index 1. """ -function UlnpCB73(l::Int64,np::Int64,r::Float64,tabPrefCB73_Ulnp::Matrix{Float64},rb::Float64=1.) +function _UlnpCB73(l::Int64,np::Int64,r::Float64,tabPrefCB73_Ulnp::Matrix{Float64},rb::Float64=1.) pref = tabPrefCB73_Ulnp[l+1,np+1] # Value of the prefactor. l,np start l=0 x = r/rb # Dimensionless radius - rho = rhoCB73(x) # Value of the rescaled parameter rho + rho = _rhoCB73(x) # Value of the rescaled parameter rho valR = ((x/(1.0+x^(2)))^(l))/(sqrt(1.0+x^(2))) # Value of the multipole factor - valC = ClnCB73(l+1.0,np,rho) # Value of the Gegenbauer polynomials + valC = _ClnCB73(l+1.0,np,rho) # Value of the Gegenbauer polynomials res = pref*valR*valC # Value of the radial function return res end function getUln(basis::CB73Basis,l::Int64,np::Int64,r::Float64) - return UlnpCB73(l,np,r,basis.tabPrefU,basis.rb) + return _UlnpCB73(l,np,r,basis.tabPrefU,basis.rb) end -"""DlnpCB73(lharmonic,np,r,prefactor_table[,rb]) +"""_DlnpCB73(lharmonic,np,r,prefactor_table[,rb]) Definition of the density radial basis elements from Clutton-Brock (1973) Be careful: l=0 is index 1, np=0 is index 1. """ -function DlnpCB73(l::Int64,np::Int64,r::Float64,tabPrefCB73_Dlnp::Matrix{Float64},rb::Float64=1.) +function _DlnpCB73(l::Int64,np::Int64,r::Float64,tabPrefCB73_Dlnp::Matrix{Float64},rb::Float64=1.) pref = tabPrefCB73_Dlnp[l+1,np+1] # Value of the prefactor. ATTENTION, l starts at l = 0 x = r/rb # Dimensionless radius - rho = rhoCB73(x) # Value of the rescaled parameter rho + rho = _rhoCB73(x) # Value of the rescaled parameter rho valR = ((x/(1.0+x^(2)))^(l))/((1.0+x^(2))^(5/2)) # Value of the multipole factor - valC = ClnCB73(l+1.0,np,rho) # Value of the Gegenbauer polynomials + valC = _ClnCB73(l+1.0,np,rho) # Value of the Gegenbauer polynomials res = pref*valR*valC return res end function getDln(basis::CB73Basis,l::Int64,np::Int64,r::Float64) - return DlnpCB73(l,np,r,basis.tabPrefD,basis.rb) + return _DlnpCB73(l,np,r,basis.tabPrefD,basis.rb) end -"""tabUlnpCB73!(lharmonic,radius,tabUlnp,nradial,prefactor_table[,rb]) +"""tab_UlnpCB73!(lharmonic,radius,tabUlnp,nradial,prefactor_table[,rb]) Compute CB73 potential table for a given l and r, and for 1 <= np <= nradial, nearly simultaneously @@ -178,14 +240,14 @@ To avoid repeated memory allocations, this overwrites a table already in the str @IMPROVE, create outs for nradial<3? """ -function tabUlnpCB73!(l::Int64,r::Float64, +function tab_UlnpCB73!(l::Int64,r::Float64, tabUlnp::Array{Float64,1}, nradial::Int64, tabPrefCB73_Ulnp::Matrix{Float64}, rb::Float64=1.) x = r/rb # Dimensionless radius - rho = rhoCB73(x) # Value of the rescaled parameter rho + rho = _rhoCB73(x) # Value of the rescaled parameter rho valR = ((x/(1.0+x^(2)))^(l))/(sqrt(1.0+x^(2))) # Value of the multipole factor alpha = l+1.0 # Value of alpha, the index of the Gegenbauer polynomials @@ -207,12 +269,12 @@ end function tabUl!(basis::CB73Basis,l::Int64,r::Float64) - tabUlnpCB73!(l,r,basis.tabUl,basis.nradial,basis.tabPrefU,basis.rb) - #return UlnpCB73(l,np,r,basis.tabPrefU,basis.rb) + tab_UlnpCB73!(l,r,basis.tabUl,basis.nradial,basis.tabPrefU,basis.rb) + #return _UlnpCB73(l,np,r,basis.tabPrefU,basis.rb) end -"""tabDlnpCB73!(lharmonic,radius,tabD,nradial,prefactor_table[,rb]) +"""tab_DlnpCB73!(lharmonic,radius,tabD,nradial,prefactor_table[,rb]) Compute CB73 density table for a given l and r, and for 1 <= np <= nradial, nearly simultaneously @@ -220,14 +282,14 @@ To avoid repeated memory allocations, this overwrites a table already in the str @IMPROVE, create outs for nradial<3? """ -function tabDlnpCB73!(l::Int64,r::Float64, +function tab_DlnpCB73!(l::Int64,r::Float64, tabDlnp::Array{Float64,1}, nradial::Int64, tabPrefCB73_Dlnp::Matrix{Float64}, rb::Float64=1.) x = r/rb # Dimensionless radius - rho = rhoCB73(x) # Value of the rescaled parameter rho + rho = _rhoCB73(x) # Value of the rescaled parameter rho valR = ((x/(1.0+x^(2)))^(l))/((1.0+x^(2))^(5/2)) # Value of the multipole factor alpha = l+1.0 # Value of alpha, the index of the Gegenbauer polynomials @@ -249,5 +311,5 @@ end function tabDl!(basis::CB73Basis,l::Int64,r::Float64) - tabDlnpCB73!(l,r,basis.tabDl,basis.nradial,basis.tabPrefD,basis.rb) + tab_DlnpCB73!(l,r,basis.tabDl,basis.nradial,basis.tabPrefD,basis.rb) end diff --git a/src/EXPeof.jl b/src/EXPeof.jl index 1cfc727..00f5798 100644 --- a/src/EXPeof.jl +++ b/src/EXPeof.jl @@ -8,10 +8,9 @@ Radial basis elements derived by SLE solution, pre-tabulated """ -struct EXPeofBasis_type <: AbstractAstroBasis +struct EXPeofBasis_type <: SphericalBasis name::String # Basis name (default CB73) - dimension::Int64 # Basis dimension (default 2) lmax::Int64 # Maximal harmonic/azimuthal index (starts at 0) nradial::Int64 # Number of radial basis elements (≥ 1) @@ -29,12 +28,32 @@ end -#isocache = "/Volumes/External1/Isochrone/.slgrid_sph_isochrone" - # set the data path for the basis prefactor elements data_path() = abspath(joinpath(@__DIR__, "tables", "slgrid_sph_isochrone")) # rename for isochrone specifically + +""" + EmpiricalIsochrone([name, dimension, lmax, nmax, G, rb, filename]) + +creates an Isochrone basis structure derived from a Sturm-Liouville solution (and fill prefactors) +""" +function HernquistBasis(;name::String="Hernquist", + lmax::Int64=0, nradial::Int64=1, + G::Float64=1., rb::Float64=1., + filename::String=data_path_hernquist()) + + tabPrefU, tabPrefD = ReadFillHernquistPrefactors(lmax,nradial,rb,G,filename) + basis = HernquistBasis(name, + lmax,nradial, + G,rb, + tabPrefU,tabPrefD, # Prefactors arrays + zeros(Int64,nradial),zeros(Float64,nradial)) # Elements value arrays + + return basis +end + + ################################################## function iso_potential(r::Float64) bc = 1. diff --git a/src/Hernquist.jl b/src/Hernquist.jl index 46c26d0..977a399 100644 --- a/src/Hernquist.jl +++ b/src/Hernquist.jl @@ -8,14 +8,29 @@ data_path_hernquist() = abspath(joinpath(@__DIR__, "tables", "data_HO92_lmax_50_ """ - HernquistBasis +HernquistBasis <: SphericalBasis + +A structure for interfacing with the radial basis elements from Hernquist & Ostriker (1992). + +# Fields +- `name::String`: Basis name (default Hernquist). +- `lmax::Int64`: Maximal harmonic/azimuthal index (starts at 0). +- `nradial::Int64`: Number of radial basis elements (≥ 1). +- `G::Float64`: Gravitational constant (default 1.0). +- `rb::Float64`: Radial extension (default 1.0). +- `tabPrefU::Array{Float64,2}`: Potential prefactors array. +- `tabPrefD::Array{Float64,2}`: Density prefactors array. +- `tabUl::Array{Float64,1}`: Potential elements value array. +- `tabDl::Array{Float64,1}`: Density elements value array. + +# Description +The `HernquistBasis` structure is used to access the radial basis elements as defined +by Hernquist & Ostriker (1992). -Radial basis elements from Hernquist & Ostriker (1992) """ -struct HernquistBasis <: AbstractAstroBasis +struct HernquistBasis <: SphericalBasis name::String # Basis name (default Hernquist) - dimension::Int64 # Basis dimension (default 2) lmax::Int64 # Maximal harmonic/azimuthal index (starts at 0) nradial::Int64 # Number of radial basis elements (≥ 1) @@ -33,17 +48,40 @@ end """ - HernquistBasis([name, dimension, lmax, nmax, G, rb, filename]) +HernquistBasis([name, dimension, lmax, nmax, G, rb, filename]) + +Create a HernquistBasis structure and fill prefactors. + +# Arguments +- `name::String="Hernquist"`: Basis name. +- `lmax::Int64=0`: Maximal harmonic/azimuthal index (starts at 0). +- `nradial::Int64=1`: Number of radial basis elements (≥ 1). +- `G::Float64=1.0`: Gravitational constant. +- `rb::Float64=1.0`: Radial extension. +- `filename::String=data_path_hernquist()`: File path to Hernquist prefactors data. + +# Description +The `HernquistBasis` constructor function creates a new instance of the `HernquistBasis` structure. +It initializes the basis with provided parameters and fills the potential and density prefactor arrays +by reading from a file specified by `filename`. + +# Example +```julia +# Creating a HernquistBasis with default parameters +hernquist_basis = HernquistBasis() + +# Creating a HernquistBasis with custom parameters +custom_basis = HernquistBasis(name="CustomBasis", lmax=10, nradial=5, G=2.0, rb=2.0) +``` -creates a HernquistBasis structure (and fill prefactors) """ -function HernquistBasis(;name::String="Hernquist", dimension::Int64=3, +function HernquistBasis(;name::String="Hernquist", lmax::Int64=0, nradial::Int64=1, G::Float64=1., rb::Float64=1., filename::String=data_path_hernquist()) tabPrefU, tabPrefD = ReadFillHernquistPrefactors(lmax,nradial,rb,G,filename) - basis = HernquistBasis(name,dimension, + basis = HernquistBasis(name, lmax,nradial, G,rb, tabPrefU,tabPrefD, # Prefactors arrays diff --git a/src/Kalnajs76.jl b/src/Kalnajs76.jl index 6a08cce..716a9e3 100755 --- a/src/Kalnajs76.jl +++ b/src/Kalnajs76.jl @@ -18,14 +18,54 @@ using SpecialFunctions using Memoize """ - K76Basis +K76Basis <: RazorThinBasis + +A structure for interfacing with the radial basis elements of Kalnajs (1976). + +# Fields +- `name::String`: Basis name (default K76). +- `lmax::Int64`: Maximal harmonic/azimuthal index (starts at 0). +- `nradial::Int64`: Number of radial basis elements (≥ 1). +- `G::Float64`: Gravitational constant (default 1.0). +- `rb::Float64`: Radial extension (default 1.0). +- `kKA::Int64`: Basis index specific to Kalnajs (1976). +- `tabPrefU::Array{Float64,2}`: Potential prefactors array. +- `tabPrefD::Array{Float64,2}`: Density prefactors array. +- `tabUl::Array{Float64,1}`: Potential elements value array. +- `tabDl::Array{Float64,1}`: Density elements value array. + +# Description +The `K76Basis` structure is used to access the radial basis elements as defined +by Kalnajs (1976). + +# Example +```julia +using Random + +# Creating an instance of K76Basis +k76 = K76Basis( + "K76", # name + 10, # lmax + 5, # nradial + 1.0, # G + 1.0, # rb + 2, # kKA + rand(10, 5), # tabPrefU + rand(10, 5), # tabPrefD + rand(5), # tabUl + rand(5) # tabDl +) + +# Accessing fields +println(k76.name) +println(k76.lmax) +println(k76.kKA) +``` -Radial basis elements from Kalnajs (1976) """ -struct K76Basis <: AbstractAstroBasis +struct K76Basis <: RazorThinBasis name::String # Basis name (default K76) - dimension::Int64 # Basis dimension (default 2) lmax::Int64 # Maximal harmonic/azimuthal index (starts at 0) nradial::Int64 # Number of radial basis elements (≥ 1) @@ -48,11 +88,11 @@ end creates a K76Basis structure """ -function K76Basis(;name::String="K76", dimension::Int64=2, +function K76Basis(;name::String="K76", lmax::Int64=0, nradial::Int64=1, G::Float64=1., rb::Float64=1., kKA::Int64=1) - basis = K76Basis(name,dimension, + basis = K76Basis(name, lmax,nradial, G,rb,kKA, zeros(Float64,lmax+1,nradial),zeros(Float64,lmax+1,nradial), # Prefactors arrays diff --git a/src/Utils/IO.jl b/src/Utils/IO.jl index d42bfd0..9027326 100644 --- a/src/Utils/IO.jl +++ b/src/Utils/IO.jl @@ -33,18 +33,57 @@ end """ GetParameters(basis) -gives all the basis parameters in a dictionnary +Returns all the basis parameters in a dictionary. The dictionary includes only the fields +that are of type `String` or `Number`, ignoring arrays and other types. + +# Arguments +- `basis::AB`: An instance of a type that is a subtype of `AbstractAstroBasis`. + +# Returns +- A dictionary where the keys are the field names as strings and the values are the corresponding field values + if they are of type `String` or `Number`. + +# Example +```julia +# Example basis type +struct ExampleBasis <: SphericalBasis + name::String + lmax::Int64 + nradial::Int64 + G::Float64 + rb::Float64 + tabPrefU::Array{Float64,2} + tabPrefD::Array{Float64,2} + tabUl::Array{Float64,1} + tabDl::Array{Float64,1} +end + +# Creating an instance of ExampleBasis +basis_instance = ExampleBasis("Example", 10, 5, 1.0, 1.0, rand(10, 5), rand(10, 5), rand(10), rand(10)) + +# Getting parameters +params = GetParameters(basis_instance) +println(params) """ -function GetParameters(basis::AB) where {AB <: AbstractAstroBasis} +function getparameters(basis::AB) where {AB <: AbstractAstroBasis} + # Initialize an empty dictionary to store the parameters paramsdict = Dict{String,Union{String,Number}}() + + # Iterate through each field in the basis type for i = 1:fieldcount(AB) - # If this field is a string or a number, it is a parameter - # (Prevent from dumping arrays) + # Check if the field type is a String or Number if fieldtype(AB,i) <: Union{String,Number} + # Convert the field name to a string varname = string(fieldname(AB,i)) + # Get the field value from the basis instance and store it in the dictionary paramsdict[varname] = getfield(basis,i) end end + + # Add the dimensionality parameter + paramsdict["dimension"] = dimension(basis) + + # Return the dictionary of parameters return paramsdict end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl new file mode 100644 index 0000000..9377ebc --- /dev/null +++ b/test/runtests.jl @@ -0,0 +1,9 @@ +using AstroBasis + +using Test + +# test spherical bases +include("test_spherical.jl") + +# test razor thin disc bases +include("test_razorthin.jl") \ No newline at end of file diff --git a/test/test_razorthin.jl b/test/test_razorthin.jl new file mode 100644 index 0000000..ed91490 --- /dev/null +++ b/test/test_razorthin.jl @@ -0,0 +1,29 @@ + +G, rb = 1., 1. +ltest, nradial = 2, 5 + + +@testset "razorthinbases" begin + @testset "CB72" begin + basis = AstroBasis.CB72Basis(lmax=ltest,nradial=nradial,G=G, rb=rb) + @test dimension(basis) == 2 + @test getparameters(basis)["name"] == "CB72" + @test getDln(basis,ltest,nradial-1,rb) ≈ 0.33126698841066865 atol=1e-6 + @test getUln(basis,ltest,nradial-1,rb) ≈ -0.3202172114362374 atol=1e-6 + tabUl!(basis,0,rb) + @test basis.tabUl[1] == -1.0 + tabDl!(basis,0,rb) + @test basis.tabDl[2] == 0 + end + @testset "Kalnajs76" begin + basis = AstroBasis.K76Basis(lmax=ltest,nradial=nradial,G=G, rb=rb) + @test dimension(basis) == 2 + @test getparameters(basis)["name"] == "K76" + @test getDln(basis,ltest,nradial-1,rb) == 0.0 + @test getUln(basis,ltest,nradial-1,rb) ≈ -0.3167679231608087 atol=1e-6 + tabUl!(basis,0,rb) + @test basis.tabUl[1] == -0.8580855308097834 + tabDl!(basis,0,rb) + @test basis.tabDl[1] == 0 + end +end \ No newline at end of file diff --git a/test/test_spherical.jl b/test/test_spherical.jl new file mode 100644 index 0000000..ad0ff98 --- /dev/null +++ b/test/test_spherical.jl @@ -0,0 +1,30 @@ + + +G, rb = 1., 1. +ltest, nradial = 2, 5 + + +@testset "sphericalbases" begin + @testset "CB73" begin + basis = AstroBasis.CB73Basis(lmax=ltest,nradial=nradial,G=G, rb=rb) + @test dimension(basis) == 3 + @test getparameters(basis)["name"] == "CB73" + @test getDln(basis,ltest,nradial-1,rb) ≈ 1.6230683210206467 atol=1e-6 + @test getUln(basis,ltest,nradial-1,rb) ≈ -0.418381088294792 atol=1e-6 + tabUl!(basis,0,rb) + @test basis.tabUl[2] == 0 + tabDl!(basis,0,rb) + @test basis.tabDl[2] == 0 + end + @testset "Hernquist" begin + basis = AstroBasis.HernquistBasis(lmax=ltest,nradial=nradial,G=G, rb=rb) + @test dimension(basis) == 3 + @test getparameters(basis)["name"] == "Hernquist" + @test getDln(basis,ltest,nradial-1,rb) ≈ 1.552003244163087 atol=1e-6 + @test getUln(basis,ltest,nradial-1,rb) ≈ -0.8668021315929388 atol=1e-6 + tabUl!(basis,0,rb) + @test basis.tabUl[2] == 0 + tabDl!(basis,0,rb) + @test basis.tabDl[2] == 0 + end +end \ No newline at end of file