# Julia spherical harmonics different from python

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.

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

1 Like

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)`.

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

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)`

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?

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

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

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).