## FastSphericalHarmonics.jl

Easy-to-use Spherical Harmonics, based on FastTransforms.jl
Author eschnett
Popularity
6 Stars
Updated Last
2 Years Ago
Started In
March 2021

# FastSphericalHarmonics.jl

East-to-use Spherical Harmonics, based on `FastTransforms.jl`

The `FastSphericalHarmonics.jl` package wraps the `FastTransforms.jl` Julia package to calculate spherical harmonics.

`FastTransforms.jl` is a powerful, efficient, and well thought out package. Unfortunately, its user interface is difficult to use for beginners, and its documentation is very technical. `FastSphericalHarmonics.jl` provides functions and documentation that are easier to use. It would be worthwhile to fold `FastSphericalHarmonics.jl` into `FastTransforms.jl` at some point.

## Features

This package provides scalar and spin spherical harmonic transforms for real and complex fields. The normalizations and conventions are chosen to be convenient for real fields, and are different from those usually used for complex spherical harmonics. This is documented by FastTransforms.jl, its underlying C implementation FastTransforms, and in the references listed there (see bottom of the page there).

The documentation lists the implemented functions as well as their Julia signatures. Most functions come in two versions, one that mutates its arguments and one that allocates its output. See also the test cases for more examples.

## Examples

### Example 1: Transform, modify coefficients, and transform back

Make the example reproducible by choosing a particular random number seed:

```julia> using Random

julia> Random.seed!(42);```

Load the package and create a random scalar field on the sphere:

```julia> using FastSphericalHarmonics

julia> lmax = 4;

julia> F = randn(lmax+1, 2lmax+1)
5×9 Matrix{Float64}:
-0.556027   -1.1449    1.08238   -0.886205  -1.05099    0.168341   0.36901      0.681085  -0.145211
-0.444383   -0.468606  0.187028   0.684565   0.502079   0.284259  -0.00761298  -1.33913    0.642896
0.0271553   0.156143  0.518149  -1.59058   -0.216248   0.569829   0.562669    -0.238284   1.81935
-0.299484   -2.64199   1.49138    0.410653  -0.706424  -1.42206    0.106869     1.01936   -0.36726
1.77786     1.00331   0.367563  -0.85635   -3.86593   -0.372401   0.569458     0.701771   0.756569```

Transform to spherical harmonics:

``````julia> C = sph_transform(F)
5×9 Matrix{Float64}:
-0.0971079  -0.480385    0.250188   -0.237836  -0.344185  -1.15475    -0.390183  -0.828615     -0.0736505
0.173019    0.399055   -0.601621    0.235944   1.15235   -0.0698683  -0.222073   0.0436882    -1.03784
-0.35473     0.28504    -0.0294237  -0.492649  -0.931321  -0.478315    0.684783  -0.000600234   0.577805
-0.309621   -0.0516482  -0.99214    -0.318465  -0.243649   0.0434071  -0.452602  -0.338297      0.332204
0.296832   -1.16363     1.65583     0.606001  -0.281404  -0.555203   -0.424356   0.21506       0.123637
``````

Note that the coefficient array `C` contains `(lmax+1) * (2lmax+1) = 45` coefficients, more than the `(lmax+1)^2 = 25` coefficients we expected. There are some higher modes (with `l > lmax`) present as well. This makes the spherical harmonic transform invertible:

```julia> F′ = sph_evaluate(C)
5×9 Matrix{Float64}:
-0.556027   -1.1449    1.08238   -0.886205  -1.05099    0.168341   0.36901      0.681085  -0.145211
-0.444383   -0.468606  0.187028   0.684565   0.502079   0.284259  -0.00761298  -1.33913    0.642896
0.0271553   0.156143  0.518149  -1.59058   -0.216248   0.569829   0.562669    -0.238284   1.81935
-0.299484   -2.64199   1.49138    0.410653  -0.706424  -1.42206    0.106869     1.01936   -0.36726
1.77786     1.00331   0.367563  -0.85635   -3.86593   -0.372401   0.569458     0.701771   0.756569```

