CODEWITHSUNDEEP
pythonjulia
Johny Dow
Johny Dow

Reputation: 61

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.

Upvotes: 3

Views: 404

Answers (1)

Johny Dow
Johny Dow

Reputation: 61

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.

Upvotes: 3

Related Questions