MSE12 版 (精华区)

发信人: chumsdock (微笑服务), 信区: MSE12       
标  题: 韦jj编程4.7的python实现
发信站: BBS 听涛站 (Tue Oct 18 21:26:08 2005), 转信

from scipy import *
from Numeric import *
#--------Depolarzition Matrix--------------------------------------------
def polar(pos):
    #-----Initialization-------------------------------------------------
    N=mat("""0., 0., 0.;
             0., 0., 0.;
             0., 0., 0.""")
    a,b,c=1.,1.,1.
    x,y,z=pos
    #-----Initialization end---------------------------------------------
    for p in (-1.,1.):
        for q in (-1.,1.):
            for w in (-1.,1.):
                R1=a/2.+p*x
                R2=b/2.+q*y
                R3=c/2.+w*z
                R=sqrt(R1*R1+R2*R2+R3*R3)
                N[0,0]+=arctan(R2*R3/R1/R)
                N[1,1]+=arctan(R3*R1/R2/R)
                N[2,2]+=arctan(R1*R2/R3/R)
                N[1,0]+=q*p*log((R-R3)/(R+R3))
                N[2,0]+=p*w*log((R-R2)/(R+R2))
                N[2,1]+=q*w*log((R-R1)/(R+R1))
    #-----To preserve Precision, calcualte 1/4pi and 1/8pi last.--------
    N[0,0]=N[0,0]/4./pi
    N[1,1]=N[1,1]/4./pi
    N[2,2]=N[2,2]/4./pi
    N[0,1]=N[1,0]=N[1,0]/8./pi
    N[0,2]=N[2,0]=N[2,0]/8./pi
    N[1,2]=N[2,1]=N[2,1]/8./pi
    return(N)
#-------------Main------------------------------------------------------
if __name__ == '__main__':
    x0, step, xend= -4, 0.2, 3.
    print trace(polar([0,0,0])),trace(polar([1,1,1]))
    while x0<xend:
        x1=0.5*exp(x0)
        pos1=[x1,0,0]
        pos2=[-x1,0,0]
        E1=-polar(pos1)*mat('1;0;0')
        E2=-polar(pos2)*mat('1;0;0')
        print x0,x1,E1[0,0],-x1,E2[0,0]
        x0=x0+step

本程序的工作量很小，因为矩阵的运算都用Scipy包实现了，polar()函数是求一个
坐标的退极化矩阵，输入参量本身就是向量(当然是通过python的列表实现的)
注意矩阵的元素是从0,0开始的，N11对角元就是N[0,0].

trace那一步是为了验证57页4.16式，取了体积内外的两点进行计算，注意
边界是+/-0.5而不是+/-1.

x轴坐标暂时取了对数形式，也可以取倒数形式，也可以随便取，用x0变量循环控制

用python的主要原因是它的语法简洁，程序直观，编写时间短，
虽然程序效率是很低的(解释型)，但是对付作业的要求还是没有问题，
如果进行比较大规模的计算，还是推荐FORTRAN。


--
    OK 5.4


※ 来源:·BBS 听涛站 tingtao.net·[FROM: 59.66.203.181]