Evaluating the modes at the grid points gives the same values back.

Let's create a new coefficient array `C′` that contains only some of the modes, and leaves all other modes set to zero:

```julia> C′ = zeros(size(C));

julia> for l in 0:lmax, m in -l:l
C′[sph_mode(l,m)] = C[sph_mode(l,m)]
end

julia> C′
5×9 Matrix{Float64}:
-0.0971079  -0.480385    0.250188   -0.237836  -0.344185  -1.15475    -0.390183  -0.828615  -0.0736505
0.173019    0.399055   -0.601621    0.235944   1.15235   -0.0698683  -0.222073   0.0        0.0
-0.35473     0.28504    -0.0294237  -0.492649  -0.931321   0.0         0.0        0.0        0.0
-0.309621   -0.0516482  -0.99214     0.0        0.0        0.0         0.0        0.0        0.0
0.296832    0.0         0.0         0.0        0.0        0.0         0.0        0.0        0.0```

We can examine which of the coefficient array elements contains what mode, and which are unused (or rather, used by the supernumerary higher modes):

```julia> lm = similar(C, Any);

julia> for l in 0:lmax, m in -l:l
lm[sph_mode(l,m)] = (l,m)
end

julia> lm
5×9 Matrix{Any}:
(0, 0)     (1, -1)     (1, 1)     (2, -2)     (2, 2)     (3, -3)     (3, 3)     (4, -4)     (4, 4)
(1, 0)     (2, -1)     (2, 1)     (3, -2)     (3, 2)     (4, -3)     (4, 3)  #undef      #undef
(2, 0)     (3, -1)     (3, 1)     (4, -2)     (4, 2)  #undef      #undef     #undef      #undef
(3, 0)     (4, -1)     (4, 1)  #undef      #undef     #undef      #undef     #undef      #undef
(4, 0)  #undef      #undef     #undef      #undef     #undef      #undef     #undef      #undef```

Evaluating these modes gives us scalar field values `F″` that contain only these modes, which are now different from the original values in `F`:

```julia> F″ = sph_evaluate(C′)
5×9 Matrix{Float64}:
-1.09395   -0.814816  -0.0626794   0.497403   0.649439    0.440388   0.0488995  -0.36331   -0.783888
-0.188812  -0.38777    0.535328   -0.220322   0.260065    0.224342  -0.309924   -0.529559   0.657753
0.290636  -0.415552   0.643938   -1.07714    0.0812705   0.622275   0.631893   -0.554523   1.38539
-1.23983   -1.27778    0.806158    0.115463  -1.27163    -1.4813     0.322699    0.908434   0.708826
0.485627   0.410934   0.464254   -0.108253  -1.02271    -1.1944    -0.337374    0.589485   0.794295```

### Example 2: Visualize some modes

Load the library and define the resolution of the images:

```julia> using FastSphericalHarmonics

julia> lmax = 40;

julia> Θ, Φ = sph_points(lmax+1);```

Load Makie to create the images:

```julia> using CairoMakie

julia> ticks = (xticks=MultiplesTicks(4, π, "π"), yticks=MultiplesTicks(4, π, "π"));```

Create images of some modes:

```julia> for l in 4:4, m in 0:l
C = zeros(lmax+1, 2lmax+1)
C[sph_mode(l,m)] = 1
F = sph_evaluate(C)

scene, layout = layoutscene(; resolution=(400, 200))
axis = layout[1, 1] = Axis(scene; xticks=MultiplesTicks(4, π, "π"), yticks=MultiplesTicks(4, π, "π"))
heatmap!(axis, Φ, Θ, F')
scale!(scene, 1, 1)
save("mode\$l\$m.png", scene)
end```

The `l=0`, `m≥0` modes of the spherical harmonics look like this:

### Example 3: Laplacian of Spherical Harmonics

We use both scalar and spin-weighted spherical harmonics to calculate derivatives. Spin-weighted spherical harmonics with spin-weight 0 are the same as scalar spherical harmonics, except with a different normalization.

Some preliminaries:

