MSE12 版 (精华区)

发信人: chumsdock (微笑服务), 信区: MSE12       
标  题: 退磁化。。
发信站: BBS 听涛站 (Thu Oct 20 01:05:30 2005), 转信

from span import *

#Depolarzition Matrix
def demag(pos):
    N=mat('0.,0.,0.;0.,0.,0.;0.,0.,0.')
    a,b,c=1.,1.,1.
    x,y,z=pos
    #Initialization
    for q in (-1.,1.):
        for w in (-1.,1.):
            R1=x
            R2=b/2.+q*y
            R3=c/2.+w*z
            R=sqrt(R1*R1+R2*R2+R3*R3)
            N[0,0]=N[0,0]+arctan(R2*R3/R1/R)
            N[1,0]=N[1,0]+q*log((R-R3)/(R+R3))
            N[2,0]=N[2,0]+w*log((R-R2)/(R+R2))
    #In order to preserve Precision, calcualte 1/4pi and 1/8pi last.
    N[0,0]=-N[0,0]/4./pi
    N[1,0]=-N[1,0]/8./pi
    N[2,0]=-N[2,0]/8./pi
    return(N)

def rotate(theta,phi):
    R = mat([[cos(theta)*cos(phi),cos(theta)*sin(phi),-sin(theta)],
             [-sin(phi),cos(phi),0],
             [sin(theta)*cos(phi),sin(theta)*sin(phi),cos(theta)]])
    return R
def rotv(pos,phi):
    x,y,z=pos
    vec=mat([[x],[y],[z]])
    vec=rotate(0,phi)*vec
    return vec.T.array[0]

def demag6(pos):
    N6=mat('0.,0.,0.;0.,0.,0.;0.,0.,0.')
    pos1=pos
    for i in range(6):
        pos1=rotv(pos,i*pi/3);pos1[0]-=sqrt(3)/2
        mag=rotate(0,i*pi/3).T*demag(pos1)
        N6=N6+mag
    return N6

--


               不主动,勤拒绝,不负全责。


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

此主题相关图片如下:b.gif (26809 字节)
[百宝箱] [返回首页] [上级目录] [根目录] [返回顶部] [刷新] [返回]
Powered by KBS BBS 2.0 (http://dev.kcn.cn)
页面执行时间:0.780毫秒