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.