Opened 2 months ago

Last modified 2 months ago

#31639 new enhancement

Spherical Harmonics: special cases and more examples

Reported by: gh-mjungmath Owned by:
Priority: major Milestone: sage-9.4
Component: misc Keywords:
Cc: egourgoulhon, tscrim, slelievre, mkoeppe Merged in:
Authors: Michael Jung Reviewers:
Report Upstream: N/A Work issues:
Branch: public/spher_harmonics_improved_31639 (Commits, GitHub, GitLab) Commit: e92db1d2f02774aafdb6a407ed24143cc9dc8b0d
Dependencies: #25034, #31637 Stopgaps:

Status badges

Description (last modified by gh-mjungmath)

This ticket is devoted to improve spherical harmonics and make them function properly. That entails treating special cases more efficiently and adding a bunch of doctests for checking correctness.

Change History (12)

comment:1 Changed 2 months ago by gh-mjungmath

  • Dependencies set to #25034

comment:2 Changed 2 months ago by gh-mjungmath

  • Authors set to Michael Jung
  • Branch set to public/spher_harmonics_improved_31639

comment:3 Changed 2 months ago by git

  • Commit set to 1ace9fa32f4dfae500aa2c285f473b33d95bb112

Branch pushed to git repo; I updated commit sha1. New commits:

4fcc574Trac #25034: formula for n=m fixed
a4edb24Trac #25034: abs(m) exceeding n gives zero immediately
1ace9faTrac #31639: special cases added + list of first spherical harmonics added

comment:4 follow-up: Changed 2 months ago by gh-mjungmath

Here is one first approach based on the current branch of #25034.

Any chance to get rid of terms like sqrt(sin^2(theta)) in a nice way?

comment:5 Changed 2 months ago by gh-mjungmath

  • Description modified (diff)

comment:6 Changed 2 months ago by git

  • Commit changed from 1ace9fa32f4dfae500aa2c285f473b33d95bb112 to e92db1d2f02774aafdb6a407ed24143cc9dc8b0d

Branch pushed to git repo; I updated commit sha1. New commits:

d5461d4Trac #31639: import statements given when needed
e92db1dTrac #31639: separate special values from general formula

comment:7 in reply to: ↑ 4 Changed 2 months ago by egourgoulhon

Replying to gh-mjungmath:

Here is one first approach based on the current branch of #25034.

Thanks!

Any chance to get rid of terms like sqrt(sin^2(theta)) in a nice way?

These ones are really annoying in Sage; this is why I introduced the functions simplify_sqrt_real and simplify_abs_trig, which are used in the default simplification chain for calculus on manifolds. They do the job here:

sage: x, y = var('x y')
sage: s = spherical_harmonic(2, 1, x, y); s
1/4*sqrt(6)*sqrt(5)*sqrt(sin(x)^2)*cos(x)*e^(I*y)/sqrt(pi)
sage: assume(0 <= x, x <= pi)                              
sage: s.simplify_full()  # no effect despite the above assumptions!
1/4*sqrt(6)*sqrt(5)*sqrt(sin(x)^2)*cos(x)*e^(I*y)/sqrt(pi)
sage: from sage.manifolds.utilities import simplify_sqrt_real, simplify_abs_trig
sage: simplify_abs_trig(simplify_sqrt_real(s))
1/4*sqrt(5)*sqrt(3)*sqrt(2)*cos(x)*e^(I*y)*sin(x)/sqrt(pi)

comment:8 Changed 2 months ago by egourgoulhon

The following results given in the EXAMPLES section:

Y_1^-1(x,y)=-1/2*sqrt(3/2)*e^(-I*y)*sin(x)/sqrt(pi)
Y_1^1(x,y)=1/4*sqrt(3)*sqrt(2)*sqrt(sin(x)^2)*e^(I*y)/sqrt(pi)
Y_2^-1(x,y)=-1/4*sqrt(15)*sqrt(2)*sqrt(sin(x)^2)*cos(x)*e^(-I*y)/sqrt(pi)
Y_2^1(x,y)=1/4*sqrt(6)*sqrt(5)*sqrt(sin(x)^2)*cos(x)*e^(I*y)/sqrt(pi)

have signs opposite to https://en.wikipedia.org/wiki/Table_of_spherical_harmonics and https://mathworld.wolfram.com/SphericalHarmonic.html. This is certainly due to the Condon-Shortley phase (-1)m. I think we should agree with Wikipedia's convention, since it is standard in mathematical physics. Problably, we should offer the option condon_shortley=True in spherical_harmonic's arguments. For instance, this is done in the package kerrgeodesic_gw:

sage: from kerrgeodesic_gw import spin_weighted_spherical_harmonic
sage: x, y = var('x y')
sage: spin_weighted_spherical_harmonic(0, 1, -1, x, y)  # weight = 0 -> spherical harmonic
0.25*sqrt(6)*e^(-I*y)*sin(x)/sqrt(pi)
sage: spin_weighted_spherical_harmonic(0, 1, 1, x, y)
-1/4*sqrt(6)*e^(I*y)*sin(x)/sqrt(pi)
sage: spin_weighted_spherical_harmonic(0, 1, 1, x, y, condon_shortley=False)
1/4*sqrt(6)*e^(I*y)*sin(x)/sqrt(pi)

By the way, there is not the sqrt(sin(x)^2) issue discussed in comment:7 in the above outputs, maybe because kerrgeodesic_gw uses the algorithm of Golberg (1967) and does not rely explicitely on associated Legendre function; see the function _compute_sw_spherical_harm in line 27 of the file spin_weighted_spherical_harm.py.

PS: kerrgeodesic_gw can be installed in Sage with sage -pip install kerrgeodesic_gw, but if you want to play with it without installing it, you can run this notebook in Binder.

Last edited 2 months ago by egourgoulhon (previous) (diff)

comment:9 Changed 2 months ago by egourgoulhon

Another motivation for using the Condon-Shortley phase: SymPy uses it:

sage: x, y = var('x y')                             
sage: from sympy import Ynm                   
sage: Ynm(1, 1, x, y).expand(func=True)
-sqrt(6)*exp(I*y)*sin(x)/(4*sqrt(pi))

We should certainly agree with SymPy (in addition to Wikipedia, MathWorld, etc.)

comment:10 Changed 2 months ago by gh-mjungmath

  • Cc mkoeppe added

Thank you for the input. I agree that a flag for the Condon-Shortley phase is the most favorable option.

I'll probably use simplify_abs_trig. Why haven't these simplifications made it to the global namespace though? I think it is definitely worth to add.

Let's wait for #25034 before we wrap this up here.

comment:11 Changed 2 months ago by mkoeppe

  • Milestone changed from sage-9.3 to sage-9.4

comment:12 Changed 2 months ago by gh-mjungmath

  • Dependencies changed from #25034 to #25034, #31637
Note: See TracTickets for help on using tickets.