杂化轨道

旧版：

        class HybridOrbital(VGroup):
            def __init__(
                self,
                hybrid_num: int,
                **kwargs
            ):
                super().__init__(**kwargs)
                property = self.hybrid_orbit_property(hybrid_num)
                for coe in property[3:]:
                    def hybrid_orbit(u, v):
                        n = 0
                        r = 0
                        for origin_orbit in property[2]:
                            r += coe[n]*self.harmonic_func(u, v, *origin_orbit)
                            n += 1
                        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)
                    orbit = Surface(
                        lambda u, v: axes.c2p(*hybrid_orbit(u, v)),
                        resolution=(20, 20),
                        v_range=[0, 2*PI],
                        u_range=[0, PI],
                        checkerboard_colors=[BLUE],
                        fill_opacity=0.8,
                        stroke_width=0,
                    )
                    self.add(orbit)

                if self.width >= self.height:
                    if self.width >= self.depth:
                        self.set(width=3)
                    else:
                        self.set(depth=3)
                else:
                    if self.height >= self.depth:
                        self.set(height=3)
                    else:
                        self.set(depth=3)
                    

                

            def hybrid_orbit_property(self, num):
                if num == 2:
                    property = ['sp','直线形',[[0,0],[1,1]],[np.sqrt(1/2), np.sqrt(1/2)],[np.sqrt(1/2), -np.sqrt(1/2)]]
                elif num == 3:
                    property = ['sp^{2}','平面正三角形',[[0,0],[1,1],[1,-1]],[np.sqrt(1/3),np.sqrt(2/3),0],[np.sqrt(1/3),-np.sqrt(1/6),np.sqrt(1/2)],[np.sqrt(1/3),-np.sqrt(1/6),-np.sqrt(1/2)]]
                elif num == 4:
                    property = ['sp^{3}','正四面体',[[0,0],[1,0],[1,1],[1,-1]],[1/2,1/2,1/2,1/2],[1/2,-1/2,-1/2,1/2],[1/2,1/2,-1/2,-1/2],[1/2,-1/2,1/2,-1/2]]
                elif num == 5:
                    property = ['dsp^{2}','平面正方形',[[0,0],[2,2],[1,1],[1,-1]],[1/2,1/2,1/2,1/2],[1/2,-1/2,1/2,-1/2],[np.sqrt(1/2),0,-np.sqrt(1/2),0],[0,np.sqrt(1/2),0,-np.sqrt(1/2)]]
                elif num == 6:
                    property = ['sp^{3}d  dsp^{3} d^{3}sp','三方双锥',[[0,0],[1,1],[1,-1],[1,0],[2,0]],[0,0,0,np.sqrt(1/2),np.sqrt(1/2)],[0,0,0,np.sqrt(1/2),-np.sqrt(1/2)],[np.sqrt(1/3),np.sqrt(1/3),np.sqrt(1/3),0,0],[np.sqrt(2/3),-np.sqrt(1/6),-np.sqrt(1/6),0,0],[0,np.sqrt(1/2),-np.sqrt(1/2),0,0]]
                elif num == 7:
                    property = ['dsp^{3}','四方锥',[[1,0],[0,0],[2,2],[1,1],[1,-1]],[1,0,0,0,0],[0,1/2,1/2,1/2,1/2],[0,1/2,-1/2,1/2,-1/2],[0,np.sqrt(1/2),0,-np.sqrt(1/2),0],[0,0,np.sqrt(1/2),0,-np.sqrt(1/2)]]
                elif num == 8:
                    property = ['sp^{3}d^{2}  d^{2}sp^{3}','八面体',[[0,0],[1,1],[1,-1],[1,0],[2,0],[2,2]],[np.sqrt(1/6),np.sqrt(1/6),np.sqrt(1/6),np.sqrt(1/6),np.sqrt(1/6),np.sqrt(1/6)],[-np.sqrt(1/12),-np.sqrt(1/12),-np.sqrt(1/12),-np.sqrt(1/12),np.sqrt(1/3),np.sqrt(1/3)],[1/2,-1/2,1/2,-1/2,0,0],[np.sqrt(1/2),0,-np.sqrt(1/2),0,0,0],[0,np.sqrt(1/2),0,-np.sqrt(1/2),0,0],[0,0,0,0,np.sqrt(1/2),-np.sqrt(1/2)]]
                elif num == 9:
                    property = ['sp^{3}d^{2}','三角锥',[[0,0],[1,0],[1,1],[1,-1],[2,1],[2,-1]],[np.sqrt(1/6),np.sqrt(1/6),np.sqrt(1/6),np.sqrt(1/6),np.sqrt(1/6),np.sqrt(1/6)],[np.sqrt(1/6),np.sqrt(1/6),np.sqrt(1/6),-np.sqrt(1/6),-np.sqrt(1/6),-np.sqrt(1/6)],[np.sqrt(1/3),np.sqrt(1/3),-np.sqrt(1/12),-np.sqrt(1/12),-np.sqrt(1/12),-np.sqrt(1/12)],[0,1/2,-1/2,0,1/2,-1/2],[np.sqrt(1/3),-np.sqrt(1/3),-np.sqrt(1/12),-np.sqrt(1/12),np.sqrt(1/12),np.sqrt(1/12)],[0,1/2,-1/2,0,-1/2,1/2]]
                return property
            
            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 harmonic_func(self, u, v, l, m):
                phi = self.wave_func_phi(m)
                theta = self.wave_func_theta(l, m)
                p = sym.symbols('p')
                t = sym.symbols('t')
                func = phi.subs(p, v)*theta.subs(t, u)  # u v 的顺序不要搞错
                return func

            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
            

