[sympy] Spherical harmonics (#1510)
raoulb
Improved implementation of spherical harmonics. Based on the additions to the
orthogonal polynomials, the "old" Ylm could handle some symbolic parameters
but only in "expanded", explicit form. This PR adds mainly real symbolic objects
for representing "Ylm". For the moment they are called "Ynm" and a pure addition
but I'd like to replace the old "Ylm" once this is finished.
Please see the docstring to get an idea of the new capabilities.
You can merge this Pull Request by running:
git pull https://github.com/raoulb/sympy spherical_harmonics
Or view, comment on, or merge it at:
https://github.com/sympy/sympy/pull/1510
Commit Summary
- Extended versions of spherical harmonics Ynm
- Add spherical harmonics to doc
- Fix doctests
File Changes
- M doc/src/modules/functions/special.rst (6)
- M sympy/functions/__init__.py (2)
- M sympy/functions/special/spherical_harmonics.py (259)
Patch Links
—
Reply to this email directly or view it on GitHub.
Julien Rioux
In sympy/functions/special/spherical_harmonics.py:
> @@ -8,6 +11,260 @@
>
> _x = Dummy("x")
>
> +class Ynm(Function):
> + r"""
> + Spherical harmonics defined as
> +
> + .. math::
> + Y_n^m(\theta, \varphi) := \sqrt{\frac{(2n+1)(n-m)!}{4\pi(n+m)}}
> + \exp(-i m \varphi)
> + \mathrm{P}_n^m\left(\cos(\theta)\right)
> +
> + Ynm() gives the spherical harmonic function of order `n` and `m`
> + in `\theta` and `\varphi`, :math:`Y_n^m(\theta, \varphi)`.
:math: is the default inline role, so it is sufficient to use backticks only, without writing out :math: in front (as you do for \theta and \varphi).
Julien Rioux
In sympy/functions/special/spherical_harmonics.py:
> + else: > + raise ArgumentIndexError(self, argindex) > + > + def _eval_rewrite_as_polynomial(self, n, m, theta, phi): > + # TODO: Make sure n \in N > + # TODO: Assert |m| <= n ortherwise we should return 0 > + return self.expand(func=True) > + > + def _eval_conjugate(self): > + # TODO: Make sure theta \in R and phi \in R > + n, m, theta, phi = self.args > + return S.NegativeOne**m * self.func(n, -m, theta, phi) > + > + > +def Ynm_c(l, m, theta, phi): > + r"""Conjugate spherical harmonics defined as
Why an extra class when Ynm().conjugate() is enough?
Julien Rioux
In sympy/functions/special/spherical_harmonics.py:
> + See Also > + ======== > + > + Ynm, Ynm_c > + > + References > + ========== > + > + .. [1] http://en.wikipedia.org/wiki/Spherical_harmonics > + .. [2] http://mathworld.wolfram.com/SphericalHarmonic.html > + .. [3] http://functions.wolfram.com/Polynomials/SphericalHarmonicY/ > + """ > + from sympy import simplify > + > + if m > 0: > + zz = (Ynm(n, m, th, ph) + C.NegativeOne()**m * Ynm(n, -m, th, ph)) / sqrt(2)
I'm not sure but it looks like this does not match the definition in the docstring: shouldn't (-1)**m come in front of Ynm(n, m, th, ph)?
raoulb
In sympy/functions/special/spherical_harmonics.py:
> @@ -8,6 +11,260 @@
>
> _x = Dummy("x")
>
> +class Ynm(Function):
> + r"""
> + Spherical harmonics defined as
> +
> + .. math::
> + Y_n^m(\theta, \varphi) := \sqrt{\frac{(2n+1)(n-m)!}{4\pi(n+m)}}
> + \exp(-i m \varphi)
> + \mathrm{P}_n^m\left(\cos(\theta)\right)
> +
> + Ynm() gives the spherical harmonic function of order `n` and `m`
> + in `\theta` and `\varphi`, :math:`Y_n^m(\theta, \varphi)`.
Oh right, thanks, I'm just used to the math too much.
raoulb
In sympy/functions/special/spherical_harmonics.py:
> + else: > + raise ArgumentIndexError(self, argindex) > + > + def _eval_rewrite_as_polynomial(self, n, m, theta, phi): > + # TODO: Make sure n \in N > + # TODO: Assert |m| <= n ortherwise we should return 0 > + return self.expand(func=True) > + > + def _eval_conjugate(self): > + # TODO: Make sure theta \in R and phi \in R > + n, m, theta, phi = self.args > + return S.NegativeOne**m * self.func(n, -m, theta, phi) > + > + > +def Ynm_c(l, m, theta, phi): > + r"""Conjugate spherical harmonics defined as
This is just a wrapper function and not a new class. The reason I introduced it is that the old
code has it, so pure compatibility. Although I could imagine a new full class to be useful as well,
f.e. for representing some formulae in more compact and elegant form. But this is for later.
raoulb
In sympy/functions/special/spherical_harmonics.py:
> + See Also > + ======== > + > + Ynm, Ynm_c > + > + References > + ========== > + > + .. [1] http://en.wikipedia.org/wiki/Spherical_harmonics > + .. [2] http://mathworld.wolfram.com/SphericalHarmonic.html > + .. [3] http://functions.wolfram.com/Polynomials/SphericalHarmonicY/ > + """ > + from sympy import simplify > + > + if m > 0: > + zz = (Ynm(n, m, th, ph) + C.NegativeOne()**m * Ynm(n, -m, th, ph)) / sqrt(2)
Maybe they got out of sync while I was experimenting.
In the mean time I worked out what I think is correct:
{ sqrt(2) * Re(Ynm) for m > 0
Znm = { Ynm for m = 0 # Yn0
{ sqrt(2) * Im(Ynm) for m > 0
next compute the parts by
Re(z) = 1/2 * (z + conj(z))
Im(z) = 1/(2I) * (z - conj(z))
and
conj(Y_{n,m}) = (-1)^m * Y_{n,-m}
If we do this we get a perfectly nice result. However I'm almost sure (except I made an error or did not catch an identity) it differs from what the old code implemented. Hence either it was buggy/wrong or used another definition
for some reason (but then, why?). We should decide how to proceed now.
Julien Rioux
In sympy/functions/special/spherical_harmonics.py:
> + See Also > + ======== > + > + Ynm, Ynm_c > + > + References > + ========== > + > + .. [1] http://en.wikipedia.org/wiki/Spherical_harmonics > + .. [2] http://mathworld.wolfram.com/SphericalHarmonic.html > + .. [3] http://functions.wolfram.com/Polynomials/SphericalHarmonicY/ > + """ > + from sympy import simplify > + > + if m > 0: > + zz = (Ynm(n, m, th, ph) + C.NegativeOne()**m * Ynm(n, -m, th, ph)) / sqrt(2)
See the wikipedia article on spherical harmonics, especially the section about the Condon–Shortley phase. That explains the different definition.
Julien Rioux
SymPy Bot Summary: 
There were test failures.
@raoulb: Please fix the test failures.
Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYi98kDA
Interpreter: /usr/bin/python (2.7.2-final-0)
Architecture: Linux (64-bit)
Cache: yes
Test command: setup.py test
master hash: 161c963
branch hash: 54089c4
Automatic review by SymPy Bot.
Julien Rioux
SymPy Bot Summary: There were test failures.
@raoulb: Please fix the test failures.
Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sY7L8kDA
Interpreter: /usr/bin/python (2.7.0-final-0)
Architecture: Linux (32-bit)
Automatic review by SymPy Bot.
—
Reply to this email directly or view it on GitHub.
raoulb
In sympy/functions/special/spherical_harmonics.py:
> + See Also > + ======== > + > + Ynm, Ynm_c > + > + References > + ========== > + > + .. [1] http://en.wikipedia.org/wiki/Spherical_harmonics > + .. [2] http://mathworld.wolfram.com/SphericalHarmonic.html > + .. [3] http://functions.wolfram.com/Polynomials/SphericalHarmonicY/ > + """ > + from sympy import simplify > + > + if m > 0: > + zz = (Ynm(n, m, th, ph) + C.NegativeOne()**m * Ynm(n, -m, th, ph)) / sqrt(2)
Yes I see. We include it into the definition of the legendre functions where P_n^{-m} = (-1)^m * ... * P_n^m.
Julien Rioux
SymPy Bot Summary: 
All tests have passed.
Test results html report: http://reviews.sympy.org/report/agZzeW1weTNyDAsSBFRhc2sYk-IjDA
Interpreter: /usr/bin/python (2.7.2-final-0)
Architecture: Linux (64-bit)
Cache: yes
Test command: setup.py test
master hash: 161c963
branch hash: 31fbd60
Automatic review by SymPy Bot.
—
Reply to this email directly or view it on GitHub.
Ondřej Čertík
@raoulb, many many thanks for this, it was long time needed. Before we merge this, let me play with this, I need to check that the we use the "right" conventions for the phase. I spent some time deriving everything for both complex and real spherical harmonics:
http://theoretical-physics.net/dev/src/math/spherical-harmonics.html#spherical-harmonics
but I forgot at the moment which exact convention Mathematica uses and so on. I want to make this right and well documented. I'll try to get to this tonight or tomorrow.
raoulb
Thanks :-)
I'd really appreciate if you double-check my work here.
Hopefully I put all necessary formulae into the docstrings
so that it is absolutely clear what the code is supposed to do.
Please note that we put the (-1)^m factor into the Legendre polys.
I'll read your chapter tomorrow.
Before we can merge this, we need to work out some failing
corner cases (see comments above) for special values of m and n.
And some more tests have to be written too.
raoulb
Oh, by the way, once this work is finished we can remove the old code
and rename everything to use l and m.
Ondřej Čertík
I sent a wrong chapter. I have to unify my notes. Here is the relevant chapter:
http://theoretical-physics.net/dev/src/math/operators.html#real-spherical-harmonics
and in particular, read this: http://theoretical-physics.net/dev/src/math/operators.html#tables
and at the end, there is a script using SymPy, that calculates the Legendre polys and complex + real spherical harmonics, and prints examples. I think sometimes there is some confusion for m < 0. Yes, The Condon & Shortley (-1)^m factor should go into Legendre polys, I actually did the same in my notes. Finally, we should document in the docstrings whether Mathematica is using the same convention. I think it is.
I'm busy tonight, but let me get to this tomorrow.
raoulb
> relevant chapter:
>
> http://theoretical-physics.net/dev/src/math/operators.html#real-spherical-harmonics
>
> and in particular, read this:
> http://theoretical-physics.net/dev/src/math/operators.html#tables
Please see attached file.
For the real ones we get incomplete simplification.
Consider the following example:
In [110]: Znm(1, 1, theta, phi)
Out[110]: sqrt(2)*(Ynm(1, 1, theta, phi) + exp(-2*I*phi)*Ynm(1, 1, theta, phi))/2
In [111]: expand_func(_)
Out[111]: sqrt(2)*(-sqrt(6)*sqrt(-cos(theta)**2 + 1)*exp(I*phi)/(4*sqrt(pi)) - sqrt(6)*sqrt(-cos(theta)**2 + 1)*exp(-I*phi)/(4*sqrt(pi)))/2
In [112]: simplify(_)
Out[112]: sqrt(3)*(-exp(2*I*phi) - 1)*exp(-I*phi)*Abs(sin(theta))/(4*sqrt(pi))
In [113]: _ / cos(phi).rewrite(exp)
Out[113]: sqrt(3)*(-exp(2*I*phi) - 1)*exp(-I*phi)*Abs(sin(theta))/(4*sqrt(pi)*(exp(I*phi)/2 + exp(-I*phi)/2))
In [114]: simplify(_)
Out[114]: -sqrt(3)*Abs(sin(theta))/(2*sqrt(pi))
We see that the exponential could be combined into
the required factor cos(phi). Note that the Abs(sin(theta))
does no harm given the range of \theta \in [0, \pi].
The negative one works too:
In [116]: Znm(1, -1, theta, phi)
Out[116]: -sqrt(2)*I*(Ynm(1, 1, theta, phi) - exp(-2*I*phi)*Ynm(1, 1, theta, phi))/2
In [117]: expand_func(_)
Out[117]: -sqrt(2)*I*(-sqrt(6)*sqrt(-cos(theta)**2 + 1)*exp(I*phi)/(4*sqrt(pi)) + sqrt(6)*sqrt(-cos(theta)**2 + 1)*exp(-I*phi)/(4*sqrt(pi)))/2
In [118]: simplify(_)
Out[118]: sqrt(3)*I*(exp(2*I*phi) - 1)*exp(-I*phi)*Abs(sin(theta))/(4*sqrt(pi))
In [119]: _ / sin(phi).rewrite(exp)
Out[119]: -sqrt(3)*(exp(2*I*phi) - 1)*exp(-I*phi)*Abs(sin(theta))/(2*sqrt(pi)*(exp(I*phi) - exp(-I*phi)))
In [120]: simplify(_)
Out[120]: -sqrt(3)*Abs(sin(theta))/(2*sqrt(pi))
I conclude that Znm(1, {1, 0, -1}, ...) agree with your list.
> and at the end, there is a script using SymPy, that calculates the
> Legendre polys and complex + real spherical harmonics, and prints
> examples. I think sometimes there is some confusion for m < 0.
> Yes, The Condon & Shortley (-1)^m factor should go into Legendre polys, I
> actually did the same in my notes.
> Finally, we should document in the docstrings whether Mathematica is
> using the same convention. I think it is.
> I'm busy tonight, but let me get to this tomorrow.
raoulb
```
from sympy import *
theta = Symbol("theta", real=True)
phi = Symbol("phi", real=True)
N = 3
for n in range(N+1):
for m in range(-n, n+1):
Y = Ynm(n, m, theta, phi)
Y = Y.expand(func=True)
Y = simplify(Y)
print((n, m))
pprint(Y)
print(80*"-")
print(80*"=")
print(5*"\n")
for n in range(N+1):
for m in range(-n, n+1):
Z = Znm(n, m, theta, phi)
Z = Z.expand(func=True)
Z = simplify(Z)
print((n, m))
pprint(Z)
print(80*"-")
print(80*"=")
```
raoulb
> is using the same convention. I think it is.
http://docs.scipy.org/doc/scipy/reference/generated/scipy.special.sph_harm.html
This may be used together with sympy more often than MMA.
Ondřej Čertík
Excellent. Should we write some regular tests as well for this? What about the old Ylm functions, should we remove them?
raoulb
Of course we need some more tests for the new functionalities.
I' like to remove old code for at least two reasons. First I is more or less redundant,
new code should be able to do everything the old could. And second it adds confusion
about which one to use.
But it uses some other convention at some places, for f.e. Zlm.
We should check what the differences are. Maybe deprecating
instead of removing right now is the better solution.
raoulb
A git grep shows that Ylm is not used in the physics code. Is this really true?
What about the hydrogen atom? I thought we have an example on this?
Julien Rioux
The old functions return the spherical harmonics expressed in terms of cosines and sines and already in the simplest form. There should be an easy way to do that with the new code, something like what .expand(func=True) does in your docstring, but we shouldn't need to call simplify on the result. Also, .expand(func=True) wasn't very intuitive; I'm not sure what sympy's convention is for this.
Julien Rioux
SymPy Bot Summary: Passed after merging raoulb/spherical_harmonics (6ab4400) into master (c49dca1).
Python 2.7.2-final-0: pass
Python 3.2.1-final-0: pass
Sphinx 1.1.3: pass
Docs build command: make clean && make html-errors && make clean && make latex && cd _build/latex && xelatex -interaction=nonstopmode sympy-*.tex
Julien Rioux
In sympy/core/tests/test_args.py:
> @@ -1106,6 +1106,14 @@ def test_sympy__functions__special__polynomials__assoc_laguerre(): > from sympy.functions.special.polynomials import assoc_laguerre > assert _test_args(assoc_laguerre(x, 0, y)) > > +def test_sympy__functions__special__spherical_harmonics__Ynm(): > + from sympy.functions.special.spherical_harmonics import Ynm > + assert _test_args(Ynm(1,1,x,y))
Add a space after each comma.
Julien Rioux
In sympy/core/tests/test_args.py:
> @@ -1106,6 +1106,14 @@ def test_sympy__functions__special__polynomials__assoc_laguerre(): > from sympy.functions.special.polynomials import assoc_laguerre > assert _test_args(assoc_laguerre(x, 0, y)) > > +def test_sympy__functions__special__spherical_harmonics__Ynm(): > + from sympy.functions.special.spherical_harmonics import Ynm > + assert _test_args(Ynm(1,1,x,y)) > + > +def test_sympy__functions__special__spherical_harmonics__Znm(): > + from sympy.functions.special.spherical_harmonics import Znm > + assert _test_args(Znm(1,1,x,y))
Julien Rioux
In sympy/functions/special/spherical_harmonics.py:
> + > + # Handle negative index m and arguments theta, phi > + if m.could_extract_minus_sign(): > + m = -m > + return S.NegativeOne**m * C.exp(-2*I*m*phi) * Ynm(n, m, theta, phi) > + if theta.could_extract_minus_sign(): > + theta = -theta > + return Ynm(n, m, theta, phi) > + if phi.could_extract_minus_sign(): > + phi = -phi > + return C.exp(-2*I*m*phi) * Ynm(n, m, theta, phi) > + > + # Move this into the expand_func routine? > + #if n.is_number and m.is_number: > + # return (sqrt((2*n+1)/(4*pi) * C.factorial(n-m)/C.factorial(n+m)) * > + # C.exp(I*m*phi) * assoc_legendre(n, m, C.cos(theta)))
This should be removed.
Julien Rioux
In sympy/functions/special/spherical_harmonics.py:
> + return Ynm(n, m, theta, phi) > + if phi.could_extract_minus_sign(): > + phi = -phi > + return C.exp(-2*I*m*phi) * Ynm(n, m, theta, phi) > + > + # Move this into the expand_func routine? > + #if n.is_number and m.is_number: > + # return (sqrt((2*n+1)/(4*pi) * C.factorial(n-m)/C.factorial(n+m)) * > + # C.exp(I*m*phi) * assoc_legendre(n, m, C.cos(theta))) > + > + # TODO Add more simplififcation here > + > + > + def _eval_expand_func(self, **hints): > + n, m, theta, phi = self.args > + return (sqrt((2*n+1)/(4*pi) * C.factorial(n-m)/C.factorial(n+m)) *
Add spaces around +.
Julien Rioux
In sympy/functions/special/spherical_harmonics.py:
> + n, m, theta, phi = self.args > + return (sqrt((2*n+1)/(4*pi) * C.factorial(n-m)/C.factorial(n+m)) * > + C.exp(I*m*phi) * assoc_legendre(n, m, C.cos(theta))) > + > + def fdiff(self, argindex=4): > + if argindex == 1: > + # Diff wrt n > + raise ArgumentIndexError(self, argindex) > + elif argindex == 2: > + # Diff wrt m > + raise ArgumentIndexError(self, argindex) > + elif argindex == 3: > + # Diff wrt theta > + n, m, theta, phi = self.args > + return (m * C.cot(theta) * Ynm(n, m, theta, phi) + > + sqrt((n-m)*(n+m+1)) * C.exp(-I*phi) * Ynm(n, m+1, theta, phi))
Add spaces around + and -.
Julien Rioux
In sympy/functions/special/spherical_harmonics.py:
> + raise ArgumentIndexError(self, argindex) > + > + def _eval_rewrite_as_polynomial(self, n, m, theta, phi): > + # TODO: Make sure n \in N > + # TODO: Assert |m| <= n ortherwise we should return 0 > + return self.expand(func=True) > + > + def _eval_conjugate(self): > + # TODO: Make sure theta \in R and phi \in R > + n, m, theta, phi = self.args > + return S.NegativeOne**m * self.func(n, -m, theta, phi) > + > + def as_real_imag(self, deep=True, **hints): > + # TODO: Handle deep and hints > + n, m, theta, phi = self.args > + re = (sqrt((2*n+1)/(4*pi) * C.factorial(n-m)/C.factorial(n+m)) *
Julien Rioux
In sympy/functions/special/spherical_harmonics.py:
> + def _eval_rewrite_as_polynomial(self, n, m, theta, phi): > + # TODO: Make sure n \in N > + # TODO: Assert |m| <= n ortherwise we should return 0 > + return self.expand(func=True) > + > + def _eval_conjugate(self): > + # TODO: Make sure theta \in R and phi \in R > + n, m, theta, phi = self.args > + return S.NegativeOne**m * self.func(n, -m, theta, phi) > + > + def as_real_imag(self, deep=True, **hints): > + # TODO: Handle deep and hints > + n, m, theta, phi = self.args > + re = (sqrt((2*n+1)/(4*pi) * C.factorial(n-m)/C.factorial(n+m)) * > + C.cos(m*phi) * assoc_legendre(n, m, C.cos(theta))) > + im = (sqrt((2*n+1)/(4*pi) * C.factorial(n-m)/C.factorial(n+m)) *
Julien Rioux
In sympy/functions/special/spherical_harmonics.py:
> + > + .. [1] http://en.wikipedia.org/wiki/Spherical_harmonics > + .. [2] http://mathworld.wolfram.com/SphericalHarmonic.html > + .. [3] http://functions.wolfram.com/Polynomials/SphericalHarmonicY/ > + """ > + > + nargs = 4 > + > + @classmethod > + def eval(cls, n, m, theta, phi): > + n, m, th, ph = [sympify(x) for x in (n, m, theta, phi)] > + > + if m.is_positive: > + zz = (Ynm(n, m, th, ph) + Ynm_c(n, m, th, ph)) / sqrt(2) > + #zz = zz.expand(complex=True) > + #zz = simplify(zz)
Should this be removed?
Julien Rioux
In sympy/functions/special/spherical_harmonics.py:
> + > + @classmethod > + def eval(cls, n, m, theta, phi): > + n, m, th, ph = [sympify(x) for x in (n, m, theta, phi)] > + > + if m.is_positive: > + zz = (Ynm(n, m, th, ph) + Ynm_c(n, m, th, ph)) / sqrt(2) > + #zz = zz.expand(complex=True) > + #zz = simplify(zz) > + return zz > + elif m.is_zero: > + return Ynm(n, m, th, ph) > + elif m.is_negative: > + zz = (Ynm(n, m, th, ph) - Ynm_c(n, m, th, ph)) / (sqrt(2)*I) > + #zz = zz.expand(complex=True) > + #zz = simplify(zz)
and this?
Julien Rioux
raoulb
The old functions return the spherical harmonics expressed in terms of cosines and sines and already in the simplest form.
Yes, but automatically "expanding" out Y to the explicit form in terms of P(cos) is a bad idea, see below.
There should be an easy way to do that with the new code, something like what .expand(func=True) does in your docstring, but we shouldn't need to call simplify on the result.
I agree that the version coming from expand is not always the simplest form. An easy solution would be to include
a call to "simplify" into this function. However, "simplify" can be pretty expensive, so I'm not sure if we want
to do this always. Sometimes it might not matter to the user which form he gets because he does further
computations on it, essentially this result is hidden within the computation chain. And at any time the user
can manually call simplify if he wishes to do so.
Also, .expand(func=True) wasn't very intuitive;
AFAIK, This is the default way to code such things.
We do not want automatic expansion, because there is no way to turn it off.
And this would prevent writing series expansions like the Laplace expansion
without getting a huge mess.
BTW: I even think that we should revert this for the orthogonal polys,
because ATM we can not write down a formal Chebyshev series
expansion, because all T_n expand out to polys for specific, given n.
raoulb
In sympy/functions/special/spherical_harmonics.py:
> + > + .. [1] http://en.wikipedia.org/wiki/Spherical_harmonics > + .. [2] http://mathworld.wolfram.com/SphericalHarmonic.html > + .. [3] http://functions.wolfram.com/Polynomials/SphericalHarmonicY/
> + """ > + > + nargs = 4 > + > + @classmethod > + def eval(cls, n, m, theta, phi): > + n, m, th, ph = [sympify(x) for x in (n, m, theta, phi)] > + > + if m.is_positive: > + zz = (Ynm(n, m, th, ph) + Ynm_c(n, m, th, ph)) / sqrt(2) > + #zz = zz.expand(complex=True) > + #zz = simplify(zz)
Depends on what we decide about simplification of expanded forms.
raoulb
In sympy/functions/special/spherical_harmonics.py:
> + > + @classmethod > + def eval(cls, n, m, theta, phi): > + n, m, th, ph = [sympify(x) for x in (n, m, theta, phi)] > + > + if m.is_positive: > + zz = (Ynm(n, m, th, ph) + Ynm_c(n, m, th, ph)) / sqrt(2) > + #zz = zz.expand(complex=True) > + #zz = simplify(zz) > + return zz > + elif m.is_zero: > + return Ynm(n, m, th, ph) > + elif m.is_negative: > + zz = (Ynm(n, m, th, ph) - Ynm_c(n, m, th, ph)) / (sqrt(2)*I) > + #zz = zz.expand(complex=True) > + #zz = simplify(zz)
Same
Julien Rioux
SymPy Bot Summary: Failed after merging raoulb/spherical_harmonics (c8af9e4) into master (423f733).
@raoulb: Please fix the test failures.
Python 2.7.2-final-0: fail
Python 3.2.1-final-0: pass
Sphinx 1.1.3: fail
Docs build command: make clean && make html-errors && make clean && make latex && cd _build/latex && xelatex sympy-*.tex
—
Reply to this email directly or view it on GitHub.
Julien Rioux
I since realised that .expand(func=True) is used across sympy; that it is not greatly documented is outside this PR. However I still think there should be a way to expand into cos. Would .rewrite(cos) do it?
Could not fetch the branch raoulb/spherical_harmonics.raoulb
Would .rewrite(cos) do it?
Added .rewrite(sin) and .rewrite(cos). Boths methods do the same simplification. I added both to make it easier
to remember. We can not do .rewrite(trigo) or similar, can we?
Please test and comment if this is what you expect.
Aaron Meurer
You can do .rewrite("whatever"), where you just use a string. But try to keep the API sane.
Julien Rioux
At least in terms of replacements for the current Ylm and Zlm it is not quite there yet:
>>> from sympy import *
>>> x, y = symbols("x y", real=True)
>>> Ynm(1, -1, x, y)
-exp(-2*I*y)*Ynm(1, 1, x, y)
>>> Ynm(1, -1, x, y).rewrite(sin)
sqrt(6)*(-I*sin(2*y) + cos(2*y))*exp(I*y)*sin(x)/(4*sqrt(pi))
>>> Ylm(1, -1, x, y)
sqrt(6)*exp(-I*y)*sin(x)/(4*sqrt(pi))
>>> Znm(1, -1, x, y)
-sqrt(2)*I*(Ynm(1, 1, x, y) - exp(-2*I*y)*Ynm(1, 1, x, y))/2
>>> Znm(1, -1, x, y).rewrite(sin)
-sqrt(2)*I*(sqrt(6)*(-I*sin(2*y) + cos(2*y))*exp(I*y)*sin(x)/(4*sqrt(pi)) - sqrt(6)*exp(I*y)*sin(x)/(4*sqrt(pi)))/2
>>> Zlm(1, -1, x, y)
sqrt(3)*sin(x)*sin(y)/(2*sqrt(pi))
but maybe that's just asking too much from simplify.
Julien Rioux
SymPy Bot Summary: Passed after merging raoulb/spherical_harmonics (9ff19ac) into master (250d764).
Python 2.7.2-final-0: pass
Python 3.2.1-final-0: pass
Sphinx 1.1.3: pass
Docs build command: make clean && make html-errors && make clean && make latex && cd _build/latex && xelatex sympy-*.tex
—
Reply to this email directly or view it on GitHub.
Julien Rioux
SymPy Bot Summary: 
There were merge conflicts (could not merge raoulb/spherical_harmonics (9ff19ac) into master (57e94e4)); could not test the branch.
@raoulb: Please rebase or merge your branch with master. See the report for a list of the merge conflicts.
raoulb
It is possible to get back the more compact form, see for example this demonstration:
In [29]: Ynm(1, -1, theta, phi).rewrite(sin)
Out[29]: sqrt(6)*(-I*sin(2*phi) + cos(2*phi))*exp(I*phi)*sin(theta)/(4*sqrt(pi))
In [30]: expand_complex(_)
Out[30]: sqrt(6)*sin(theta)*re((I*sin(phi) + cos(phi))*(-I*sin(2*phi) + cos(2*phi)))/(4*sqrt(pi)) + sqrt(6)*I*sin(theta)*im((I*sin(phi) + cos(phi))*(-I*sin(2*phi) + cos(2*phi)))/(4*sqrt(pi))
In [31]: expand_trig(_)
Out[31]: sqrt(6)*sin(theta)*re((I*sin(phi) + cos(phi))*(-2*I*sin(phi)*cos(phi) + 2*cos(phi)**2 - 1))/(4*sqrt(pi)) + sqrt(6)*I*sin(theta)*im((I*sin(phi) + cos(phi))*(-2*I*sin(phi)*cos(phi) + 2*cos(phi)**2 - 1))/(4*sqrt(pi))
In [32]: simplify(_)
Out[32]: sqrt(6)*(-I*sin(phi) + cos(phi))*sin(theta)/(4*sqrt(pi))
In [33]: Ylm(1, -1, theta, phi).rewrite(sin)
Out[33]: sqrt(6)*(-I*sin(phi) + cos(phi))*sin(theta)/(4*sqrt(pi))
The question is how much of these calls we want to hardcode into the class.
I'd rather leave this up to the user. There are several reasons for this:
- Simplification is mostly a heuristic process. I'm not sure in which order the calls yield the best result, and if the do so always.
- The user knows what he wants. Undoing built-in processes is non-trivial if not impossible while manually calling any of these functions is easy.
This has also technical advantages:
- We can avoid potentially expensive simplification calls as long as we don't really need or want them.
- We can easily apply more powerful algorithms that might get implemented in the future without having to change this code.
Maybe we should put a hint into the docstring, saying that it might be necessary to call this and that function
in order to obtain maximally simplified results. What do you think?
Julien Rioux
Passed after merging raoulb/spherical_harmonics (Ondřej Čertík
What about the old Ylm? Should those be removed? Should Ynm be renamed to Ylm? Your code seems ok to me, but we should not have two kinds of spherical harmonics. Is there any reason to keep the old ones around?
Julien Rioux
The new functions are still not useful as drop-in replacements of the old ones. Quoting:
I'm not sure in which order the calls yield the best result, and if the do so always.
[...]
Maybe we should put a hint into the docstring, saying that it might be necessary
to call this and that function in order to obtain maximally simplified results.
I think we should figure out exactly which simplifying functions are needed, and provide a single simple interface to apply them (be it in .rewrite or whatever). As stated, this is not trivial, that's why I wouldn't want to require the user to jump through hoops with the various simplifying functions to figure it out. Let's do it once and have this easily accessible already.
raoulb
What about the old Ylm? Should those be removed?
Yes, at some point.
Should Ynm be renamed to Ylm?
Yes. However this breaks compatibility.
Your code seems ok to me, but we should not have two kinds of spherical harmonics.
Sure. But maybe we need a transition period?
Is there any reason to keep the old ones around?
Not really.
The new functions are still not useful as drop-in replacements of the old ones.
They never will be. They can not because we do no evaluation by default.
I think we should figure out exactly which simplifying functions are needed
If you really want to do this, please go ahead.
You can take my example above as starting point.
I still think it is not the responsibility of Ynm to provide fully simplified results.
Please see my comments above.
Ondřej Čertík
Here is what I propose:
1) Implement Ylm using Ynm + proper function calls in order to reproduce the old results. Is that possible? Or are there some cases that actually cannot even be calculated using the new Ynm? Implement a new method, something like Ynm.simple_formula() or something, that would do this. Then implement Ylm as a single call to Ynm.simple_formula(). Does this make sense?
2) Document in the docstring of Ynm, that the user can call either Ylm which returns a formula, or Ynm, which is a class and can be used to represent the spherical harmonic symbolically.
3) Merge this PR
4) Later, we can make the incompatibility switch and remove Ylm, and rename Ynm to Ylm.
Ondřej Čertík
Alternatively, why cannot you reuse the code from Ylm and simply put it into Ynm, as a method, something like simple_formula() method or something.
My whole point is that I don't want to have two incompatible implementations of spherical harmonics. Let's put all the functionality into one Ynm class and then we can merge this PR.
Brian E. Granger
>
>
--
Brian E. Granger
Cal Poly State University, San Luis Obispo
bgra...@calpoly.edu and elli...@gmail.com
raoulb
We can also remove the old one right now. I would not want to mess up the new class with the old code.
(If you want to do it yourself, then please go on.)