球协函数图像的绘制

本文利用 Python 编程语言实现归一化球协函数图像的数据准备,并用 Mathematica 绘制图像。

归一化球协函数为:

 Y_{lm}(\theta,\phi)=\sqrt{\frac{2l+1}{4\pi}\frac{(l-m)!}{(l+m)!}}P_l^m(\cos(\theta))e^{im\phi}

其中 l=0,1,2,\cdots,m=0,\pm 1,\pm 2, \cdots, \pm l,\quad P_l^m(\cos \theta) 为连带勒让德函数

m>0 时,把 m 递归到零:

\sqrt{1-x^2}P_l^m(x)=\frac{1}{2l+1}[(l-m+1)(l-m+2)P_{l+1}^{m-1}(x)-(l+m-1)(l+m)P_{l-1}^{m-1}(x)]

m<0 时,把 m 递归的零:

\sqrt{1-x^2}P_l^m(x)=\frac{1}{2l+1}[P_{l-1}^{m+1}(x)-P_{l+1}^{m+1}(x)]

m=0 时,把 l 递归到零:

(l-m)P_l^m(x)=(2l-1)xP_{l-1}^m(x)-(l+m-1)P_{l-2}^m(x)

利用以上三式就可以通过递归程序定义连带勒让德函数。

下面是 Python 程序代码

#!/usr/bin/python
from math import sin,cos,pi

#递归定义连带勒让德函数
def pfun(m,l,x):
    if m==0:
        if l==0:
            return 1
        elif l==1:
            return x
        elif l==2:
            return 0.5*(3*x**2-1)
        else:
            return ((2*l-1)*x*pfun(m,l-1,x)-(l+m-1)*pfun(m,l-2,x))/(l-m)
    elif m>0:
        return ((l-m+1)*(l-m+2)*pfun(m-1,l+1,x)-(l+m-1)*(l+m)*pfun(m-1,l-1,x))/((2*l+1)*((1-x**2)**(1./2.)))
    elif m<0:
        return (pfun(m+1,l-1,x)-pfun(m+1,l+1,x))/((2*l+1)*((1-x**2)**(1./2.)))
#阶乘
def fac(n):
    s=1
    for i in range(1,n+1):
        s*=i
    return s

#球协函数实部
def yfunre(m,l,q,f):
    return ((((2*l+1)*fac(l-m)/(4*pi*fac(l+m)))**(1./2.))*pfun(m,l,cos(q)))*cos(m*f)

#球协函数虚部
def yfunim(m,l,q,f):
    return ((((2*l+1)*fac(l-m)/(4*pi*fac(l+m)))**(1./2.))*pfun(m,l,cos(q)))*sin(m*f)

#数据计算
l=int(raw_input("l= "))                   #输入数据 l:0,1,2,...
m=int(raw_input("m= "))                   #输入数据 m:-l,...,-2,-1,0,1,2,...,l
step=float(raw_input("step= "))           #输入迭代步长: 0.03 到 0.05 较合适

print "computing the data ..."

fh=step
qh=step
data=[]
q=qh
while True:
    f=fh
    while True:
        #r=yfunim(m,l,q,f)
        #r=yfunre(m,l,q,f)
        r=(yfunre(m,l,q,f)**2+yfunim(m,l,q,f)**2)**(1./2.)
        data.append([r*sin(q)*cos(f),r*sin(q)*sin(f),r*cos(q)])         #注意球坐标系到笛卡尔坐标系的转换
        f+=fh
        if f>2*pi:break
    q+=qh
    if q>=pi:break

#数据储存到文件
fil=raw_input("file name is: ")                   #输入文件名:任意
print "writing the data ..."
out=open(fil,"w")
out.write("{\n")                                  #按照 Mathematica 列表格式储存
for item in data[:-1]:
    out.write("{")
    out.write("%f"%item[0]+",")
    out.write("%f"%item[1]+",")
    out.write("%f"%item[2]+"}")
    out.write(",\n")

out.write("{"+"%f"%data[-1][0]+",")
out.write("%f"%data[-1][1]+",")
out.write("%f"%data[-1][2]+"}")
out.write("}")

out.close()
print "done"

Python 代码下载: pycode.py [1.67 kB] -

如果数据储存在文件:“data.txt“,利用下面两行 Mathematica 代码即可绘制出图像。

data=ToExpression[Import["data.txt"]];
ListPointPlot3D[data, PlotRange->All,BoxRatios->Automatic,ImageSize->400]

下面是 m=0 时, l=2l=3 的球协函数三维图像:

l2l3

发表评论

电子邮件地址不会被公开。 必填项已用*标注

*