但是使用波谐函数计算耗时太长，故直接使用球坐标下的角函数

改进后：

        class HybridOrbital(VGroup):
            def __init__(
                self,
                hybrid_num: int,
                **kwargs
            ):
                super().__init__(**kwargs)
                self.num = hybrid_num
                property = self.hybrid_orbit_property(hybrid_num)
                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)
                orbits = VGroup()
                labels = VGroup()
                for coe in property[4:]:
                    def hybrid_orbit(u, v):
                        r = 0
                        n = 0
                        for origin_orbit in property[3]:
                            r += coe[n]*self.harmonic_func(u, v, *origin_orbit)
                            n += 1
                        r = r ** 2
                        return [r*np.sin(u)*np.cos(v), r*np.sin(u)*np.sin(v), r*np.cos(u)]
                    orbit = Surface(
                        lambda u, v: axes.c2p(*hybrid_orbit(u, v)),
                        resolution=(40, 40),
                        v_range=[0, 2*PI],
                        u_range=[0, PI],
                        checkerboard_colors=[BLUE],
                        fill_opacity=0.8,
                        stroke_width=0,
                    )
                    orbits.add(orbit)
                self.scale_to_fit(orbits)
                geometry = self.geometry()
                self.add(orbits, geometry)
                self.arrange(buff=2)
                label1 = TextFlower(property[0], color_type=1, use_latex=True).to_corner(UL)
                label2 = TextFlower(property[1], color_type=2, use_latex=True).next_to(label1, RIGHT)
                label3 = TextFlower(property[2], color_type=3).next_to(label2, RIGHT)
                labels.add(label1, label2, label3)
                self.add(labels)  

            def scale_to_fit(self, mob):
                if mob.width >= mob.height:
                    if mob.width >= mob.depth:
                        mob.set(width=4.5)
                    else:
                        mob.set(depth=4.5)
                else:
                    if mob.height >= mob.depth:
                        mob.set(height=4.5)
                    else:
                        mob.set(depth=4.5)

            def hybrid_orbit_property(self, num):
                if num == 2:
                    property = ['AB_{2}','sp','直线形',[[0,0],[1,1]],[np.sqrt(1/2), np.sqrt(1/2)],[np.sqrt(1/2), -np.sqrt(1/2)]]
                elif num == 3:
                    property = ['AB_{3}','sp^{2}','平面正三角形',[[0,0],[1,1],[1,-1]],[np.sqrt(1/3),np.sqrt(2/3),0],[np.sqrt(1/3),-np.sqrt(1/6),np.sqrt(1/2)],[np.sqrt(1/3),-np.sqrt(1/6),-np.sqrt(1/2)]]
                elif num == 4:
                    property = ['AB_{4}','sp^{3}','正四面体',[[0,0],[1,0],[1,1],[1,-1]],[1/2,1/2,1/2,1/2],[1/2,-1/2,-1/2,1/2],[1/2,1/2,-1/2,-1/2],[1/2,-1/2,1/2,-1/2]]
                elif num == 5:
                    property = ['AB_{4}','dsp^{2}','平面正方形',[[0,0],[2,2],[1,1],[1,-1]],[1/2,1/2,np.sqrt(1/2),0],[1/2,-1/2,0,np.sqrt(1/2)],[1/2,1/2,-np.sqrt(1/2),0],[1/2,-1/2,0,-np.sqrt(1/2)]]
                elif num == 6:
                    property = ['AB_{5}','dsp^{3}/d^{3}sp','三角双锥',[[0,0],[2,0],[1,1],[1,-1],[1,0]],[np.sqrt(1/2),0,0,0,np.sqrt(1/2)],[np.sqrt(1/2),0,0,0,-np.sqrt(1/2)],[0,np.sqrt(1/3),np.sqrt(2/3),0,0],[0,np.sqrt(1/3),-np.sqrt(1/6),np.sqrt(1/2),0],[0,np.sqrt(1/3),-np.sqrt(1/6),-np.sqrt(1/2),0]]
                elif num == 7:
                    property = ['AB_{6}','d^{4}sp','正三棱柱',[[0,0],[1,1],[1,-1],[1,0],[2,1],[2,-1]],[np.sqrt(1/6),np.sqrt(1/3),0,-np.sqrt(1/6),-np.sqrt(1/3),0],[np.sqrt(1/6),np.sqrt(1/3),0,np.sqrt(1/6),np.sqrt(1/3),0],[np.sqrt(1/6),-np.sqrt(1/12),1/2,-np.sqrt(1/6),np.sqrt(1/12),-1/2],[np.sqrt(1/6),-np.sqrt(1/12),1/2,np.sqrt(1/6),-np.sqrt(1/12),1/2],[np.sqrt(1/6),-np.sqrt(1/12),-1/2,np.sqrt(1/6),-np.sqrt(1/12),-1/2],[np.sqrt(1/6),-np.sqrt(1/12),-1/2,-np.sqrt(1/6),np.sqrt(1/12),1/2]]
                elif num == 8:
                    property = ['AB_{6}','sp^{3}d^{2}/d^{2}sp^{3}','正八面体',[[0,0],[2,0],[2,2],[1,1],[1,-1],[1,0]],[np.sqrt(1/6),-np.sqrt(1/12),1/2,np.sqrt(1/2),0,0],[np.sqrt(1/6),-np.sqrt(1/12),1/2,-np.sqrt(1/2),0,0],[np.sqrt(1/6),-np.sqrt(1/12),-1/2,0,np.sqrt(1/2),0],[np.sqrt(1/6),-np.sqrt(1/12),-1/2,0,-np.sqrt(1/2),0],[np.sqrt(1/6),np.sqrt(1/3),0,0,0,np.sqrt(1/2)],[np.sqrt(1/6),np.sqrt(1/3),0,0,0,-np.sqrt(1/2)]]
                elif num == 9:
                    sin1=np.sin(0.4*PI)
                    sin2=np.sin(0.8*PI)
                    cos1=np.cos(0.4*PI)
                    cos2=np.cos(0.8*PI)
                    s1 = sin1/np.sqrt(2*(sin1**2+sin2**2))
                    s2 = sin2/np.sqrt(2*(sin1**2+sin2**2))
                    c1 = cos1/np.sqrt(1+2*(cos1**2+cos2**2))
                    c2 = cos2/np.sqrt(1+2*(cos1**2+cos2**2))
                    constant = 1/np.sqrt(1+2*(cos1**2+cos2**2))
                    property = ['AB_{7}','d^{3}sp^{3}','五角双锥',[[2,0],[0,0],[1,1],[1,-1],[2,2],[2,-2],[1,0]],[np.sqrt(1/7),np.sqrt(1/2),0,0,0,0,np.sqrt(1/2)],[np.sqrt(1/7),np.sqrt(1/2),0,0,0,0,-np.sqrt(1/2)],[np.sqrt(1/7),0,0,constant,constant,0,0],[np.sqrt(1/7),0,s2,c2,c1,s1,0],[np.sqrt(1/7),0,-s1,c1,c2,s2,0],[np.sqrt(1/7),0,s1,c1,c2,-s2,0],[np.sqrt(1/7),0,-s2,c2,c1,-s1,0]]
                elif num == 10:
                    c = np.sqrt(1/8)
                    property = ['AB_{8}','d^{3}fsp^{3}','立方体',[[0,0],[2,1],[2,-1],[2,-2],[3,2],[1,1],[1,-1],[1,0]],[c,c,c,c,c,c,c,c],[c,-c,c,-c,-c,-c,c,c],[c,-c,-c,c,c,-c,-c,c],[c,c,-c,-c,-c,c,-c,c],[c,-c,-c,c,-c,c,c,-c],[c,c,-c,-c,c,-c,c,-c],[c,c,c,c,-c,-c,-c,-c],[c,-c,c,-c,c,c,-c,-c]]
                elif num == 11:
                    c = np.sqrt(1/8)
                    property = ['AB_{8}','d^{4}sp^{3}','正四角反棱柱',[[0,0],[1,0],[1,1],[1,-1],[2,2],[2,-2],[2,1],[2,-1]],[c,-c,c,-c,0,-1/2,-c,c],[c,-c,c,c,0,1/2,-c,-c],[c,-c,-c,c,0,-1/2,c,-c],[c,-c,-c,-c,0,1/2,c,c],[c,c,1/2,0,1/2,0,1/2,0],[c,c,0,1/2,-1/2,0,0,1/2],[c,c,-1/2,0,1/2,0,-1/2,0],[c,c,0,-1/2,-1/2,0,0,-1/2]]
                elif num == 12:
                    property = ['AB_{8}','d^{4}sp^{3}','十二面体',[[0,0],[2,0],[2,2],[1,0],[1,1],[1,-1],[2,1],[2,-1]],[1/2,0,1/2,0,np.sqrt(1/2),0,0,0],[1/2,0,-1/2,0,0,np.sqrt(1/2),0,0],[1/2,0,1/2,0,-np.sqrt(1/2),0,0,0],[1/2,0,-1/2,0,0,-np.sqrt(1/2),0,0],[0,1/2,0,1/2,0,0,np.sqrt(1/2),0],[0,1/2,0,1/2,0,0,-np.sqrt(1/2),0],[0,1/2,0,-1/2,0,0,0,np.sqrt(1/2)],[0,1/2,0,-1/2,0,0,0,-np.sqrt(1/2)]]
                return property

            def harmonic_func(self, u, v, l, m):
                if l == 0 & m == 0:
                    r = np.sqrt(1/(4*PI))
                elif l == 1:
                    if m == 0:
                        r = np.sqrt(3/(4*PI))*np.cos(u)
                    elif m == 1:
                        r = np.sqrt(3/(4*PI))*np.sin(u)*np.cos(v)
                    elif m == -1:
                        r = np.sqrt(3/(4*PI))*np.sin(u)*np.sin(v)
                elif l == 2:
                    if m == 0:
                        r = np.sqrt(5/(16*PI))*(3*np.cos(u)**2 - 1)
                    if m == 1:
                        r = np.sqrt(15/(16*PI))*np.sin(2*u)*np.cos(v)
                    if m == -1:
                        r = np.sqrt(15/(16*PI))*np.sin(2*u)*np.sin(v)
                    if m == 2:
                        r = np.sqrt(15/(16*PI))*np.sin(u)**2*np.cos(2*v)
                    if m == -2:
                        r = np.sqrt(15/(16*PI))*np.sin(u)**2*np.sin(2*v)
                elif l == 3:
                    if m == 0:
                        r = 1/4*np.sqrt(7/PI)*(5*np.cos(u)**3 - 3*np.cos(u))
                    if m == 1:
                        r = 1/8*np.sqrt(42/PI)*np.sin(u)*(5*np.cos(u)**2 - 1)*np.cos(v)
                    if m == -1:
                        r = 1/8*np.sqrt(42/PI)*np.sin(u)*(5*np.cos(u)**2 - 1)*np.sin(v)
                    if m == 2:
                        r = 1/4*np.sqrt(105/PI)*np.sin(u)**2*np.cos(u)*np.sin(2*v)
                    if m == -2:
                        r = 1/4*np.sqrt(105/PI)*np.sin(u)**2*np.cos(u)*np.cos(2*v)
                    if m == 3:
                        r = 1/8*np.sqrt(70/PI)*np.sin(u)**3*np.cos(3*v)
                    if m == -3:
                        r = 1/8*np.sqrt(70/PI)*np.sin(u)**3*np.sin(3*v)
                else:
                    raise Exception('l or m is too large')
                return r

            def geometry(self):
                num = self.num
                if num == 2:
                    vertex_coords = [[1,0,0],[-1,0,0]]
                    faces_list = [[0,1]]
                elif num == 3:
                    vertex_coords = [[2,0,0],[-1,np.sqrt(3),0],[-1,-np.sqrt(3),0]]
                    faces_list = [[0,1,2]]
                elif num == 4:
                    vertex_coords = [[1,1,1],[-1,-1,1],[1,-1,-1],[-1,1,-1]]
                    faces_list = [[0,1,2],[0,1,3],[0,2,3],[1,2,3]]
                elif num == 5:
                    vertex_coords = [[1,0,0],[0,1,0],[-1,0,0],[0,-1,0]]
                    faces_list = [[0,1,2,3]]
                elif num == 6:
                    vertex_coords = [[-2,0,0],[1,np.sqrt(3),0],[1,-np.sqrt(3),0],[0,0,2],[0,0,-2]]
                    faces_list = [[0,1,3],[1,2,3],[0,2,3],[0,1,4],[1,2,4],[0,2,4]]
                elif num == 7:
                    vertex_coords = [[-1,np.sqrt(3),np.sqrt(3)],[2,0,np.sqrt(3)],[-1,-np.sqrt(3),np.sqrt(3)],[-1,np.sqrt(3),-np.sqrt(3)],[2,0,-np.sqrt(3)],[-1,-np.sqrt(3),-np.sqrt(3)]]
                    faces_list = [[0,1,2],[0,2,5,3],[2,1,4,5],[0,1,4,3],[3,4,5]]
                elif num == 8:
                    vertex_coords = [[0,-1,0],[-1,0,0],[0,1,0],[1,0,0],[0,0,1],[0,0,-1]]
                    faces_list = [[0,1,4],[1,2,4],[2,3,4],[0,3,4],[0,1,5],[1,2,5],[2,3,5],[0,3,5]]
                elif num == 9:
                    vertex_coords = [[0,-1,0],[np.cos(0.1*PI),-np.sin(0.1*PI),0],[np.cos(0.3*PI),np.sin(0.3*PI),0],[-np.cos(0.3*PI),np.sin(0.3*PI),0],[-np.cos(0.1*PI),-np.sin(0.1*PI),0],[0,0,1],[0,0,-1]]
                    faces_list = [[5,4,0],[5,0,1],[5,1,2],[5,2,3],[5,3,4],[5,4,0],[6,4,0],[6,0,1],[6,1,2],[6,2,3],[6,3,4],[6,4,0]]
                elif num == 10:
                    vertex_coords = [[1,1,1],[1,-1,1],[-1,-1,1],[-1,1,1],[1,1,-1],[1,-1,-1],[-1,-1,-1],[-1,1,-1]]
                    faces_list = [[0,1,2,3],[3,0,4,7],[0,1,5,4],[2,1,5,6],[2,3,7,6],[4,5,6,7]]
                elif num == 11:
                    vertex_coords = [[0,np.sqrt(2),1],[np.sqrt(2),0,1],[0,-np.sqrt(2),1],[-np.sqrt(2),0,1],[1,1,-1],[1,-1,-1],[-1,-1,-1],[-1,1,-1]]
                    faces_list = [[0,1,2,3],[3,0,7],[0,1,4],[1,4,5],[1,2,5],[2,5,6],[2,3,6],[3,6,7],[4,5,6,7]]
                elif num == 12:
                    vertex_coords = [[0,2,0],[2,0,0],[0,-2,0],[-2,0,0],[-1,0,np.sqrt(3)],[1,0,np.sqrt(3)],[0,1,-np.sqrt(3)],[0,-1,-np.sqrt(3)]]
                    faces_list = [[0,3,4],[0,1,5],[0,4,5],[2,3,4],[1,2,5],[2,4,5],[0,3,6],[0,1,6],[2,3,7],[1,2,7],[3,6,7],[1,6,7]]
                geo = Polyhedron(vertex_coords, faces_list)
                self.scale_to_fit(geo)
                for dot in geo.graph:
                    dot.set(width=0.3)
                
                return geo