Skip to content

Commit fd73a70

Browse files
haakoekhaakoek
committed
Fix compability sph_harm to sph_harm_y for SciPy versions >= 1.15.0.
1 parent 713adb2 commit fd73a70

1 file changed

Lines changed: 84 additions & 51 deletions

File tree

‎grid_lib/spherical_coordinates/angular_matrix_elements.py‎

Lines changed: 84 additions & 51 deletions
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,47 @@
11
import numpy as np
22
import time
33
from numpy.polynomial import legendre
4-
from scipy.special import eval_legendre as Legendre, sph_harm_y, spherical_jn
4+
from scipy.special import spherical_jn
55
from pathlib import Path
66

7+
import scipy.special
8+
from packaging import version
9+
10+
11+
def sph_harm_y(m, n, phi, theta):
12+
"""
13+
Args:
14+
m: azimuthal quantum number
15+
n: Degree of the harmonic, commonly denoted as l in quantum mechanics
16+
theta: polar angle in [0,pi]
17+
phi: azimuthal angle in [0, 2pi)
18+
--------------------------------------------------------------------------------------------------------
19+
--------------------------------------------------------------------------------------------------------
20+
In Scipy < 1.15.0
21+
sph_harm: sph_harm(m, n, theta, phi, out=None)
22+
- m: azimuthal quantum number
23+
- n: Degree of the harmonic, commonly denoted as l in quantum mechanics
24+
- theta in [0, 2pi): azimuthal angle
25+
- phi in [0,pi]: polar angle
26+
27+
In SciPy >= 1.15.0 sph_harm is deprecated (and will be removed in SciPy 1.17.0) called sph_harm_y where
28+
sph_harm_y: sph_harm_y(n, m, theta, phi, *, diff_n=0)
29+
- m: azimuthal quantum number
30+
- n: Degree of the harmonic, commonly denoted as l in quantum mechanics
31+
- theta in [0,pi]: polar angle
32+
- phi in [0, 2pi): azimuthal angle
33+
--------------------------------------------------------------------------------------------------------
34+
--------------------------------------------------------------------------------------------------------
35+
"""
36+
scipy_version = version.parse(scipy.__version__)
37+
38+
if scipy_version >= version.parse("1.15.0"):
39+
# New: sph_harm_y(degree, order, polar, azimuthal)
40+
return scipy.special.sph_harm_y(n, m, theta, phi)
41+
else:
42+
# Old: sph_harm(order, degree, azimuthal, polar)
43+
return scipy.special.sph_harm(m, n, phi, theta)
44+
745

846
class AngularMatrixElements:
947
def __init__(self, l_max=5, N=101):
@@ -34,9 +72,7 @@ def kron_delta(self, i, j):
3472

3573
def a_lm(self, l, m):
3674
if l >= 0 and abs(m) <= l:
37-
return np.sqrt(
38-
((l + 1) ** 2 - m**2) / ((2 * l + 1) * (2 * l + 3))
39-
)
75+
return np.sqrt(((l + 1) ** 2 - m**2) / ((2 * l + 1) * (2 * l + 3)))
4076
else:
4177
return 0
4278

