Coverage for /usr/lib/python3/dist-packages/sympy/functions/special/spherical_harmonics.py: 34%
79 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1from sympy.core.expr import Expr
2from sympy.core.function import Function, ArgumentIndexError
3from sympy.core.numbers import I, pi
4from sympy.core.singleton import S
5from sympy.core.symbol import Dummy
6from sympy.functions import assoc_legendre
7from sympy.functions.combinatorial.factorials import factorial
8from sympy.functions.elementary.complexes import Abs, conjugate
9from sympy.functions.elementary.exponential import exp
10from sympy.functions.elementary.miscellaneous import sqrt
11from sympy.functions.elementary.trigonometric import sin, cos, cot
13_x = Dummy("x")
15class Ynm(Function):
16 r"""
17 Spherical harmonics defined as
19 .. math::
20 Y_n^m(\theta, \varphi) := \sqrt{\frac{(2n+1)(n-m)!}{4\pi(n+m)!}}
21 \exp(i m \varphi)
22 \mathrm{P}_n^m\left(\cos(\theta)\right)
24 Explanation
25 ===========
27 ``Ynm()`` gives the spherical harmonic function of order $n$ and $m$
28 in $\theta$ and $\varphi$, $Y_n^m(\theta, \varphi)$. The four
29 parameters are as follows: $n \geq 0$ an integer and $m$ an integer
30 such that $-n \leq m \leq n$ holds. The two angles are real-valued
31 with $\theta \in [0, \pi]$ and $\varphi \in [0, 2\pi]$.
33 Examples
34 ========
36 >>> from sympy import Ynm, Symbol, simplify
37 >>> from sympy.abc import n,m
38 >>> theta = Symbol("theta")
39 >>> phi = Symbol("phi")
41 >>> Ynm(n, m, theta, phi)
42 Ynm(n, m, theta, phi)
44 Several symmetries are known, for the order:
46 >>> Ynm(n, -m, theta, phi)
47 (-1)**m*exp(-2*I*m*phi)*Ynm(n, m, theta, phi)
49 As well as for the angles:
51 >>> Ynm(n, m, -theta, phi)
52 Ynm(n, m, theta, phi)
54 >>> Ynm(n, m, theta, -phi)
55 exp(-2*I*m*phi)*Ynm(n, m, theta, phi)
57 For specific integers $n$ and $m$ we can evaluate the harmonics
58 to more useful expressions:
60 >>> simplify(Ynm(0, 0, theta, phi).expand(func=True))
61 1/(2*sqrt(pi))
63 >>> simplify(Ynm(1, -1, theta, phi).expand(func=True))
64 sqrt(6)*exp(-I*phi)*sin(theta)/(4*sqrt(pi))
66 >>> simplify(Ynm(1, 0, theta, phi).expand(func=True))
67 sqrt(3)*cos(theta)/(2*sqrt(pi))
69 >>> simplify(Ynm(1, 1, theta, phi).expand(func=True))
70 -sqrt(6)*exp(I*phi)*sin(theta)/(4*sqrt(pi))
72 >>> simplify(Ynm(2, -2, theta, phi).expand(func=True))
73 sqrt(30)*exp(-2*I*phi)*sin(theta)**2/(8*sqrt(pi))
75 >>> simplify(Ynm(2, -1, theta, phi).expand(func=True))
76 sqrt(30)*exp(-I*phi)*sin(2*theta)/(8*sqrt(pi))
78 >>> simplify(Ynm(2, 0, theta, phi).expand(func=True))
79 sqrt(5)*(3*cos(theta)**2 - 1)/(4*sqrt(pi))
81 >>> simplify(Ynm(2, 1, theta, phi).expand(func=True))
82 -sqrt(30)*exp(I*phi)*sin(2*theta)/(8*sqrt(pi))
84 >>> simplify(Ynm(2, 2, theta, phi).expand(func=True))
85 sqrt(30)*exp(2*I*phi)*sin(theta)**2/(8*sqrt(pi))
87 We can differentiate the functions with respect
88 to both angles:
90 >>> from sympy import Ynm, Symbol, diff
91 >>> from sympy.abc import n,m
92 >>> theta = Symbol("theta")
93 >>> phi = Symbol("phi")
95 >>> diff(Ynm(n, m, theta, phi), theta)
96 m*cot(theta)*Ynm(n, m, theta, phi) + sqrt((-m + n)*(m + n + 1))*exp(-I*phi)*Ynm(n, m + 1, theta, phi)
98 >>> diff(Ynm(n, m, theta, phi), phi)
99 I*m*Ynm(n, m, theta, phi)
101 Further we can compute the complex conjugation:
103 >>> from sympy import Ynm, Symbol, conjugate
104 >>> from sympy.abc import n,m
105 >>> theta = Symbol("theta")
106 >>> phi = Symbol("phi")
108 >>> conjugate(Ynm(n, m, theta, phi))
109 (-1)**(2*m)*exp(-2*I*m*phi)*Ynm(n, m, theta, phi)
111 To get back the well known expressions in spherical
112 coordinates, we use full expansion:
114 >>> from sympy import Ynm, Symbol, expand_func
115 >>> from sympy.abc import n,m
116 >>> theta = Symbol("theta")
117 >>> phi = Symbol("phi")
119 >>> expand_func(Ynm(n, m, theta, phi))
120 sqrt((2*n + 1)*factorial(-m + n)/factorial(m + n))*exp(I*m*phi)*assoc_legendre(n, m, cos(theta))/(2*sqrt(pi))
122 See Also
123 ========
125 Ynm_c, Znm
127 References
128 ==========
130 .. [1] https://en.wikipedia.org/wiki/Spherical_harmonics
131 .. [2] https://mathworld.wolfram.com/SphericalHarmonic.html
132 .. [3] https://functions.wolfram.com/Polynomials/SphericalHarmonicY/
133 .. [4] https://dlmf.nist.gov/14.30
135 """
137 @classmethod
138 def eval(cls, n, m, theta, phi):
139 # Handle negative index m and arguments theta, phi
140 if m.could_extract_minus_sign():
141 m = -m
142 return S.NegativeOne**m * exp(-2*I*m*phi) * Ynm(n, m, theta, phi)
143 if theta.could_extract_minus_sign():
144 theta = -theta
145 return Ynm(n, m, theta, phi)
146 if phi.could_extract_minus_sign():
147 phi = -phi
148 return exp(-2*I*m*phi) * Ynm(n, m, theta, phi)
150 # TODO Add more simplififcation here
152 def _eval_expand_func(self, **hints):
153 n, m, theta, phi = self.args
154 rv = (sqrt((2*n + 1)/(4*pi) * factorial(n - m)/factorial(n + m)) *
155 exp(I*m*phi) * assoc_legendre(n, m, cos(theta)))
156 # We can do this because of the range of theta
157 return rv.subs(sqrt(-cos(theta)**2 + 1), sin(theta))
159 def fdiff(self, argindex=4):
160 if argindex == 1:
161 # Diff wrt n
162 raise ArgumentIndexError(self, argindex)
163 elif argindex == 2:
164 # Diff wrt m
165 raise ArgumentIndexError(self, argindex)
166 elif argindex == 3:
167 # Diff wrt theta
168 n, m, theta, phi = self.args
169 return (m * cot(theta) * Ynm(n, m, theta, phi) +
170 sqrt((n - m)*(n + m + 1)) * exp(-I*phi) * Ynm(n, m + 1, theta, phi))
171 elif argindex == 4:
172 # Diff wrt phi
173 n, m, theta, phi = self.args
174 return I * m * Ynm(n, m, theta, phi)
175 else:
176 raise ArgumentIndexError(self, argindex)
178 def _eval_rewrite_as_polynomial(self, n, m, theta, phi, **kwargs):
179 # TODO: Make sure n \in N
180 # TODO: Assert |m| <= n ortherwise we should return 0
181 return self.expand(func=True)
183 def _eval_rewrite_as_sin(self, n, m, theta, phi, **kwargs):
184 return self.rewrite(cos)
186 def _eval_rewrite_as_cos(self, n, m, theta, phi, **kwargs):
187 # This method can be expensive due to extensive use of simplification!
188 from sympy.simplify import simplify, trigsimp
189 # TODO: Make sure n \in N
190 # TODO: Assert |m| <= n ortherwise we should return 0
191 term = simplify(self.expand(func=True))
192 # We can do this because of the range of theta
193 term = term.xreplace({Abs(sin(theta)):sin(theta)})
194 return simplify(trigsimp(term))
196 def _eval_conjugate(self):
197 # TODO: Make sure theta \in R and phi \in R
198 n, m, theta, phi = self.args
199 return S.NegativeOne**m * self.func(n, -m, theta, phi)
201 def as_real_imag(self, deep=True, **hints):
202 # TODO: Handle deep and hints
203 n, m, theta, phi = self.args
204 re = (sqrt((2*n + 1)/(4*pi) * factorial(n - m)/factorial(n + m)) *
205 cos(m*phi) * assoc_legendre(n, m, cos(theta)))
206 im = (sqrt((2*n + 1)/(4*pi) * factorial(n - m)/factorial(n + m)) *
207 sin(m*phi) * assoc_legendre(n, m, cos(theta)))
208 return (re, im)
210 def _eval_evalf(self, prec):
211 # Note: works without this function by just calling
212 # mpmath for Legendre polynomials. But using
213 # the dedicated function directly is cleaner.
214 from mpmath import mp, workprec
215 n = self.args[0]._to_mpmath(prec)
216 m = self.args[1]._to_mpmath(prec)
217 theta = self.args[2]._to_mpmath(prec)
218 phi = self.args[3]._to_mpmath(prec)
219 with workprec(prec):
220 res = mp.spherharm(n, m, theta, phi)
221 return Expr._from_mpmath(res, prec)
224def Ynm_c(n, m, theta, phi):
225 r"""
226 Conjugate spherical harmonics defined as
228 .. math::
229 \overline{Y_n^m(\theta, \varphi)} := (-1)^m Y_n^{-m}(\theta, \varphi).
231 Examples
232 ========
234 >>> from sympy import Ynm_c, Symbol, simplify
235 >>> from sympy.abc import n,m
236 >>> theta = Symbol("theta")
237 >>> phi = Symbol("phi")
238 >>> Ynm_c(n, m, theta, phi)
239 (-1)**(2*m)*exp(-2*I*m*phi)*Ynm(n, m, theta, phi)
240 >>> Ynm_c(n, m, -theta, phi)
241 (-1)**(2*m)*exp(-2*I*m*phi)*Ynm(n, m, theta, phi)
243 For specific integers $n$ and $m$ we can evaluate the harmonics
244 to more useful expressions:
246 >>> simplify(Ynm_c(0, 0, theta, phi).expand(func=True))
247 1/(2*sqrt(pi))
248 >>> simplify(Ynm_c(1, -1, theta, phi).expand(func=True))
249 sqrt(6)*exp(I*(-phi + 2*conjugate(phi)))*sin(theta)/(4*sqrt(pi))
251 See Also
252 ========
254 Ynm, Znm
256 References
257 ==========
259 .. [1] https://en.wikipedia.org/wiki/Spherical_harmonics
260 .. [2] https://mathworld.wolfram.com/SphericalHarmonic.html
261 .. [3] https://functions.wolfram.com/Polynomials/SphericalHarmonicY/
263 """
264 return conjugate(Ynm(n, m, theta, phi))
267class Znm(Function):
268 r"""
269 Real spherical harmonics defined as
271 .. math::
273 Z_n^m(\theta, \varphi) :=
274 \begin{cases}
275 \frac{Y_n^m(\theta, \varphi) + \overline{Y_n^m(\theta, \varphi)}}{\sqrt{2}} &\quad m > 0 \\
276 Y_n^m(\theta, \varphi) &\quad m = 0 \\
277 \frac{Y_n^m(\theta, \varphi) - \overline{Y_n^m(\theta, \varphi)}}{i \sqrt{2}} &\quad m < 0 \\
278 \end{cases}
280 which gives in simplified form
282 .. math::
284 Z_n^m(\theta, \varphi) =
285 \begin{cases}
286 \frac{Y_n^m(\theta, \varphi) + (-1)^m Y_n^{-m}(\theta, \varphi)}{\sqrt{2}} &\quad m > 0 \\
287 Y_n^m(\theta, \varphi) &\quad m = 0 \\
288 \frac{Y_n^m(\theta, \varphi) - (-1)^m Y_n^{-m}(\theta, \varphi)}{i \sqrt{2}} &\quad m < 0 \\
289 \end{cases}
291 Examples
292 ========
294 >>> from sympy import Znm, Symbol, simplify
295 >>> from sympy.abc import n, m
296 >>> theta = Symbol("theta")
297 >>> phi = Symbol("phi")
298 >>> Znm(n, m, theta, phi)
299 Znm(n, m, theta, phi)
301 For specific integers n and m we can evaluate the harmonics
302 to more useful expressions:
304 >>> simplify(Znm(0, 0, theta, phi).expand(func=True))
305 1/(2*sqrt(pi))
306 >>> simplify(Znm(1, 1, theta, phi).expand(func=True))
307 -sqrt(3)*sin(theta)*cos(phi)/(2*sqrt(pi))
308 >>> simplify(Znm(2, 1, theta, phi).expand(func=True))
309 -sqrt(15)*sin(2*theta)*cos(phi)/(4*sqrt(pi))
311 See Also
312 ========
314 Ynm, Ynm_c
316 References
317 ==========
319 .. [1] https://en.wikipedia.org/wiki/Spherical_harmonics
320 .. [2] https://mathworld.wolfram.com/SphericalHarmonic.html
321 .. [3] https://functions.wolfram.com/Polynomials/SphericalHarmonicY/
323 """
325 @classmethod
326 def eval(cls, n, m, theta, phi):
327 if m.is_positive:
328 zz = (Ynm(n, m, theta, phi) + Ynm_c(n, m, theta, phi)) / sqrt(2)
329 return zz
330 elif m.is_zero:
331 return Ynm(n, m, theta, phi)
332 elif m.is_negative:
333 zz = (Ynm(n, m, theta, phi) - Ynm_c(n, m, theta, phi)) / (sqrt(2)*I)
334 return zz