发信人: 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 字节)