jishnub / sphericalharmonics.jl

Associated Legendre Polynomials and Spherical Harmonics in Julia

Spherical Harmonics

For a full description of the code, please see:

Associated Legendre Polynomials and Spherical Harmonics Computation for Chemistry Applications (2014). Taweetham Limpanuparb and Josh Milthorpe. arXiv: 1410.1748 [physics.chem-ph]

Quick start

The normalized associated Legendre polynomials for an angle θ for all l in 0 <= l <= lmax and all m in -l <= m <= l may be generated using the signature computePlm(θ; lmax), eg.

julia> P = computePlmcostheta(pi/2, lmax = 1)
3-element SHArray(::Vector{Float64}, (ML(0:1, 0:1),)):
  0.3989422804014327
  4.231083042742082e-17
 -0.4886025119029199

The polynomials are ordered with m increasing faster than l, and the returned array may be indexed using (l,m) Tuples as

julia> P[(0,0)]
0.3989422804014327

julia> P[(1,1)] == P[3]
true

Spherical harmonics for a colatitude θ and azimuth ϕ may be generated using the signature computeYlm(θ, ϕ; lmax), eg.

julia> Y = computeYlm(pi/3, 0, lmax = 1)
4-element SHArray(::Vector{Complex{Float64}}, (ML(0:1, -1:1),)):
  0.2820947917738782 + 0.0im
  0.2992067103010745 - 0.0im
 0.24430125595146002 + 0.0im
 -0.2992067103010745 - 0.0im

The returned array may be indexed using (l,m) Tuples as well, as

julia> Y[(1,-1)]
0.2992067103010745 - 0.0im

julia> Y[(1,-1)] == Y[2]
true

Special angles SphericalHarmonics.NorthPole() and SphericalHarmonics.SouthPole() may be passed as θ to use efficient algorithms.

Increased precision

Arguments of BigInt and BigFloat types would increase the precision of the result.

julia> Y = computeYlm(big(pi)/2, big(0), lmax = 1) # Arbitrary precision
4-element SHArray(::Vector{Complex{BigFloat}}, (ML(0:1, -1:1),)):
    0.2820947917738781434740397257803862929220253146644994284220428608553212342207478 + 0.0im
    0.3454941494713354792652446460318896831393773703262433134867073548945156550201567 - 0.0im
 2.679783085063171668225419916118067917387251852939708540164955895366691604430101e-78 + 0.0im
   -0.3454941494713354792652446460318896831393773703262433134867073548945156550201567 - 0.0im

Semi-positive harmonics

For real functions it might be sufficient to compute only the functions for m >= 0. These may be computed by passing the flag m_range = SphericalHarmonics.ZeroTo.

julia> computeYlm(pi/3, 0, lmax = 1, m_range = SphericalHarmonics.ZeroTo)
3-element SHArray(::Vector{Complex{Float64}}, (ML(0:1, 0:1),)):
  0.2820947917738782 + 0.0im
 0.24430125595146002 + 0.0im
 -0.2992067103010745 - 0.0im

Real harmonics

It's also possible to compute real spherical harmonics by passing the flag SHType = SphericalHarmonics.RealHarmonics(), eg.

julia> Y = computeYlm(pi/3, pi/3, lmax = 1, SHType = SphericalHarmonics.RealHarmonics())
4-element SHArray(::Vector{Float64}, (ML(0:1, -1:1),)):
  0.2820947917738782
 -0.3664518839271899
  0.24430125595146002
 -0.21157109383040865

These are faster to evaluate and require less memory to store.

See also

FastTransforms.jl: The function FastTransforms.sphevaluate is faster at evaluating real spherical harmonics for a single mode.

julia> @btime FastTransforms.sphevaluate($(big(pi)/3), $(big(pi)/3), 100, 100)
  153.142 μs (1336 allocations: 72.64 KiB)
-3.801739606943941485088961175328961189010502022112528054751517912248264631529766e-07

