Julia spherical harmonics different from python

General Usage
1

I would like to calculate the Spherical Harmonics with Julia. I have done this with the following code:

using GSL

function radius(x, y, z)
    return sqrt(x^2 + y^2 + z^2)
end

function theta(x, y, z)
    return acos(z / radius(x, y, z))
end

function phi(x, y, z)
    return atan(y, x)
end

function harmonics(l, m, x, y, z)
   return (-1)^(m) * GSL.sf_legendre_sphPlm(l, m, cos(theta(x,y,z)))*ℯ^(im*m*phi(x,y,z)) 
end

harmonics(1, 1, 11.66, -35, -35)
harmonics(1, 1, -35, -35, -35)

The output is the following:

0.07921888327321648 - 0.23779253126608726im
-0.1994711402007164 - 0.19947114020071643im

But doing the same with the following python code:

import scipy.special as spe
import numpy as np

def radius(x, y, z):
    return np.sqrt(x**2 + y**2 + z**2)

def theta(x, y, z):
    return np.arccos(z / radius(x, y, z))

def phi(x, y, z):
    return np.arctan(y / x)

def harmonics(l, m, x, y, z):
    return spe.sph_harm(m, l, phi(x, y, z), theta(x, y, z))

harmonics(1, 1, 11.66, -35, -35)
harmonics(1, 1, -35, -35, -35)

Results in the following output:

(-0.07921888327321645+0.23779253126608718j)
(-0.19947114020071638-0.19947114020071635j)

So the sign of the first result is different. But since only one of the results has a different sign, the cause cannot be in the prefactor (-1)^m. I can’t see through this anymore and can’t explain why the results are different.

2

Wasn’t this solved a few days ago by Julia spherical harmonics different from python - Stack Overflow?

1 Like
3

One thing that wasn’t mentioned in the StackOverflow thread (which I assumed the OP had already figured out) is that the difference in atan behaviour actually affects the second call, not the first. So the prefactor (-1)^(m) does seem to be wrong. After fixing atan, to get parity with the numpy results, you also need to change the prefactor to something like (-1)^(m+1).

4

Put the correctness aside, I think the code in Julia is not as elegant as in Python.

5

One minor thing that will slightly improve accuracy and IMO elegence is to use cis(m*phi(x,y,z)) instead of ℯ^(im*m*phi(x,y,z)

6

What do you have in mind here? The codes look pretty much identical to me. Or do you mean that GSL.sf_legendre_sphPlm could have a simpler interface?

7

OK, let me be more explicit.

In Python, we have

spe.sph_harm(m, l, phi(x, y, z), theta(x, y, z))

spe is the package, sph_harm is a function name, which takes four arguments, m, l, phi, and theta, looking very natural.

While in Julia, we have

GSL.sf_legendre_sphPlm(l, m, cos(theta(x,y,z)))*ℯ^(im*m*phi(x,y,z))

It at least has several downsides:

  • It is lengthy.
  • It is cluttering, with all those cos, ℯ, and im*m things.
  • It is indirect. Users have to write out intermediate steps. (it seems OP has an extra “)” in his code, which clearly demonstrate the harm made by these intermediate steps)
  • The package name GSL though short but hardly reveals it has anything to do with spherical harmonics (especially for someone like me not familiar with the topic).
  • The method name is too long.

I would propose to write a more high-level method provided by GSL like GSL.sph_harm:

sph_harm(m, l, ϕ, θ) = sf_legendre_sphPlm(l, m, cos(θ)*cis(m*ϕ))

Then we can use it like

GSL.sph_harm(m, l, theta(x,y,z), phi(x,y,z))
3 Likes
8

Can you elaborate why this is more accurate? Thanks! =)

9

After more thought, I don’t think it is more accurate. Just faster (because knowing the imaginary part of the input is zero lets you skip a decent amount of work).