Quantenmechanik/ Magnetfeld/ Skript
Skript zum Zeeman-Effekt. Schreibt die Datei zeeman.svg
from math import sqrt
# computations for angular momentum
def colsalign(s) : # align string array s columns, strings start with ' '
n=len(s); m=32; cols=[0]*m; kmax=0
for i in range(n) :
t=s[i]; h=0; j=t.find(' ',h+1); k=0
while j>=0 :
if ((j-h)>cols[k]) : cols[k]=j-h-1
h=j; j=t.find(' ',h+1); k+=1
j=len(t)
if ((j-h)>cols[k]) : cols[k]=j-h-1
if (k+1)>kmax : kmax=k+1
# dbg print('Cols: '+str(cols[:kmax]))
for i in range(n) :
t=s[i]; h=0; j=t.find(' ',h+1); lt=len(t); k=0; v=''
while j<=lt :
u=t[h:j]+(' '*(cols[k]-(j-h-1))); v+=u
h=j; j= (lt+10) if (h==lt) else t.find(' ',h+1)
k+=1; j= lt if (j<0) else j
s[i]= v
return s
def hasfact(x,z) : return (x>0) and ((x % z)==0)
def havefact(x,y,z) : # x,y have common factor z?
fx=hasfact(x,z); fy=hasfact(y,z)
return (fx and fy) or (fx and (y==0)) or (fy and (x==0))
def normalize(x) : # rational-square root number x
pn=[2,3,5,7,11,13,17,19,23]; ln= len(pn) # small prime nbrs
ps=[0,1,4,9,16,25,36,49,64,81,100]; ls=len(ps) # small squares
pr,pi,q,sp,sq = tuple(x)
pr,sgr= (-pr,-1) if (pr<0) else (pr,1)
pi,sgi= (-pi,-1) if (pi<0) else (pi,1)
for i in range(2,len(ps)) :
z=ps[i]
while hasfact(sp,z) : sp=div(sp,z); pr *= i; pi *= i
while hasfact(sq,z) : sq=div(sq,z); q *=i
for i in range(ln) :
z=pn[i]
while hasfact(sp,z) and hasfact(sq,z) : sp=div(sp,z); sq=div(sq,z)
while havefact(pr,pi,z) and hasfact(q,z) :
pr=div(pr,z); pi=div(pi,z); q=div(q,z)
if (pr==0)and(pi==0) : q=1;sp=1;sq=1
x[0],x[1],x[2],x[3],x[4] = (pr*sgr,pi*sgi,q,sp,sq)
def real(x) : # rational square root number, get real part
pr,pi,q,sp,sq = tuple(x); return sqrt(sp/(sq+0.0))*pr/(q+0.0)
def encodenb(x) :
pr,pi,q,sp,sq = tuple(x); s=''; w=''
if (sp!=1)or(sq!=1) :
w='\\sqrt{';
w += str(sp) if(sq==1) else ('\\frac{'+str(sp)+'}{'+str(sq)+'}')
w +='}'
if pr != 0 :
u=str(pr); s += u if (q==1) else ('\\frac{'+u+'}{'+str(q)+'}')
s = w if((s=='1')and(w!='')) else (s if (s=='1') else (s+w))
if pi != 0 :
if (pr!=0) and (pi>0) : s +='+i'
elif (pi<0) : s += '-i'; pi=-pi
else : s += 'i'
u=str(pi); u='' if ((pi==1)and(q==1)) else u
s += u if (q==1) else ('\\frac{'+u+'}{'+str(q)+'}')
s = s+w
if (pr==0)and(pi==0) : s='0'
return s
def displaynb(x) : # simple readable, root in square brackets
pr,pi,q,sp,sq = tuple(x); s=''; w=''
if (sp!=1)or(sq!=1) :
w='[';
w += str(sp) if(sq==1) else (str(sp)+'/'+str(sq))
w +=']'
if pr != 0 :
u=str(pr); s += u if (q==1) else (u+'/'+str(q))
s = w if((s=='1')and(w!='')) else (s if (s=='1') else (s+w))
if pi != 0 :
if (pr!=0) and (pi>0) : s +='+i'
elif (pi<0) : s += '-i'; pi=-pi
else : s += 'i'
u=str(pi); u='' if ((pi==1)and(q==1)) else u
s += u if (q==1) else (u+'/'+str(q))
s = s+w
if (pr==0)and(pi==0) : s='0'
return s
def reformat(x) : # int, rational, rcomplex to root numbers
if type(x)==type(1) : z=[x,0,1,1,1]; n=-1
elif type(x)==type([]) : n=len(x)
if n==1 : z=[x[0],0,1,1,1]
if n==2 : z=[x[0],0,x[1],1,1]
if n==3 : z=[x[0],x[1],x[2],1,1]
if n==5 : z=x
if (n>5)or(n==0)or(n==4): print('reformat '+str(x)+' not supported.'); exit()
return z
def multip(x,y) : # format (u,v,w,p,q) = (u+iv)/w sqrt(p/q)
xx=reformat(x); yy=reformat(y)
xpr,xpi,xq,xsp,xsq= tuple(xx); ypr,ypi,yq,ysp,ysq= tuple(yy)
pr= xpr*ypr-xpi*ypi; pi= xpr*ypi+ypr*xpi; q=xq*yq; sp=xsp*ysp;sq=xsq*ysq
z= [pr,pi,q,sp,sq]; normalize(z); return z
def addit(x,y) :
xx=reformat(x); yy=reformat(y)
return lincomb(xx,yy,1,1)
def sqroot(x) : # input should be (p,0,q,1,1)
pr,pi,q,sp,sq = tuple(x)
if (pi!=0)or(sp!=1)or(sq!=1) : print('sqroot undefined.'); exit()
sp=pr;sq=q;pr=1;q=1; z=[pr,pi,q,sp,sq]; normalize(z); return z
def imagi(x) : pr,pi,q,sp,sq = tuple(x); return [-pi,pr,q,sp,sq]
def scale(x,f) : pr,pi,q,sp,sq = tuple(x); return [pr*f,pi*f,q,sp,sq]
def lincomb(x,y,f,g) : # f,g integer
# a sqrt(r/s) + b sqrt(p/q) ?= sqrt(rp/qs) ?
# can only be combined if r/s and p/q differ by square factor T/U ?
# first write rq/qs, ps/qs and test if rq,ps have square factor
ps=[0,1,4,9,16,25,36,49,64,81,100]; ls=len(ps) # small squares
xpr,xpi,xq,xsp,xsq= tuple(x); xpr *= f; xpi *=f
ypr,ypi,yq,ysp,ysq= tuple(y); ypr *= g; ypi *=g
if (xpr==0)and(xpi==0) : z=[ypr,ypi,yq,ysp,ysq]
elif (ypr==0)and(ypi==0) : z= [xpr,xpi,xq,xsp,xsq]
else :
sq=xsq*ysq; xp=xsp*ysq; yp=ysp*xsq
for i in range(2,len(ps)) :
z=ps[i]
while hasfact(xp,z) : xp=div(xp,z); xpr *= i; xpi *= i
while hasfact(yp,z) : yp=div(yp,z); ypr *= i; ypi *= i
# here xp/sq new sqrt arg? if yp=n xp, rewrite x ?
if xp != yp : print('lincomb:sqrt clash'+str([xp,yp])); exit()
sp=xp; q=xq*yq; pr=xpr*yq+ypr*xq; pi=xpi*yq+ypi*xq
z= [pr,pi,q,sp,sq]
normalize(z); return z
def matshow(s,x,new=0) :
# new= 1 simplified, columns. new=0 raw.
# new= 2 use encodenb, latex-like, check square hermitian
n=len(x); m=len(x[0]); print(s); t=''; ok=True; tt=['']*n
for i in range(n) :
for k in range(m) :
if (new==2) and (k>0) : t +=' &'
if (new==1) : t +=' '+displaynb(x[i][k])
if (new==2) : t +=' '+encodenb(x[i][k])
if new==0 : t +=' '+str(x[i][k]) # beware about v, we need a copy!
if new==2 :
u=x[i][k]; v= [] + x[k][i]; v[1]= -v[1]; ok= ok and isequal(u,v)
if (new==2) and (i<(n-1)) : t +=' \\\\';
tt[i]=t; t=''
if new==1 : colsalign(tt)
for i in range(n) : print(tt[i])
if new==2 : print('Test Hermitian: '+str(ok))
return tt
def vecshow(s,x) :
n=len(x); print(s); t=''
for k in range(n) : t +=' '+displaynb(x[k])
print(t)
def paulimats() :
z=[0,0,1,1,1]; u=[1,0,1,1,1]; i=[0,1,1,1,1]; j=[0,-1,1,1,1]; v=[-1,0,1,1,1]
sig=[[],[],[],[]]
sig[0]= [[u,z],[z,u]]
sig[1]= [[z,u],[u,z]] ; matshow('S1',sig[1])
sig[2]= [[z,j],[i,z]] ; matshow('S2',sig[2])
sig[3]= [[u,z],[z,v]] ; matshow('S3',sig[3])
for h in range(1,4) :
for k in range(1,4) :
a= matmul(sig[h],sig[k]); matshow('S'+str(h)+'*S'+str(k),a)
def div(p,q) :
if p>=0 : return int(p/q)
else : return -int(-p/q)
def getprimes(n) : # want n prime numbers
p=[2]; i=1; plast=2
while i<=n :
q=plast+1; found=False
while not found:
j=0
while ((q%p[j])!=0) and ((p[j]*p[j])<=q) : j +=1
if (q%p[j] != 0) : p += [q]; plast=q; found=True; i+=1
else : q+=1
return p
def greatestdiv(x) : # greatest divisor, skip 0s in list
n=len(x); y=[0]*n; w=[0]*n; h=0; ymin=-1; gcd=1
pn=[2,3,5,7,11,13,17,19,23,29,31,37,41,43,47]; ln= len(pn) # small primes
for i in range(n) :
z=abs(x[i])
if z>0 : ymin= z if(ymin<0) else min(ymin,z); w[h]=z; h+=1
w=w[:h]; ip=0
while (ip<ln)and(pn[ip]<=ymin) :
z=pn[ip]; ok=True
while ok :
for i in range(h) : ok= ok and ((w[i]%z)==0)
if ok :
for i in range(h) : w[i]= div(w[i],z)
ymin= div(ymin,z); gcd *= z
ip +=1
return gcd
def smallestmult(x) : # smallest multiple of nonzero numbers
n=len(x); y=abs(x[0])
for i in range(1,n) :
z=abs(x[i]); y = y if (z==0) else div(y*z,greatestdiv([y,z]))
return y
def divide(p,q) : # improve integer fractions
def hasfact(x,z) : return (x>0) and ((x % z)==0)
pn=[2,3,5,7,11,13,17,19,23,29,31,37]; ln= len(pn) # small prime nbrs
p,q,s = (abs(p),abs(q),-1) if ((p*q)<0) else (abs(p),abs(q),1)
if p==0 : q=1
for i in range(ln) :
z=pn[i]
while hasfact(p,z) and hasfact(q,z) : p=div(p,z); q=div(q,z)
return [s*p,q]
'''
Basic fine-struct & zeeman effect model.
Old basis: L3 and S3 diag. New: 2L*S=[j(j+1)-l(l+1)-s(s+1)] diag.
S3 does not commute with J^2, but with J3,L3,S3,S^2,L^2.
matrix elements of L*S, L3 and S3 in new basis Ci, from old basis Bk:
(Ci,Op,Cj)= (sum k Mik, sum h Mjh) (Bk,Op,Bh) with matrix transform M
matrices l3,s3 have same nonzero patterns
scalar products <j,j3(L3,S3)j',j3'> wont stay in j-block, but j3-blocks.
'''
def fullmatrix(diag,mat) : # diag + x * mat, latex-style
t='\\begin{pmatrix}{llllll}\n'; n=len(diag)
for i in range(n) :
for k in range(n) :
mik=mat[i][k]; zero= (mik[0]==0); dg=diag[i]
v=encodenb(mik); neg=(v[0]=='-'); zero=(v[0]=='0')
w= ('-x\cdot'+v[1:]) if neg else ('0' if zero else ('x\cdot'+v))
if i==k : u=encodenb(dg); u+= w if neg else ('+'+w)
else : u=w
if k>0 : t += ' & '
t+= u
if (i<(n-1)) : t +=' \\\\';
t+='\n'
t+='\\end{pmatrix}\n'
print(t)
def jj3(j,j3, l,l3, s,s3) : # get coeffs of vector (j,j3). numbers halfints
j2=div(j[0]*2,j[1]); j32=div(j3[0]*2,j3[1])
l2=div(l[0]*2,l[1]); l32=div(l3[0]*2,l3[1])
s2=div(s[0]*2,s[1]); s32=div(s3[0]*2,s3[1])
ok=True; jmax=abs(l2)+abs(s2); jmin=abs(l2-s2); z=[0,1]; sign=0
if (s2 !=1) : print('spin s '+str(s)+' not supported'); ok=False
if (j2>jmax)or(j2<jmin) : print('sum j '+str(j)+' out of range'); ok=False
if (j32>j2)or(j32<-j2)or(l32>l2)or(l32<-l2)or(s32>s2)or(s32<-s2) :
print('one of '+str([j3,l3,s3])+' out of range'); ok=False
if (j2>l2) : # contribute l3,s3= (j3-0.5,+) and (j3+0.5,-)
if (l32==(j32-1))and(s32==1) : z=[l2+j32+1,2*l2+2]; sign=1
if (l32==(j32+1))and(s32==-1) : z=[l2-j32+1,2*l2+2]; sign=1
else : # contribute l3,s3= (j3+0.5,-) and (j3-0.5,+)
if (l32==(j32+1))and(s32==-1) : z=[l2+j32+1,2*l2+2]; sign=1
if (l32==(j32-1))and(s32==1) : z=[l2-j32+1,2*l2+2]; sign=-1
if (not ok)or(z[0]==0) : z= [0,0,1,1,1]
else : z=[sign,0,1,z[0],z[1]]
return z
def testcleb(j,l) : # Clebsch-Gordan case spin 1/2 only
# make a matrix of coeffs for (j,j3) in (l,l3)(x)(s,s3)
noisy=False # example l=[1,1]; j=[3,2]; s fixed for spin 1/2
nl= div(2*l[0],l[1])+1
s=[1,2]; ns= div(2*s[0],s[1])+1
nj= div(2*j[0],j[1])+1; mx=[[]]*nj
for ij in range(nj) :
j3=j if (ij==0) else [j3[0]-j3[1],j3[1]]; hj=0; mx[ij]=[[]]*(nl*ns)
if ij==0 : s3val=[[]]*(nl*ns); l3val=[[]]*(nl*ns)
for il in range(nl) :
l3=l if (il==0) else [l3[0]-l3[1],l3[1]]
for ks in range(ns) :
s3=s if (ks==0) else [s3[0]-s3[1],s3[1]]
if ij==0 : s3val[hj]=s3; l3val[hj]=l3
z= jj3(j,j3, l,l3,s,s3); normalize(z); mx[ij][hj]=z; hj+=1
if noisy : print(str([ij,il,ks,z]))
# matshow('',mx,1); return mx,s3val,l3val
return mx,s3val,l3val
def matpermute(m,s) : # matrix line-column permutation
n=len(m); mp=[[]]*n; p=[0]*n
for i in range(n) : p[i]=int(s[i]); mp[i]=[[]]*n
for i in range(n) :
for k in range(n) : mp[i][k]=m[p[i]][p[k]]
return mp
def vecpermute(v,s) : # vector permutation
n=len(v); vp=[[]]*n
for i in range(n) : p=int(s[i]); vp[i]= v[p]
return vp
def spectralcurve(dg,m) : # vary Dg+xM and compute 6 values for each
# dg diagonal matrix, M 2x2-block matrix
# wrong, should get 3 not 5 asymptotics ? 2 triplets with diff g-factors?
nx=100; fx=5.0/nx; curve=[[]]*6; y=[0]*6
for j in range(6) : curve[j]=[0.0]*(2*nx+2)
for ix in range(nx+1) :
x=fx*ix; y[0]=real(dg[0])+real(m[0][0])*x
y[1]=real(dg[1])+real(m[1][1])*x # 1st block trivial
for j in range(2) :
k=2+2*j; a=real(dg[k])+real(m[k][k])*x;
b=real(m[k][k+1])*x; c=real(m[k+1][k])*x
d=real(dg[k+1])+real(m[k+1][k+1])*x
# solve (a-z)(d-z)=bc=zz-(a+d)z+ad; (z-p)^2=...
q=b*c-a*d; p=0.5*(a+d); root=sqrt(q+p*p)
y[k]=p+root; y[k+1]=p-root
for j in range(6) : curve[j][2*ix]=x; curve[j][2*ix+1]=y[j]
return curve
def zeeman() :
# use special multip and addit for square-root numbers
l=[1,1]; s=[1,2]; ja=[3,2]; jb=[1,2] # combine ang-momenta l,s to ja,jb
mx4,s3val,l3val= testcleb(ja,l)
mx2,s3val,l3val= testcleb(jb,l)
mx= mx4+mx2; n=len(mx); diag=[[]]*n; noisy=False # mx is 6x6 matrix
print('Dimension of mx='+str(len(mx))+','+str(len(mx[0])))
# diag vector has (LS) = [j(j+1)-l(l+1)-s(s+1)]/2 entries
z=addit(ja,1); ua=multip(ja,z)
z=addit(jb,1); ub=multip(jb,z)
z=addit(l,1); v=multip(l,z)
z=addit(s,1); w=multip(s,z)
x=lincomb(ua,v,1,-1); x=lincomb(x,w,1,-1); ua=multip(x,[1,2])
x=lincomb(ub,v,1,-1); x=lincomb(x,w,1,-1); ub=multip(x,[1,2])
for i in range(len(mx4)) : diag[i]= ua
for i in range(len(mx4),n) : diag[i]= ub
if noisy : print('ua,ub='+str(ua)+str(ub))
if noisy : matshow('Matrix mx',mx,1)
s3=[[]]*n; l3=[[]]*n; l2s=[[]]*n
for i in range(n) :
s3[i]=[[]]*n; l3[i]=[[]]*n; l2s[i]=[[]]*n
for j in range(n) : # compute op[i][j]
s3sum=0; l3sum=0
for k in range(n) :
u= s3val[k]; v=l3val[k] # the diagonal op(k,h)
u1= multip(u,mx[i][k]); v1=multip(v,mx[i][k])
for h in (k,) : # not 'in range(n)' because diagonal
u2=multip(u1,mx[j][h]); s3sum= addit(s3sum,u2)
v2=multip(v1,mx[j][h]); l3sum= addit(l3sum,v2)
s3[i][j]=s3sum; l3[i][j]=l3sum
l2s[i][j]=lincomb(l3sum,s3sum,1,2) # L+2S matrix
if noisy : matshow('Matrix s3',s3,1); matshow('Matrix l3',l3,1)
if noisy : matshow('Matrix l2s',l2s,1)
l2sp= matpermute(l2s,'031425'); diagp= vecpermute(diag,'031425')
matshow('Matrix l2s permuted',l2sp,1)
vecshow('Vector diag permuted', diagp)
data= spectralcurve(diagp, l2sp); multiplot(data) # SVG plot
print('\nPlot zeeman.svg written, for Perturbation Matrix\n')
fullmatrix(diagp, l2sp) # Latex output
# svg sketch 6 levels vs fieldstrength
def multiplot(d) :
n=len(d); xmin=1e12;xmax=-xmin;ymin=xmin;ymax=xmax
txt=[]; fn='zeeman'; ix=[0,1,2,3,4,5,6,7,8,9]
for i in range(n) :
k=0; kmax=len(d[i]); #print('Curve '+str(i)+' length '+str(kmax))
while k<kmax :
x=d[i][k]; y=d[i][k+1]; k+=2
xmin= x if(x<xmin) else xmin; xmax= x if (x>xmax) else xmax
ymin= y if(y<ymin) else ymin; ymax= y if (y>ymax) else ymax
#print('xmin,ymin, xmax,ymax'+str([xmin,ymin, xmax,ymax]))
xyf=[ xmin,ymin, xmax,ymin, xmax,ymax, xmin,ymax, xmin,ymin]
xys=[xyf]+d
xs,fx,xpix= xmin, 750/(xmax-xmin),25 # scale user units, min in pixels
ys,fy,ypix= ymin,-550/(ymax-ymin),575 ; scale=(xs,fx,xpix, ys,fy,ypix)
colo=['#aaa','#00a','#077','#0a0','#770','#a00','#707','#0ff','#f0f']
sd=svgdump(); buf=sd.plotbegin(800,600,1); i=0; dy=(ymax-ymin)/20.0
for j in range(len(txt)) :
buf += label(sd,scale, txt[j], 0.2*xmin+0.8*xmax,ymax-dy-j*dy)
for t in xys : buf += curve(sd,scale, t, fill=colo[ix[i]]); i+=1
g=open(fn+'.svg','w'); g.write(buf+sd.plotend()); g.close()
class svgdump :
def plotbegin(sd, w,h,f) : # args: width height zoomfactor
W=str(w); H=str(h); u=''; fw=str(f*w); fh=str(f*h)
s=( '<?xml version="1.0" encoding="UTF-8" standalone="no"?>\n'+
'<!DOCTYPE svg>\n'+
'<svg xmlns="http://www.w3.org/2000/svg"\n'+
' xmlns:xlink="http://www.w3.org/1999/xlink"\n'+
' width="'+W+'" height="'+H+'" viewBox="0 0 '+fw+' '+fh+'">\n')
return s
def plotend(sd) : return '</svg>\n'
def move(sd,x,y) : return 'M'+str(x)+','+str(y)
def line(sd,x,y) : return 'L'+str(x)+','+str(y)
def labl(sd,s,x,y, size=14,fill='#000') :
f=' font-family="Arial"'
t='<tspan x="'+str(x)+'" y="'+str(y)+'">'+s+'</tspan>'; fsz= str(size)
return ('<text fill="'+fill+'"'+f+' font-size="'+fsz+'px">'+t+'</text>\n')
def path(sd,w,s,lc='#000') : # path w=with lc=linecolor
t='<path fill="none" stroke="'+lc+'" stroke-width="'+str(w)+'" '
t+= ' stroke-linecap="round" '
return t + 'd="\n'+s+'" />\n'
# end class svgdump
def aprx(x) : # 3-digit approximation
y= (-x) if(x<0) else x; z= int(1e3*y+1e-3); z= (-z) if (x<0) else z
return z/1000.0
def curve(sd,scale,xy, fill='#000', width=1) :
(xs,fx,xmin, ys,fy,ymin) = scale; s=''; n= div(len(xy),2)
for i in range(n) :
x= aprx(xmin+fx*(xy[2*i]-xs)); y= aprx(ymin+fy*(xy[2*i+1]-ys))
s+= sd.move(x,y) if(i==0) else sd.line(x,y)
return sd.path(width,s,lc=fill)
def label(sd,scale,txt,tx,ty) : # tx,ty plot coordinates not pixels
(xs,fx,xmin, ys,fy,ymin) = scale
x= aprx(xmin+fx*(tx-xs)); y= aprx(ymin+fy*(ty-ys))
return sd.labl(txt,x,y)
zeeman()