3

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.

3
  • 1
    Some special functions have multiple different conventions. Have you checked the docs to make sure they are using the same conventions? Commented Oct 11, 2022 at 12:20
  • 1
    @OscarSmith I checked this doc for GSL gnu.org/software/gsl/doc/html/specfunc.html and this docs.scipy.org/doc/scipy/reference/generated/… for scipy. The definition seems to match. Commented Oct 11, 2022 at 12:41
  • 2
    The other potential issue is whether cos and exp are taking the branch cuts that you expect. Commented Oct 11, 2022 at 12:47

1 Answer 1

3

The comment by @Oscar Smith got me started on the solution. Julia uses a different convention for the angles returned by atan, provided two arguments are passed [Julia, Numpy]. If we use atan(y / x) instead of atan(y, x) in Julia we get the same result.

Sign up to request clarification or add additional context in comments.

2 Comments

For Julia you can link directly to the function's doc: Link. As the docs say, the 2 argument form corresponds to atan2 or np.arctan2 (so, atan(-35, -35) in Julia and np.arctan2(-35, -35) are equal).
For correctness you should switch both to a two-argument function. The one-argument form will give incorrect results for odd m over half the range.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.