update documentation

Author: fgerick

README.md

@@ -8,16 +8,45 @@
 
 A Julia wrapper around the [SHTns](https://bitbucket.org/nschaeff/shtns) spherical harmonic transform library. 
 
-# This is WORK IN PROGRESS and will certainly have BREAKING changes.
-
-For the moment not a lot is tested. There are wrappers generated by Clang.jl to call most functions exported by the C library. The *more Julian* API is under development. The documentation will come, for the time being I refer to the [documentation of the C library](https://nschaeff.bitbucket.io/shtns/). 
+> [!WARNING]
+> This is WORK IN PROGRESS and will certainly have BREAKING changes.
 
 # Installation
 
 ```julia
 import Pkg; Pkg.add("SHTns")
 ```
 
+# Example usage
+```julia
+using SHTns
+
+cfg = SHTnsCfg(64) #define configuration
+
+θ,ϕ = SHTns.grid(cfg) #get (co-)latitude and longitude
+d = @. 0.3*sin(ϕ') * sin(θ) #create some spatial data (0.3 Y₁¹)
+
+q = SHTns.analys(cfg,d) #transform to spectral space
+q[SHTns.LM(cfg, 1,1)] ≈ 0.3*sqrt(2π/3)*im #true
+
+d2 = SHTns.synth(cfg,q) #transform back to spatial space
+d2 ≈ d #true
+```
+
+Or we can do the same on a CUDA-enabled GPU!
+```julia
+using CUDA
+
+cfg_gpu = SHTnsCfg(64; shtype=SHTns.QuickInit(;gpu=true))
+
+d_gpu = CuArray(d)
+
+q_gpu = SHTns.analys(cfg_gpu,d_gpu)
+Array(q_gpu)[SHTns.LM(cfg, 1,1)] ≈ 0.3*sqrt(2π/3)*im
+
+d2_gpu = SHTns.synth(cfg_gpu,q_gpu)
+(d_gpu ≈ d2_gpu) && (Array(d2_gpu) ≈ d) #true
+```
 
 # License 
 

docs/Project.toml

@@ -1,2 +1,3 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+SHTns = "5dc4ea12-f621-4038-8458-a22eae100703"

docs/make.jl

@@ -5,6 +5,8 @@ makedocs(;
     # format=Documenter.HTML(),
     pages=[
         "Home" => "index.md",
+        "Quickstart" => "quickstart.md",
+        "Advanced usage" => "advanced.md",
         "API" => "api.md",
     ],
     repo="https://github.com/fgerick/SHTns.jl/blob/{commit}{path}#L{line}",

docs/src/advanced.md

@@ -0,0 +1,108 @@
+# Advanced usage
+
+## Configuring the transform
+
+The [`SHTnsCfg`](@ref) configuration object takes several keywords allowing the user to tweak the spherical harmonic transform.
+
+The simplest configuration takes only the maximum spherical harmonic degree `lmax` as an argument
+```@example advanced
+using SHTns
+
+lmax = 100
+cfg = SHTnsCfg(lmax)
+```
+
+We can also define the grid resolution (`nlat` and `nphi`) first, and then `lmax`, as well as the maximum azimuthal degree `mmax`, as well as the _resolution_ in ``m`` through `mres`.
+```@example advanced
+nlat = 100
+nphi = 2nlat+1
+lmax = 2nlat÷3
+mmax = 5
+mres = 1
+cfg = SHTnsCfg(lmax,mmax,mres,nlat,nphi)
+```
+
+The transform type [`SHTnsType`](@ref) can be varied through the `shtype` keyword
+```@example advanced
+cfg = SHTnsCfg(lmax; shtype = SHTns.Gauss())
+```
+
+See also the documentation of the C library on [shtns_type](https://nschaeff.bitbucket.io/shtns/shtns_8h.html#abdccbb8fbce176dbe189d494c94f0f7b) and [grid layouts](https://nschaeff.bitbucket.io/shtns/spat.html).
+
+We can also change the normalization of the spherical harmonics, for example to the Schmidt semi-normalization without [Condon-Shortley phase](https://mathworld.wolfram.com/Condon-ShortleyPhase.html) (default enabled)
+```@example advanced
+cfg = SHTnsCfg(lmax; norm = SHTns.Schmidt(; cs_phase=false))
+```
+See also [Spherical Harmonics storage and normalization](https://nschaeff.bitbucket.io/shtns/spec.html) of the C library.
+
+
+`SHTns` also supports complex to complex transforms (instead of real spatial space to complex spectral space, which is the default).
+
+```@example advanced
+cfg = SHTnsCfg(lmax; transform=Complex)
+```
+
+## GPU
+
+Using `SHTns` on a CUDA-enabled GPU is straightforward.
+Simply install and import [`CUDA.jl`](https://github.com/JuliaGPU/CUDA.jl):
+
+```julia
+import Pkg; Pkg.add("CUDA")
+using CUDA
+```
+
+Then, configure your [`SHTnsCfg`](@ref) using an [`SHTnsType`](@ref) with the keyword `gpu=true`.
+```julia
+cfg_gpu = SHTnsCfg(64; shtype=SHTns.QuickInit(;gpu=true))
+```
+
+The rest remains almost identical to the CPU transform
+```julia
+θ,ϕ = SHTns.grid(cfg) #get (co-)latitude and longitude
+d = @. 0.3*sin(ϕ') * sin(θ) #create some spatial data (0.3 Y₁¹) - on CPU!
+
+d_gpu = CuArray(d) #send spatial data to GPU
+
+q_gpu = SHTns.analys(cfg_gpu,d_gpu) #transform to spectral space
+
+#scalar indexing not allowed on GPU arrays, therefore call Array(q_gpu) to send back to CPU
+y11_coeff = Array(q_gpu)[SHTns.LM(cfg, 1,1)] 
+y11_coeff ≈ 0.3*sqrt(2π/3)*im #true
+
+d2_gpu = SHTns.synth(cfg_gpu,q_gpu) #transform backto spatial space
+d2_gpu ≈ d_gpu #true
+```
+
+
+> [!NOTE]
+> ROCm GPUs may be supported in future releases of SHTns.jl (already available in the C library).
+
+## Batched transform
+
+`SHTns` supports batched transforms, i.e. multiple spherical harmonic transforms together. For example transforms on spherical surfaces at different radii at the same time.
+
+One can exploit this feature using the keyword `howmany`, to give the configuration the number of transforms to be performed.
+
+```@example advanced
+howmany = 10
+cfg = SHTnsCfg(lmax; howmany)
+
+θ,ϕ = SHTns.grid(cfg)
+
+# create spatial array (nlat_padded x nphi x howmany)
+d = zeros(cfg.nlat_padded, cfg.nphi, howmany)
+d1 = @. 0.3*sin(ϕ') * sin(θ)
+for i=1:howmany
+	@. d[:,:,i] = i*d1
+end
+
+q = SHTns.analys(cfg,d)
+d2 = SHTns.synth(cfg,q)
+d2 ≈ d #true
+```
+
+> [!NOTE]
+> Complex to complex transforms are currently not supported in the batched transform.
+
+

docs/src/api.md

@@ -1,8 +1,51 @@
 # API
 
-```@index
+
+## Configuration
+
+```@docs
+SHTns.SHTnsCfg
+```
+### Transform types (SHTnsType)
+
+```@docs
+SHTns.SHTnsType
+```
+
+```@docs
+SHTns.Gauss
+SHTns.GaussFly
+SHTns.QuickInit
+SHTns.RegDCT
+SHTns.RegFast
+SHTns.RegPoles
+```
+
+### Spherical harmonic normalizations (SHTnsNorm)
+
+```@docs
+SHTns.ForRotations
+SHTns.FourPi
+SHTns.Orthonormal
+SHTns.Schmidt
 ```
 
-```@autodocs
-Modules = [SHTns]
+## Transforms
+
+The transforms from spatial space to spectral space, also called _analysis_, are available through [`SHTns.analys`](@ref) and [`SHTns.analys!`](@ref).
+The transforms from spectral space to spatial space, also called _synthesis_, are available through [`SHTns.synth`](@ref) and [`SHTns.synth!`](@ref).
+
+```@docs
+SHTns.analys
+SHTns.analys!
+SHTns.synth
+SHTns.synth!
+```
+
+## Miscellaneous 
+
+```@docs
+SHTns.LM
+SHTns.gauss_weights
+SHTns.grid
 ```

docs/src/quickstart.md

@@ -0,0 +1,54 @@
+# Quickstart
+
+Load the package
+```@example quickstart
+using SHTns
+```
+
+Define your [`SHTnsCfg`](@ref) configuration for some maximum spherical harmonic degree `lmax`
+```@example quickstart
+lmax = 64
+cfg = SHTnsCfg(lmax)
+```
+
+
+Create some spatial data `d`, using the latitude and longitude from [`SHTns.grid`](@ref)
+```@example quickstart
+function spat_func_test(theta, phi)
+	return 1.0 + 0.01*cos(theta)  + 0.1*(3cos(theta)*cos(theta) - 1.0) +	#// Y00, Y10, Y20
+	(cos(phi) + 0.3*sin(phi)) * sin(theta)	+ #// Y11
+	(cos(2phi) + 0.1*sin(2phi)) * sin(theta)*sin(theta) * (7.0* cos(theta)*cos(theta) - 1.0) * 3/8 	#// Y42
+end
+
+θ,ϕ = SHTns.grid(cfg)
+d = spat_func_test.(θ, ϕ')
+nothing #hide
+```
+
+Transform the data `d` to the spectral coefficients `q` (i.e. perform an analysis)
+```@example quickstart
+q = SHTns.analys(cfg, d)
+nothing #hide
+```
+
+Transform back to spatial data (i.e. perform a synthesis)
+```@example quickstart
+d2 = SHTns.synth(cfg, q)
+d2 ≈ d
+```
+
+These transforms can also be performed in place
+```@example quickstart
+q2 = similar(q)
+SHTns.analys!(cfg, d2, q2)
+q2 ≈ q
+```
+!!! warning 
+	`d2` is overwritten during [`SHTns.analys!`](@ref).
+
+```@example quickstart
+SHTns.synth!(cfg, q2, d2)
+d2 ≈ d
+```
+
+

src/SHTns.jl

@@ -73,6 +73,18 @@ for (type, enumtype) in [(:Orthonormal, :sht_orthonormal), (:FourPi, :sht_fourpi
     end
 end
 
+"""
+    abstract type SHTnsType
+
+SHTnsType is an abstract type for the spherical harmonic transform types. All subtypes contain the following keyword arguments:
+
+- `contiguous_lat::Bool=true`
+- `contiguous_phi::Bool=false`
+- `padding::Bool=false`
+- `gpu::Bool=false`
+- `southpolefirst::Bool=false`
+- `float32::Bool=false`
+"""
 abstract type SHTnsType end
 
 for (type, enumtype) in [(:Gauss, :sht_gauss), (:RegFast, :sht_reg_fast), (:RegDCT, :sht_reg_dct), (:QuickInit, :sht_quick_init), (:RegPoles, :sht_reg_poles), (:GaussFly, :sht_gauss_fly) ] #(:Auto, :sht_auto), 
@@ -117,7 +129,27 @@ function _init_checks(shtype, lmax, mmax, mres, nlat, nphi)
 end
 
 """
-    mutable struct SHTnsCfg{N<:SHTnsNorm, T<:SHTnsType}
+    mutable struct SHTnsCfg{TR<:Union{Real,Complex}, N<:SHTnsNorm, T<:SHTnsType}
+    
+- `cfg::Ptr{shtns_info}`
+- `norm::N`
+- `shtype::T`
+- `robert_form::Bool`
+- `nlm::Int`
+- `lmax::Int`
+- `mmax::Int`
+- `mres::Int`
+- `nlat_2::Int`
+- `nlat::Int`
+- `nphi::Int`
+- `nspat::Int`
+- `li::Vector{Int}`
+- `mi::Vector{Int}`
+- `ct::Vector{Float64}`
+- `st::Vector{Float64}`
+- `nlat_padded::Int`
+- `nlm_cplx::Int`
+- `howmany::Int`
 
 Configuration of spherical harmonic transform.
 """

src/analys.jl

@@ -70,5 +70,24 @@ function analys!(cfg::SHTnsCfg{Complex,T,N}, ur::Tv, utheta::Tv, uphi::Tv, qlm,
     return spat_cplx_to_SHqst(cfg.cfg, ur, utheta, uphi, qlm, slm, tlm)
 end
 
+"""
+    analys!(cfg::SHTnsCfg, v, qlm)
+    analys!(cfg::SHTnsCfg, utheta, uphi, slm, tlm)
+    analys!(cfg::SHTnsCfg, ur, utheta, uphi, qlm, slm, tlm)
+
+In-place transforms of the spatial data into spherical harmonics coefficients for scalar, 2D or 3D fields.
+!!! warning 
+    This function modifies the input arrays `v`, `ur`, `utheta` and `uphi`.
+"""
+function analys! end
+
+"""
+    analys(cfg::SHTnsCfg, v)
+    analys(cfg::SHTnsCfg, utheta, uphi)
+    analys(cfg::SHTnsCfg, ur, utheta, uphi)
+
+Transforms the spatial data into spherical harmonics coefficients `qlm`; `slm` and `tlm`; `qlm`, `slm` and `tlm` for scalar; 2D; 3D fields, respectively.
+"""
+function analys end
 
 export analys, analys!