python实现完全规格化缔合勒让德（Legendre）函数

说在前面

关于怎样实现完全正规化勒让德函数的计算，在网上算法是有的，但是我发现用python具体实现的情况比较少，所以上传我的代码给有需求的人参考一下。

代码实例

#用Phi和阶数计算完全规格化勒让德函数，Phi=pi/2-Theta，这里的索引是1-181对应P0-180，故用182阶方阵
def Pnm(Phi,Degree):    
    P=np.zeros([Degree+2,Degree+2])
    #跨阶次正规化勒让德系数
    P[1][1] = 1
    P[2][1] = sin(Phi)*3**0.5
    P[2][2] = sqrt(3*(1-sin(Phi)**2))

    for j in range(1,3):
        for i in range(3,Degree+2):
            l=i-1
            m=j-1
            a = sqrt((4*l**2-1)/(l**2-m**2))
            b = sqrt((2*l+1)/(2*l-3))*sqrt(((l-1)**2-m**2)/(l**2-m**2))
            P[i][j] = a*sin(Phi)*P[i-1][j]-b*P[i-2][j]

    for j in range(3,Degree+1):
        for i in range(j,j+2):
            l=i-1
            m=j-1
            alpha = sqrt((2*l+1)*(l-m)*(l-m-1)/(2*l-3)/(l+m)/(l+m-1))
            if(m==2):
                beta = sqrt(2*(2*l+1)*(l+m-2)*(l+m-3)/(2*l-3)/(l+m)/(l+m-1))
                gama = sqrt(2*(l-m+1)*(l-m+2)/(l+m)/(l+m-1))
            else:
                beta = sqrt((2*l+1)*(l+m-2)*(l+m-3)/(2*l-3)/(l+m)/(l+m-1))
                gama = sqrt((l-m+1)*(l-m+2)/(l+m)/(l+m-1))
            P[i][j] = beta*P[i-2][j-2]-gama*P[i][j-2]

        if ((j+2)<Degree+2):
            for i in range(j+2,Degree+2):
                l = i-1
                m = j-1
                alpha = sqrt((2*l+1)*(l-m)*(l-m-1)/(2*l-3)/(l+m)/(l+m-1))
                if (m ==2):
                    beta = sqrt(2*(2*l+1)*(l+m-2)*(l+m-3)/(2*l-3)/(l+m)/(l+m-1))
                    gama = sqrt(2*(l-m+1)*(l-m+2)/(l+m)/(l+m-1))
                else:
                    beta = sqrt((2*l+1)*(l+m-2)*(l+m-3)/(2*l-3)/(l+m)/(l+m-1))
                    gama = sqrt((l-m+1)*(l-m+2)/(l+m)/(l+m-1))        
                P[i][j] = alpha*P[i-2][j]+beta*P[i-2][j-2]-gama*P[i][j-2]
    l = Degree
    m = Degree
    beta = sqrt((2*l+1)*(l+m-2)*(l+m-3)/(2*l-3)/(l+m)/(l+m-1))
    gama = sqrt((l-m+1)*(l-m+2)/(l+m)/(l+m-1))
    P[l+1][m+1] = beta*P[l+1-2][m+1-2]-gama*P[l+1][m+1-2]  
    return P

#用于看Pnm，不建议大规模使用，非常慢
def P_final(theta,n,m,Degree=180):
    Phi=pi/2-theta
    P=Pnm(Phi,Degree)
    return P[n+1][m+1]

值得注意的事情就是，θ（theta）是地心余纬，和φ（phi）的关系是θ=pi/2-φ；还有就是序号问题。因为代码是从matlab修改过来的，所以索引要加1。比如需要P21，在矩阵中的值是P[2+1,1+1]。这一点可以看到后面的函数P_final。Degree就是展开阶数，输入几就是几。