# Julia spherical harmonics different from python

**URL:** <https://discourse.julialang.org/t/julia-spherical-harmonics-different-from-python/88638>\
**Category:** General Usage\
**Tags:** question\
**Created:** [October 12, 2022, 9:50pm UTC](https://discourse.julialang.org/t/julia-spherical-harmonics-different-from-python/88638 "2022-10-12T21:50:13Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![Juliapix](https://avatars.discourse-cdn.com/v4/letter/j/71e660/32.png) [@Juliapix](https://discourse.julialang.org/u/Juliapix)\
**Post date:** [October 12, 2022, 9:50pm UTC](https://discourse.julialang.org/t/julia-spherical-harmonics-different-from-python/88638/1 "2022-10-12T21:50:13Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [October 12, 2022, 10:00pm UTC](https://discourse.julialang.org/t/julia-spherical-harmonics-different-from-python/88638/2 "2022-10-12T22:00:07Z")

</div>

Wasn’t this solved a few days ago by [Julia spherical harmonics different from python - Stack Overflow](https://stackoverflow.com/questions/74026454/julia-spherical-harmonics-different-from-python)?

---

<div class="post-metadata">

**Author:** ![digital\_carver](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/digital_carver/32/33818_2.png) [@digital\_carver](https://discourse.julialang.org/u/digital_carver)\
**Post date:** [October 12, 2022, 11:16pm UTC](https://discourse.julialang.org/t/julia-spherical-harmonics-different-from-python/88638/3 "2022-10-12T23:16:37Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![liuyxpp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/liuyxpp/32/9870_2.png) [@liuyxpp](https://discourse.julialang.org/u/liuyxpp)\
**Post date:** [October 13, 2022, 12:47am UTC](https://discourse.julialang.org/t/julia-spherical-harmonics-different-from-python/88638/4 "2022-10-13T00:47:12Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [October 13, 2022, 12:59am UTC](https://discourse.julialang.org/t/julia-spherical-harmonics-different-from-python/88638/5 "2022-10-13T00:59:44Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [October 13, 2022, 7:40am UTC](https://discourse.julialang.org/t/julia-spherical-harmonics-different-from-python/88638/6 "2022-10-13T07:40:05Z")

</div>

> [@liuyxpp](#):
>
> Put the correctness aside, I think the code in Julia is not as elegant as in Python.

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?

---

<div class="post-metadata">

**Author:** ![liuyxpp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/liuyxpp/32/9870_2.png) [@liuyxpp](https://discourse.julialang.org/u/liuyxpp)\
**Post date:** [October 13, 2022, 8:03am UTC](https://discourse.julialang.org/t/julia-spherical-harmonics-different-from-python/88638/7 "2022-10-13T08:03:17Z")

</div>

OK, let me be more explicit.

In Python, we have

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

```julia
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`:

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

```

Then we can use it like

```julia
GSL.sph_harm(m, l, theta(x,y,z), phi(x,y,z))

```

---

<div class="post-metadata">

**Author:** ![Juliapix](https://avatars.discourse-cdn.com/v4/letter/j/71e660/32.png) [@Juliapix](https://discourse.julialang.org/u/Juliapix)\
**Post date:** [October 14, 2022, 4:49am UTC](https://discourse.julialang.org/t/julia-spherical-harmonics-different-from-python/88638/8 "2022-10-14T04:49:08Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [October 14, 2022, 4:11pm UTC](https://discourse.julialang.org/t/julia-spherical-harmonics-different-from-python/88638/9 "2022-10-14T16:11:46Z")

</div>

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