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()