julia> @btime SphericalHarmonics.sphericalharmonic($(big(pi)/3), $(big(pi)/3), 100, 100, SphericalHarmonics.RealHarmonics())
  165.932 μs (1439 allocations: 78.01 KiB)
-3.801739606943941485088961175328961189010502022112528054751517912248264631529107e-07

This difference might reduce in the future.

sphericalharmonics.jl's Issues

Why are my spherical harmonics plots not matching with scipy?

Consider the following attempt to plot $Y_{1,0}$

begin
	θ = collect(0:0.1:2π+0.1)
	ϕ = collect(0:0.05:π+0.05)
	x = zeros(64, 64)
	y = zeros(64, 64)
	z = zeros(64, 64)
	f = zeros(64, 64)
	for i = 1:64
		for j = 1:64
			Y = computeYlm(θ[i],ϕ[j],lmax=1)
			f[i, j] = abs(Y[(1,0)])
			r = f[i, j]
			x[i, j] = r*sin(ϕ[i])*cos(θ[j])
			y[i, j] = r*sin(ϕ[i])*sin(θ[j])
			z[i, j] = r*cos(ϕ[i])
		end
	end
end

It yields the following plot:

compare it with equivalent scipy sph_harm plot,

theta = np.linspace(0,2*np.pi,100)
phi = np.linspace(0,np.pi,100)
x = np.zeros((100,100))
y = np.zeros((100,100))
z = np.zeros((100,100))
f = np.zeros((100,100))
for i in range(100):
    for j in range(100):
        Y10 = sph_harm(0,1,theta[j],phi[i])
        f[i,j] = np.abs(Y10)
        r = f[i,j]
        x[i,j] = r*np.sin(phi[i])*np.cos(theta[j])
        y[i,j] = r*np.sin(phi[i])*np.sin(theta[j])
        z[i,j] = r*np.cos(phi[i])

What I am doing wrong here?

Condon-Shortley Phase Factor

According to the documentation for Plm (https://docs.juliahub.com/SphericalHarmonics/NjDk0/0.1.14/#SphericalHarmonics.computePlmcostheta-Tuple{Any}) and for the spherical harmonics Ylm (https://docs.juliahub.com/SphericalHarmonics/NjDk0/0.1.14/#SphericalHarmonics.computeYlm ) the C-S factor (-1)^m is not included in the definitions. However, that's not what I am observing. For example, the formula for Y21 with the C-S factor is given by

Y21 (theta, phi) = (-)sqrt(15/32pi) sin(2*theta) exp(i phi).

For theta=pi/3, phi=0, it gives a value of -0.3345232717786446.

Now following code produces exactly the same value:
ynm = computeYlm(pi/3,0.0,lmax=2); ynm[(2,1)]
-0.33452327177864466 - 0.0im

I have verified the same behavior for Y(1,1) as well.

Is it possible that the documentation is incorrect? Or am I missing something? [I can think of another possibility, the branch used in defining (1-x^2)^{m/2}) in the definition of Pnm.]

Could you please clarify?

Thanks,
Sanjay Velamparambil

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Support SN3D normalization

Hi @jishnub, thanks for providing this library. I've been using it for Ambisonics 3D audio processing. I saw you adding new normalizations. I was wondering if you could add an option for SN3D normalization with no CS phase term? Like the AmbiX format here. Thanks!

Real matrix of Spherical Harmonics

For multiple points on the sphere what's the best way to get a matrix where the columns are the points and the rows are the spherical harmonic modes for each of the points?

I did this but it does not work if only one point is passed in.

function Y(N::Integer, θ, ϕ)
    n, m = acnorder(N)
    Yres = computeYlm.(θ, ϕ, lmax = N, SHType = SphericalHarmonics.RealHarmonics())
    Ymat = reduce(hcat, Yres)
    # remove the Condon-Shortley phase
    condon = (-1).^m
    condon .* Ymat
end