@@ -530,19 +566,16 @@ def setup_l_matrix_elements(self, arr_to_calc_dict, m=0):
530566
and arr_to_calc_dict["z_Omega"] == False
531567
):
532568
if abs(m) <= l1 and abs(m) <= l2:
533-
self.arr["H_z_beta"][
534-
l1, l2
535-
] = -self.l1m1_sinth_ddtheta_l2m2(
536-
l1, m, l2, m
537-
) - self.l1m1_costh_l2m2(
538-
l1, m, l2, m
569+
self.arr["H_z_beta"][l1, l2] = (
570+
-self.l1m1_sinth_ddtheta_l2m2(l1, m, l2, m)
571+
- self.l1m1_costh_l2m2(l1, m, l2, m)
539572
)
540573
if arr_to_calc_dict["H_Bz_Omega"]:
541574
if abs(m) <= l1 and abs(m) <= l2:
542-
self.arr["H_Bz_Omega"][
543-
l1, l2
544-
] = self.l1m1_sinth_sq_l2m2_Lebedev(
545-
Yl1m_cc, Yl2m, l1, m, l2, m
575+
self.arr["H_Bz_Omega"][l1, l2] = (
576+
self.l1m1_sinth_sq_l2m2_Lebedev(
577+
Yl1m_cc, Yl2m, l1, m, l2, m
578+
)
546579
)
547580

548581

@@ -594,60 +627,60 @@ def setup_matrix_elements(self, arr_to_calc_dict):
594627
Yl2m2 = sph_harm_y(m2, l2, self.phi, self.theta)
595628

596629
if arr_to_calc_dict["x_Omega"]:
597-
self.arr["x_Omega"][
598-
I, J
599-
] = self.l1m1_sinth_cosph_l2m2(l1, m1, l2, m2)
630+
self.arr["x_Omega"][I, J] = (
631+
self.l1m1_sinth_cosph_l2m2(l1, m1, l2, m2)
632+
)
600633

601634
if arr_to_calc_dict["x_x_Omega"]:
602-
self.arr["x_x_Omega"][
603-
I, J
604-
] = self.l1m1_sinth_sq_cosph_sq_l2m2_Lebedev(
605-
Yl1m1_cc, Yl2m2, l1, m1, l2, m2
635+
self.arr["x_x_Omega"][I, J] = (
636+
self.l1m1_sinth_sq_cosph_sq_l2m2_Lebedev(
637+
Yl1m1_cc, Yl2m2, l1, m1, l2, m2
638+
)
606639
)
607640

608641
if arr_to_calc_dict["y_Omega"]:
609-
self.arr["y_Omega"][
610-
I, J
611-
] = self.l1m1_sinth_sinph_l2m2(l1, m1, l2, m2)
642+
self.arr["y_Omega"][I, J] = (
643+
self.l1m1_sinth_sinph_l2m2(l1, m1, l2, m2)
644+
)
612645

613646
if arr_to_calc_dict["y_y_Omega"]:
614-
self.arr["y_y_Omega"][
615-
I, J
616-
] = self.l1m1_sinth_sq_sinph_sq_l2m2_Lebedev(
617-
Yl1m1_cc, Yl2m2, l1, m1, l2, m2
647+
self.arr["y_y_Omega"][I, J] = (
648+
self.l1m1_sinth_sq_sinph_sq_l2m2_Lebedev(
649+
Yl1m1_cc, Yl2m2, l1, m1, l2, m2
650+
)
618651
)
619652

620653
if arr_to_calc_dict["z_Omega"]:
621-
self.arr["z_Omega"][
622-
I, J
623-
] = self.l1m1_costh_l2m2(l1, m1, l2, m2)
654+
self.arr["z_Omega"][I, J] = (
655+
self.l1m1_costh_l2m2(l1, m1, l2, m2)
656+
)
624657

625658
if arr_to_calc_dict["z_z_Omega"]:
626-
self.arr["z_z_Omega"][
627-
I, J
628-
] = self.l1m1_costh_sq_l2m2_Lebedev(
629-
Yl1m1_cc, Yl2m2, l1, m1, l2, m2
659+
self.arr["z_z_Omega"][I, J] = (
660+
self.l1m1_costh_sq_l2m2_Lebedev(
661+
Yl1m1_cc, Yl2m2, l1, m1, l2, m2
662+
)
630663
)
631664

632665
if arr_to_calc_dict["z_x_Omega"]:
633-
self.arr["z_x_Omega"][
634-
I, J
635-
] = self.l1m1_costh_sinth_cosph_l2m2_Lebedev(
636-
Yl1m1_cc, Yl2m2, l1, m1, l2, m2
666+
self.arr["z_x_Omega"][I, J] = (
667+
self.l1m1_costh_sinth_cosph_l2m2_Lebedev(
668+
Yl1m1_cc, Yl2m2, l1, m1, l2, m2
669+
)
637670
)
638671

639672
if arr_to_calc_dict["z_y_Omega"]:
640-
self.arr["z_y_Omega"][
641-
I, J
642-
] = self.l1m1_costh_sinth_sinph_l2m2_Lebedev(
643-
Yl1m1_cc, Yl2m2, l1, m1, l2, m2
673+
self.arr["z_y_Omega"][I, J] = (
674+
self.l1m1_costh_sinth_sinph_l2m2_Lebedev(
675+
Yl1m1_cc, Yl2m2, l1, m1, l2, m2
676+
)
644677
)
645678

646679
if arr_to_calc_dict["y_x_Omega"]:
647-
self.arr["y_x_Omega"][
648-
I, J
649-
] = self.l1m1_sinph_cosph_sinthsq_l2m2(
650-
Yl1m1_cc, Yl2m2, l1, m1, l2, m2
680+
self.arr["y_x_Omega"][I, J] = (
681+
self.l1m1_sinph_cosph_sinthsq_l2m2(
682+
Yl1m1_cc, Yl2m2, l1, m1, l2, m2
683+
)
651684
)
652685

653686
if arr_to_calc_dict["y_px_beta"]:
@@ -723,10 +756,10 @@ def setup_matrix_elements(self, arr_to_calc_dict):
723756
)
724757

725758
if arr_to_calc_dict["H_Bz_Omega"]:
726-
self.arr["H_Bz_Omega"][
727-
I, J
728-
] = self.l1m1_sinth_sq_l2m2_Lebedev(
729-
Yl1m1_cc, Yl2m2, l1, m1, l2, m2
759+
self.arr["H_Bz_Omega"][I, J] = (
760+
self.l1m1_sinth_sq_l2m2_Lebedev(
761+
Yl1m1_cc, Yl2m2, l1, m1, l2, m2
762+
)
730763
)
731764
toc = time.time()
732765
print(f"Time setup angular matrix elements: {toc-tic}")

0 commit comments

Comments
 (0)