```julia> using FastSphericalHarmonics

julia> chop(x) = abs2(x) < 100eps(x) ? zero(x) : x;

julia> chop(x::Complex) = Complex(chop(real(x)), chop(imag(x)));```

Choose a function (here `z + 2x` with `z = cos θ` and `x = sin θ cos ϕ`:

```julia> lmax = 4;

julia> Θ, Φ = sph_points(lmax+1);

julia> F = [cos(θ) + sin(θ)*cos(ϕ) for θ in Θ, ϕ in Φ]
5×9 Matrix{Float64}:
1.26007    1.18778     1.00472   …   0.796548   1.00472    1.18778
1.3968     1.20753     0.72827       0.183277   0.72827    1.20753
1.0        0.766044    0.173648     -0.5        0.173648   0.766044
0.221232   0.0319577  -0.447301     -0.992294  -0.447301   0.0319577
-0.64204   -0.714336   -0.897396     -1.10557   -0.897396  -0.714336```

We transform to scalar spherical harmonics, calculate the Laplacian (which is efficient in spectral space), and convert back to point values:

```julia> C = sph_transform(F); chop.(C)
5×9 Matrix{Float64}:
0.0      0.0  2.04665  0.0  0.0  0.0  0.0  0.0  0.0
2.04665  0.0  0.0      0.0  0.0  0.0  0.0  0.0  0.0
0.0      0.0  0.0      0.0  0.0  0.0  0.0  0.0  0.0
0.0      0.0  0.0      0.0  0.0  0.0  0.0  0.0  0.0
0.0      0.0  0.0      0.0  0.0  0.0  0.0  0.0  0.0

julia> ΔC = sph_laplace(C); chop.(ΔC)
5×9 Matrix{Float64}:
0.0      0.0  -4.09331  0.0  0.0  0.0  0.0  0.0  0.0
-4.09331  0.0   0.0      0.0  0.0  0.0  0.0  0.0  0.0
0.0      0.0   0.0      0.0  0.0  0.0  0.0  0.0  0.0
0.0      0.0   0.0      0.0  0.0  0.0  0.0  0.0  0.0
0.0      0.0   0.0      0.0  0.0  0.0  0.0  0.0  0.0

julia> ΔF = sph_evaluate(ΔC)
5×9 Matrix{Float64}:
-2.52015   -2.37555    -2.00943   …  -1.5931    -2.00943   -2.37555
-2.7936    -2.41506    -1.45654      -0.366554  -1.45654   -2.41506
-2.0       -1.53209    -0.347296      1.0       -0.347296  -1.53209
-0.442463  -0.0639154   0.894602      1.98459    0.894602  -0.0639154
1.28408    1.42867     1.79479       2.21113    1.79479    1.42867```

Since we chose `F` to consist of two `l=1` modes, we know its Laplacian: `ΔF = -l(l+1) F`:

```julia> ΔF ≈ -2F
true```

Now let's do the same calculation with spin-weighted spherical harmonics. These are defined for complex functions, so we first convert the real array to a complex array, and then transform to spin-weighted spherical harmonics (of spin-weight `0`).

```julia> F⁰ = Complex.(F);

julia> C⁰ = spinsph_transform(F⁰, 0); chop.(C⁰)
5×9 Matrix{ComplexF64}:
0.0+0.0im  1.4472+0.0im  …  0.0+0.0im  0.0+0.0im  0.0+0.0im
2.04665+0.0im     0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im
0.0+0.0im     0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im
0.0+0.0im     0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im
0.0+0.0im     0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im```

Due to the way in which `FastTransforms.jl` defines spherical harmonics, and because we started with a real function `F`, the imaginary part of `C⁰` is actually zero.

We test that evaluating the spin-weighted spherical harmonics gives us back the original function `F`:

```julia> spinsph_evaluate(C⁰, 0) ≈ F
true```

We then apply the ð ("eth") operator, which is the gradient when applied to a spin-0 function, yielding a spin-1 function.

```julia> ðC¹ = spinsph_eth(C⁰, 0); chop.(ðC¹)
5×9 Matrix{ComplexF64}:
0.0+0.0im  -2.04665+0.0im  …  0.0+0.0im  0.0+0.0im  0.0+0.0im
2.89441+0.0im       0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im
0.0+0.0im       0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im
0.0+0.0im       0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im
0.0+0.0im       0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im```

We can evaluate this spin-1 function on the points on the sphere and thus read off the gradient of `F`:

```julia> ðF¹ = spinsph_evaluate(ðC¹, 1);

julia> ∂θF = real(ðF¹)
5×9 Matrix{Float64}:
0.657164      0.28605    0.0885849  …   1.15716    1.22574     1.02828
1.06331       0.6922     0.494734       1.56331    1.63189     1.43443
1.37051e-16  -0.371114  -0.568579       0.5        0.568579    0.371114
-1.06331      -1.43443   -1.63189       -0.563314  -0.494734   -0.6922
-0.657164     -1.02828   -1.22574       -0.157164  -0.0885849  -0.28605

julia> cscθ∂ϕF = imag(ðF¹)
5×9 Matrix{Float64}:
8.66754e-16  0.642788  0.984808  …  -0.866025  -0.984808  -0.642788
-5.4187e-17   0.642788  0.984808     -0.866025  -0.984808  -0.642788
-2.88227e-16  0.642788  0.984808     -0.866025  -0.984808  -0.642788
1.87818e-16  0.642788  0.984808     -0.866025  -0.984808  -0.642788
-4.87207e-16  0.642788  0.984808     -0.866025  -0.984808  -0.642788```

We thus have `julia> ∇F = (∂θF, sin(θ) * cscθ∂ϕF)`. However, since `sin(θ)` has a coordinate singularity at the pole, it's best not to actually perform this multiplication.

Next we apply the ð̄ ("eth-bar") operator, which is the divergence when applied to a spin-1 function, yielding a spin-0 function again, which we evaluate on the grid points.

```julia> ð̄ðC⁰ = spinsph_ethbar(ðC¹, 1); chop.(ð̄ðC⁰)
5×9 Matrix{ComplexF64}:
0.0+0.0im  -2.89441+0.0im  …  0.0+0.0im  0.0+0.0im  0.0+0.0im
-4.09331+0.0im       0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im
0.0+0.0im       0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im
0.0+0.0im       0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im
0.0+0.0im       0.0+0.0im     0.0+0.0im  0.0+0.0im  0.0+0.0im

julia> ð̄ðF⁰ = spinsph_evaluate(ð̄ðC⁰, 0); chop.(ð̄ðF⁰)
5×9 Matrix{ComplexF64}:
-2.52015+0.0im    -2.37555+0.0im  …   -2.00943+0.0im    -2.37555+0.0im
-2.7936+0.0im    -2.41506+0.0im      -1.45654+0.0im    -2.41506+0.0im
-2.0+0.0im    -1.53209+0.0im     -0.347296+0.0im    -1.53209+0.0im
-0.442463+0.0im  -0.0639154+0.0im      0.894602+0.0im  -0.0639154+0.0im
1.28408+0.0im     1.42867+0.0im       1.79479+0.0im     1.42867+0.0im```

This function `ð̄ðF⁰` is the Laplacian of our original function `F` above. It is complex, but since we started with a real function `F`, `ð̄ðF⁰` has a zero imaginay part (up to round-off):

```julia> maximum(abs.(imag.(ð̄ðF⁰)))
1.1102230246251565e-16```

Of course, both ways of evaluating the Laplacian give the same result:

```julia> real.(ð̄ðF⁰) ≈ ΔF
true```

We can also apply the ð̄ operator first, and then the ð operator:

```julia> ð̄C⁻¹ = spinsph_ethbar(C⁰, 0);

julia> ðð̄C⁰ = spinsph_eth(ð̄C⁻¹, -1);

julia> ðð̄F⁰ = spinsph_evaluate(ðð̄C⁰, 0);

julia> ðð̄F⁰ ≈ ΔF
true```

### Used By Packages

No packages found.