原子轨道（角向部分）

        class WaveFunc(VGroup):
            def __init__(
                self,
                l: int,
                m: int,
                **kwargs
            ):
                super().__init__(**kwargs)
                self.l = l
                self.m = m
                phi = self.wave_func_phi(m)
                theta = self.wave_func_theta(l, m)
                p = sym.symbols('p')
                t = sym.symbols('t')
                def harmonic_func(u, v):
                    r = phi.subs(p, v)*theta.subs(t, u)  # u v 的顺序不要搞错
                    r = r ** 2
                    return [r*np.sin(u)*np.cos(v), r*np.sin(u)*np.sin(v), r*np.cos(u)]
                axes = ThreeDAxes(x_range=(-5, 5, 1), y_range=(-5, 5, 1), z_range=(-5, 5, 1),
                                  x_length=5, y_length=5, z_length=5)
                nucleus = Sphere(fill_opacity=0.7, stroke_width=0, radius=0.5).set_color("#9AFF9A")
                orbit = Surface(
                    lambda u, v: axes.c2p(*harmonic_func(u, v)),
                    resolution=(40, 40),
                    v_range=[0, 2*PI],
                    u_range=[0, PI],
                    checkerboard_colors=[BLUE],
                    fill_opacity=0.8,
                    stroke_width=0,
                )
                nucleus_radius = nucleus.width
                if orbit.width >= orbit.height:
                    if orbit.width >= orbit.depth:
                        orbit.set(width=6*nucleus_radius)
                    else:
                        orbit.set(depth=6*nucleus_radius)
                else:
                    if orbit.height >= orbit.depth:
                        orbit.set(height=6*nucleus_radius)
                    else:
                        orbit.set(depth=6*nucleus_radius)

                self.orbit = orbit
                self.l_and_m = f'l = {str(l)}, m = {str(m)}'
                self.add(orbit)
            
            def wave_func_phi(self, m):
                p = sym.symbols('p')
                if m >= 0:
                    result = sym.cos(m*p)  # 省略了系数
                elif m < 0:
                    result = sym.sin(np.abs(m)*p)  # 省略了系数
                return result

            def wave_func_theta(self, l, m):
                m = np.abs(m)
                x = sym.symbols('x')
                t = sym.symbols('t')
                def associated_legendre_polynomial(l, m):
                    def legendre_polynomial(n):
                        if n == 0:
                            return 1
                        elif n == 1:
                            return x
                        elif n >= 2:
                            result = (2*n-1)/n*x*legendre_polynomial(n-1) - (n-1)/n*legendre_polynomial(n-2)
                            return sym.simplify(result)
                    
                    result = (1-x**2)**(m/2)*sym.diff(legendre_polynomial(l), x, m)
                    return sym.simplify(result.subs(x, sym.cos(t)))

                result = associated_legendre_polynomial(l, m)  # 省略了系数
                return result

            def wave_func_radius(self, n, l):
                x = sym.symbols('x')
                def associated_Laguerre_polynomial(n, l):
                    result = x**(-l)*sym.exp(x)*sym.diff(sym.exp(-x)*x**(n+l), x, n)  # 这里的 x 实际上也要乘以一个系数
                    return sym.simplify(result)

                result = sym.exp(-x/2)*x**l*associated_Laguerre_polynomial(n, l)  # 省略了系数
                return result
            
            def axes(self):
                axes = SVGMobject(r"D:\manimSVG\axes.svg").set_height(5.5)
                axes_distance_deviation = ORIGIN - axes[-1].get_center()
                axes.shift(axes_distance_deviation)

                return axes