Created
May 5, 2013 04:15
-
-
Save ZhenghengLi/5519671 to your computer and use it in GitHub Desktop.
qixiehanshu
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #!/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" |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment