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

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 

12 

13_x = Dummy("x") 

14 

15class Ynm(Function): 

16 r""" 

17 Spherical harmonics defined as 

18 

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) 

23 

24 Explanation 

25 =========== 

26 

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]$. 

32 

33 Examples 

34 ======== 

35 

36 >>> from sympy import Ynm, Symbol, simplify 

37 >>> from sympy.abc import n,m 

38 >>> theta = Symbol("theta") 

39 >>> phi = Symbol("phi") 

40 

41 >>> Ynm(n, m, theta, phi) 

42 Ynm(n, m, theta, phi) 

43 

44 Several symmetries are known, for the order: 

45 

46 >>> Ynm(n, -m, theta, phi) 

47 (-1)**m*exp(-2*I*m*phi)*Ynm(n, m, theta, phi) 

48 

49 As well as for the angles: 

50 

51 >>> Ynm(n, m, -theta, phi) 

52 Ynm(n, m, theta, phi) 

53 

54 >>> Ynm(n, m, theta, -phi) 

55 exp(-2*I*m*phi)*Ynm(n, m, theta, phi) 

56 

57 For specific integers $n$ and $m$ we can evaluate the harmonics 

58 to more useful expressions: 

59 

60 >>> simplify(Ynm(0, 0, theta, phi).expand(func=True)) 

61 1/(2*sqrt(pi)) 

62 

63 >>> simplify(Ynm(1, -1, theta, phi).expand(func=True)) 

64 sqrt(6)*exp(-I*phi)*sin(theta)/(4*sqrt(pi)) 

65 

66 >>> simplify(Ynm(1, 0, theta, phi).expand(func=True)) 

67 sqrt(3)*cos(theta)/(2*sqrt(pi)) 

68 

69 >>> simplify(Ynm(1, 1, theta, phi).expand(func=True)) 

70 -sqrt(6)*exp(I*phi)*sin(theta)/(4*sqrt(pi)) 

71 

72 >>> simplify(Ynm(2, -2, theta, phi).expand(func=True)) 

73 sqrt(30)*exp(-2*I*phi)*sin(theta)**2/(8*sqrt(pi)) 

74 

75 >>> simplify(Ynm(2, -1, theta, phi).expand(func=True)) 

76 sqrt(30)*exp(-I*phi)*sin(2*theta)/(8*sqrt(pi)) 

77 

78 >>> simplify(Ynm(2, 0, theta, phi).expand(func=True)) 

79 sqrt(5)*(3*cos(theta)**2 - 1)/(4*sqrt(pi)) 

80 

81 >>> simplify(Ynm(2, 1, theta, phi).expand(func=True)) 

82 -sqrt(30)*exp(I*phi)*sin(2*theta)/(8*sqrt(pi)) 

83 

84 >>> simplify(Ynm(2, 2, theta, phi).expand(func=True)) 

85 sqrt(30)*exp(2*I*phi)*sin(theta)**2/(8*sqrt(pi)) 

86 

87 We can differentiate the functions with respect 

88 to both angles: 

89 

90 >>> from sympy import Ynm, Symbol, diff 

91 >>> from sympy.abc import n,m 

92 >>> theta = Symbol("theta") 

93 >>> phi = Symbol("phi") 

94 

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) 

97 

98 >>> diff(Ynm(n, m, theta, phi), phi) 

99 I*m*Ynm(n, m, theta, phi) 

100 

101 Further we can compute the complex conjugation: 

102 

103 >>> from sympy import Ynm, Symbol, conjugate 

104 >>> from sympy.abc import n,m 

105 >>> theta = Symbol("theta") 

106 >>> phi = Symbol("phi") 

107 

108 >>> conjugate(Ynm(n, m, theta, phi)) 

109 (-1)**(2*m)*exp(-2*I*m*phi)*Ynm(n, m, theta, phi) 

110 

111 To get back the well known expressions in spherical 

112 coordinates, we use full expansion: 

113 

114 >>> from sympy import Ynm, Symbol, expand_func 

115 >>> from sympy.abc import n,m 

116 >>> theta = Symbol("theta") 

117 >>> phi = Symbol("phi") 

118 

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

121 

122 See Also 

123 ======== 

124 

125 Ynm_c, Znm 

126 

127 References 

128 ========== 

129 

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 

134 

135 """ 

136 

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) 

149 

150 # TODO Add more simplififcation here 

151 

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

158 

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) 

177 

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) 

182 

183 def _eval_rewrite_as_sin(self, n, m, theta, phi, **kwargs): 

184 return self.rewrite(cos) 

185 

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

195 

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) 

200 

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) 

209 

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) 

222 

223 

224def Ynm_c(n, m, theta, phi): 

225 r""" 

226 Conjugate spherical harmonics defined as 

227 

228 .. math:: 

229 \overline{Y_n^m(\theta, \varphi)} := (-1)^m Y_n^{-m}(\theta, \varphi). 

230 

231 Examples 

232 ======== 

233 

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) 

242 

243 For specific integers $n$ and $m$ we can evaluate the harmonics 

244 to more useful expressions: 

245 

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

250 

251 See Also 

252 ======== 

253 

254 Ynm, Znm 

255 

256 References 

257 ========== 

258 

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/ 

262 

263 """ 

264 return conjugate(Ynm(n, m, theta, phi)) 

265 

266 

267class Znm(Function): 

268 r""" 

269 Real spherical harmonics defined as 

270 

271 .. math:: 

272 

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} 

279 

280 which gives in simplified form 

281 

282 .. math:: 

283 

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} 

290 

291 Examples 

292 ======== 

293 

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) 

300 

301 For specific integers n and m we can evaluate the harmonics 

302 to more useful expressions: 

303 

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

310 

311 See Also 

312 ======== 

313 

314 Ynm, Ynm_c 

315 

316 References 

317 ========== 

318 

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/ 

322 

323 """ 

324 

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