eschnett / FastSphericalHarmonics.jl

Easy-to-use Spherical Harmonics, based on FastTransforms.jl
MIT License
19 stars 3 forks source link

Phase error in spin_sph ? #1

Open dubosipsl opened 2 years ago

dubosipsl commented 2 years ago

Hello,

Thanks for writing and sharing FastSphericalHarmonics. I am trying to use it to solve some PDEs on the sphere but have problems with the eth/ethbar operators for scalars. For example :

using Plots, FastSphericalHarmonics
pcolor(args... ; kwargs...) = display(heatmap(args... ; kwargs...));
lmax = 32
Θ, Φ = sph_points(lmax+1);
f = [sin(θ)*cos(ϕ) for θ in Θ, ϕ in Φ]
pcolor(f ; title="f")
ff=Complex.(f)
spinsph_transform!(ff, 0)
ff = spinsph_ethbar(ff,0)
gradff = spinsph_evaluate(ff,-1)
pcolor(abs.(gradff) ; title="abs(gradff)")

The scalar function f has a maximum and minimum at antipodal points on the Equator. index

I expect abs(grad(ff))) to be zero at these extremal points. At first glance abs(grad(ff)) is as expected : it does vanish at two antipodal points on the Equator.

index

However these points are shifted by 90 degrees of longitude relative to the extrema of f.

Am I missing something ? I am running Julia 1.7.1.

Thanks,

Thomas

eschnett commented 2 years ago

I can confirm your results. Yes, this looks wrong.

The spin-weighted spherical harmonics of the FastTransforms C library might be defined differently than usual. This page lists the definition (see near the bottom of the page).

Incidentally, I encountered a similar problem as you with this library yesterday (!). I am now comparing to the ssht library. Any help or insight from your side will be appreciated!

eschnett commented 2 years ago

I confirm that FastTransforms.jl and, by implication, FastSphericalHarmonics.jl use an unexpected sign convention for spin-weighted spherical harmonics. To the best of my knowledge, this function

function change_signs!(flm::AbstractArray{<:Complex}, s::Integer)
    for rowcol in CartesianIndices(flm)
        row, col = Tuple(rowcol)
        m = col == 1 ? 0 : (col % 2 == 0 ? -1 : 1) * (col ÷ 2)
        l = row + max(abs(s), abs(m)) - 1
        sign = (isodd(s) && m < -s) || (m >= -s && isodd(m)) ? -1 : +1
        flm[row, col] = flipsign(flm[row, col], sign)
    end
    return flm
end

converts between this and the "standard" Wikipedia sign convention. With this transformation, the gradient is zero on the x axis.

I will updated/document the code.

eschnett commented 2 years ago

I created a new package ssht.jl based on the ssht library. It uses, of course, a different API. I ran thorough checks and I believe that the sign conventions used there are less surprising.