Quantenmechanik/ Drehimpuls-Kopplung/ CG-Skript
Skript zu Clebsch-Gordan-Koeffizienten
[Bearbeiten]Das Skript listet Clebsch-Gordan-Koeffizienten zu Paaren von Drehimpulsen L und S. Sie sollen zum Gesamt-Drehimpuls mit dem Operator J=L+S kombiniert werden. Nur die nichtverschwindenden Werte für die möglichen Eigenvektoren des Gesamt-Drehimpulses werden aufgezählt.
Alle Rechnerei ist Ganzzahl-Arithmetik, mit der "Rationale Wurzelzahlen"
aufgebaut werden. Das sind hier Quintupel von Ganzzahlen [r,p,q,u,v] mit
q>0, v>0 und der Bedeutung: .
Der Imaginärteil p wird
in diesem Skript nicht benutzt, er ist immer Null.
Die Funktion reformat() verwandelt ganze und einfach-rationale, d.h. Paare
[p,q], in solche Quintupel. Rechenschritte erfolgen mit den Funktionen
multip,addit,sqroot,recip,scale,lincomb,quotient,iszero.
Die Addition geht nur, wenn beide Argument rationale Vielfache derselben Wurzel
sind, sonst bricht das Programm ab.
Zu Eingabe-Paaren (L,S) berechnet das Skript alle Clebsch-Gordan-Koeffizienten.
L und S sind Vielfache von (1/2). Also (2L)=1,2,3,... und (2S)=1,2,3,...
Eine Clebsch-Gordan-Matrix enthält (2L+1)*(2S+1) zum Quadrat rationale
Wurzelzahlen.
De Matrix ist hier eine Tabelle mit 4 Indizes .
Quantenzahl J ist die absteigende Folge (L+S),(L+S-1)...|L-S|.
J,L und S sind entweder halb- oder ganzzahlig.
i,j,k,l fangen bei Null an, meinen aber folgende physikalischen Werte:
- Index i adressiert den Wert J, vom Maximum (L+S) abwärts.
- Index j adressiert den Wert J_3, in absteigender Reihenfolge J,J-1,...,-J.
- Index k adressiert den Wert L_3, Reihenfolge abfallend L,L-1,...,-L
- Index l adressiert S_3, Reihenfolge S,S-1,...,-S.
In physikalischer Notation ist .
Die Elemente mit festem und variabel nennen wir
etwas willkürlich eine Zeile der Matrix, bei umgekehrter Variation eine Spalte.
Für L=3,S=2 hat eine Zeile die Länge (2*3+1)(2*2+1)=35 und eine Spalte
(2*1+1)+(2*2+1)+(2*3+1)+(2*4+1)+(2*5+1)=3+5+7+9+11=35, ganz wie geplant.
Der Algorithmus folgt sklavisch dem Standard-Rezept. Angefangen wird mit
dem maximalen J,J_3 und einem Koeffizienten 1. Mit den Absteige-Operatoren werden
die Zeilenvektoren für die Folge der J_3 errechnet. Beim nächst niedrigeren
J wird angefangen mit dem orthogonalen Komplement zu den bis dahin bereits
vorhandenen Vektoren zur Quantenzahl J_3=J. Der Orthogonalisierungs-Algorithmus
bekommt n-tupel von Wurzelzahlen als Eingabe zu sehen, jeweils (n-1) davon.
Ist eine Menge von schon orthogonalen Vektoren und y ein Kandidat für
ihr Komplement, dann werden einfach alle seine Projektionen abgezogen:
.
Für werden brutal die
Standard-Einheitsvektoren durchprobiert bis ein orthonormierbares da ist.
Es gibt mindestens eine der (n) Achsen, die von den angelieferten (n-1) Vektoren
linear unabhängig ist, daher wird ein Komplement-Vektor gefunden.
Konvention: Der Eintrag zum maximalen Wert von J_3 und L_3 bei jeder
neuen Treppenstufe von J wird positiv gemacht.
Halbzahlige Variablen werden vermieden zugunsten ihres Doppels. Wenn ein Variablenname wie j2 oder m2 mit der Ziffer 2 endet, enthält er zweimal die physikalische Größe j oder m.
Das Skript wurde durchprobiert mit den Quantenzahlen L und S = (1/2),1,(3/2),2,(5/2),3. Keine Garantie darüber hinaus. Mit L=S=3 hat die CG-Matrix 2401=(7*7)*(7*7) Elemente, davon 216 nicht Null!
Der Code hat (langsame) Testsequenzen dazu, dass alle Zeilen und Spalten der Clebsch-Gordan-Matrix orthonormal sind, dass die Symmetrie beim Umklappen von stimmt und dass die Kopfzeile jedes j-Blocks auch mit den Racah-Formeln genau gleich herauskommt.
Definierende Gleichung der CG-Matrix:
Symmetrien der Koeffizienten:
Auf- und Absteigeregeln für festes J:
Beziehung zwischen J-Nachbarn bei festem :
Racah-Formel für den Spezialfall
Es existiert eine alleswissende, abschreckend komplizierte Racah-Formel für die Gesamtheit der CG-Koeffizienten.
Wigner-3j-Symbole sind eine alternative Form der CG-Koeffizienten.
- .
Kostprobe der Ausgabe
[Bearbeiten]Clebsch-Gordan L=(1/2) S=(1/2)
Clebsch-Gordan L=(1/2) S=1
Clebsch-Gordan L=(1/2) S=(3/2)
Clebsch-Gordan L=(1/2) S=2
Clebsch-Gordan L=1 S=(1/2)
Clebsch-Gordan L=1 S=1
Python-Listing
[Bearbeiten]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,121,144]; 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 finetune(x) : pn=[2,3,5,7,11,13,17,19,23]; ln= len(pn) # small primes xx=reformat(x); pr,pi,q,sp,sq = tuple(xx) nmax=max(sp,sq);imax=0 # minimize values under root for i in range(2,5) : tp=sp*i*i; tq=sq for k in range(ln) : z=pn[k] while hasfact(tp,z) and hasfact(tq,z) : tp=div(tp,z); tq=div(tq,z) if max(tp,tq)<=nmax: nmax=max(tp,tq);imax=i; pmax=tp;qmax=tq if imax>0 : q *=imax; sp=pmax; sq=qmax return [pr,pi,q,sp,sq] def displaynb(x) : # x is 5-tuple. to readable string, root in square brackets # bug: ambiguous output rt(2/3)=2*rt(1/6)= (1/3)*sqrt(6) # minimize max(p,q) in sqrt(p/q)? reduce {4,9,1/4,1/9)*p/q= p'/q'? # prefer sqrt(p) over sqrt(1/q) xx=reformat(x); pr,pi,q,sp,sq = tuple(xx); 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) xx=reformat(x); pr,pi,q,sp,sq = tuple(xx) 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 recip(x) : xx=reformat(x); pr,pi,q,sp,sq = tuple(xx); qq=pr*pr+pi*pi; y= [pr*q,-pi*q,qq,sq,sp]; normalize(y); return y def imagi(x) : pr,pi,q,sp,sq = tuple(x); return [-pi,pr,q,sp,sq] def scale(x,f) : xx=reformat(x); pr,pi,q,sp,sq = tuple(xx); return [pr*f,pi*f,q,sp,sq] def lincomb(x,y,f,g) : # f,g integer, x,y fractional root numbers # 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,121,144]; 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 # test common sqroot-denom 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,x,y])); 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 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] 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('j='+str(j)+' outside |l-s|,(l+s)'); 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,s=[1,2]) : # Clebsch-Gordan case spin 1/2 only # make a matrix of coeffs for (j,j3) in (l,l3)(x)(s,s3) # return matrix mx, rows in order j3max...j3min # mx rows ordered: (lmax,smax)...(lmax..smin) (lmax-1,smax)... noisy=2 # example l=[1,1]; j=[3,2]; s fixed for spin 1/2 nl= div(2*l[0],l[1])+1 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==1 : print(str([ij,il,ks,z])) if noisy==2 : matshow('',mx,1) return mx,s3val,l3val def rootclass(x) : # x is list of real 5-tuple numbers, are they additive? noisy=False; # maybe root-incompatible subsets cancel? if noisy : print(' rootclass'+str(x)) n=len(x); q=[0]*n; p=[0]*n; r=[0]*n; clas=[0]*n; pclas=[0]*n; ncl=0 for i in range(n) : p[i]=x[i][3];q[i]=x[i][4];r[i]=x[i][0] qq= smallestmult(q); ls=12 for i in range(n) : p[i]= div(qq*p[i],q[i]); j=0; found=False for k in range(2,ls) : # take squares out of root z= k*k while hasfact(p[i],z) : p[i]=div(p[i],z); r[i] *= k while (not found) and (j<ncl) : found=(p[j]==pclas[j]); j+=1 if not found : pclas[j]=p[i]; clas[i]=j; ncl +=1 else : clas[i]= j-1 # equivalence class of same p factor if noisy : print('p-classes: '+str(pclas[:ncl])); sp=[0]*ncl; sq=[1]*ncl for k in range(ncl) : for i in range(n) : if clas[i]==k : sp[k]= r[i]*sq[k]+ sp[k]*x[i][2]; sq[k]=sq[k]*x[i][2] if noisy : print('p-sums: '+str(sp)) # only one of them non-zero! for i in range(n) : x[i]=[r[i],0,x[i][2],p[i],qq] return clas,pclas[:ncl] def scalprod(x,y) : # square-root compatible subset n=len(x); m=len(y); w=[0]*n; z=0; bigbug=False if bigbug : print('x= '+str(x)+'\ny= '+str(y)) if n != m : print('scalprod length mismatch.') for i in range(n) : w[i]=multip(x[i],y[i]); z=addit(z,w[i]) if bigbug : print('w= '+str(w)) return z def quotient(x,y) : z=recip(y); return multip(z,x) def iszero(x) : if type(x)==type(1) : return (x==0) else : return (x[0]==0)and(x[1]==0) def orthogonize(y,x) : # make vector y orthonormal to vectorlist x[i] # as startvector y try maximal n linear indep ones, y=delta(i,k) n=len(y);m=len(x); ok=True;sok=True; z=[0]*n; norm=[[]]*m;pxy=[[]]*m for i in range(m) : ok= ok and (len(x[i])==n) if not ok : print('Vector Length Error') for i in range(m) : for j in range(i,m) : sp=scalprod(x[i],x[j]) # root-non-compatible here? if i==j : norm[i]=sp; pxy[i]=scalprod(y,x[i]) if i != j : sok= sok and iszero(sp) if not sok : print('Orthogonality error') for k in range(n) : z[k]= reformat(y[k]) for i in range(m) : u= quotient(pxy[i],norm[i]) for k in range(n) : v=multip(u,x[i][k]); z[k]=lincomb(z[k],v,1,-1) nz= scalprod(z,z); ok=True if iszero(nz) : print('Orthogonize result zero.'); ok=False else : for i in range(m) : sp=scalprod(x[i],z); ok= ok and iszero(sp) if not ok : print('Result z not orthogonal.') nz= sqroot(nz); qz= recip(nz) for k in range(n) : z[k] = multip(z[k],qz) sp=scalprod(z,z); sp=lincomb(sp,reformat(1),1,-1); nok= iszero(sp) if not nok : print('Result z not normalized.') return z,ok def signflip(z) : # make first entry positive n=len(z); flip= z[0][0]<0; nf= n if flip else 0 for i in range(nf) : z[i][0]= -z[i][0] if flip : print('Signflip') def orthoplement(x,noisy) : # orthogonal complement # if all of x root-compatible, forget the roots # solve linear equ on rationals, then normalize n=len(x); m=len(x[0]); ok=True if noisy : print(' orthoplem '+str(n)+' vectors of length '+str(m)) if (n+1)!=m : print('Bad number of of vectors!') if m==2 : return [scale(x[0][1],1),scale(x[0][0],-1)] # other sign? for i in range(n) : clas,pclas=rootclass(x[i]); ok= ok and (len(pclas)==1) if not ok : print('Roots not compatible.') if noisy : print('Orthoplement '+str(n)) n2=n; n2=0 #dbg, n2=n is fatally bad for k in range(n2) : print(str(x[k])) for j in range(len(x[k])) : if len(x[k][j])>4 : x[k][j][3]=1;x[k][j][4]=1 # roots away is fatal y=[0]*m; found=False; i=0 while (not found) and (i<m) : # try all unit vector startpoints y[i]=1; z,found=orthogonize(y,x); y[i]=0; i+=1 if not found : print('Cannot orthoplement m='+str(m)+', exit.'); exit() else : signflip(z) return z # faster? guess int-only startvector y; Fi=prod(k,k!=i,(xk,xk))(xi,y); # G=prod(k,(xk,xk)); y'= Gy- sum(i,Fi xi) . if zero, modify y. int-solution? def initclebsch(ja2,jb2,noisy) : # ja,jb half-integral, ja2,jb2 their double # matrices quad-index: [n][d][ja2+1][jb2+1], d=(da+db-x)...|da-db|(2) max2=ja2+jb2; min2=abs(ja2-jb2); n=div(max2-min2,2)+1; da=ja2+1;db=jb2+1 cg=[[]]*n; sd=0; d=max2+1; s='CG-Matrices Cols='+str(da*db)+' Rows=' for i in range(n) : cg[i]=[[]]*d ; s+=' '+str(d); sd+=d for j in range(d) : cg[i][j]=[[]]*da for k in range(da) : cg[i][j][k]=[0]*db # defines cg[i,j,k,h] d -= 2 j2=ja2+jb2; i=0; cg[i][0][0][0]=1 if noisy : print(s+' Sd='+str(sd)) return cg # quadruple-index array: j,m, ma,mb def gathervectors(cgm,jmax2,m2, ja2,jb2, noisy) : # cgm vectors filled at m2? # extract relevant row vectors from da*db-size submatrices # put together the row vectors for j2.3==m2 # to start (j2=m2), find an orthogonal vector to those found here n=len(cgm); mydim=m2+1; k=0; x=[[]]*n # all cgm with higher dims contribute terms=[0]*n ; nj=[0]*n; j2=m2; h=0; jm2=jmax2 for i in range(n) : cg=cgm[i]; cdim=len(cg); okdim= cdim>mydim # if ok, add a vector x[k] if (noisy and okdim) : print('gather from submatrix dim='+str(cdim)) # extract the da*db matrix entries to m2 # cg submatrix i, row that counts belongs to m2: ma2=jmax2;j=0; y=[0]*(jmax2+1) # items in cg start with max=ma2, min=mb2 while ma2>ja2 : ma2 -=1 k= div(jm2-m2,2) # relevant row in matrix cg while okdim and (ma2>=(-jmax2)) : # scan diagonal ma2+mb2=m2 mb2=m2-ma2; ok= (abs(ma2)<=ja2)and(abs(mb2)<=jb2); z=0 if ok : ia=div(ja2-ma2,2);ib=div(jb2-mb2,2); z=cg[k][ia][ib]; terms[i]+=1 if noisy : print(' ma2,mb2,ok='+str([ma2,mb2,ok])) if ok : y[j]=z; j+=1; dn= 1 if (not iszero(z)) else 0; nj[i]+=dn ma2 -=2 if (terms[i]>0) and (nj[i]>0) : x[h]=y[:j]; h+=1 jm2 -= 2 if noisy : print('Gather Terms: '+str(terms)+' NonZero: '+str(nj)) return x[:h] def downsweep(cgm, z, j2, ja2,jb2,noisy) : # one diagonal z is ready: for m2=j2, all ma,mb with ma+mb=m2 are defined. # m index ordering is top down. input vector z = diagonal m2=j2. nz=len(z); da=ja2+1; db=jb2+1; ix= div(ja2+jb2-j2,2); cg=cgm[ix]; d=j2+1; k=0 if noisy : print('downsweep ix,nz,j2,ja2,jb2='+str([ix,nz,j2,ja2,jb2])) if noisy : print(' input vector='+str(z)); jmax2=ja2+jb2 ma2=jmax2 # items in z ordered to start with maximal ma2, minimal mb2? while ma2>ja2 : ma2 -=1 # j2=2 ja2=jb2=1, must start ma2=1 not 2! while ma2>=(-jmax2) : # insert start-diagonal from vector z to cg[0] row mb2=j2-ma2; ok= (abs(ma2)<=ja2)and(abs(mb2)<=jb2); i=0;ia=-1;ib=-1 if ok : ia=div(ja2-ma2,2); ib=div(jb2-mb2,2); cg[i][ia][ib]=z[k]; k+=1 if noisy: print(' ins ma2,mb2, i,ia,ib,z='+str([ma2,mb2,i,ia,ib,z[k-1]])) ma2 -=2 m2=j2; i=0; nput=0; imax=j2 # reach source row cg[2j] while m2>(-jmax2) : # fetch data from index m2 ma2=jmax2 if noisy : print(' Pass i,j2,m2,ja2='+str([i,j2,m2,ja2])) while ma2>ja2 : ma2 -=1 while ma2>=(-jmax2) : mb2=m2-ma2; ok= (abs(ma2)<=ja2)and(abs(mb2)<=jb2)and(i<imax) if ok : # add stuff to cg[ia+1][ib] and cg[ia][ib+1] if in range ia=div(ja2-ma2,2); ib=div(jb2-mb2,2); x=cg[i][ia][ib] q=j2*(j2+2)-m2*(m2-2) # divisor for norm if ok and ((ia+1)<da) : if noisy : print(' put1 i,ia,ib'+str([i,ia,ib, i+1,ia+1,ib,x])) y=ja2*(ja2+2)-ma2*(ma2-2); y=[y,q]; y= sqroot(y); nput+=1 y=multip(y,x); z=cg[i+1][ia+1][ib]; cg[i+1][ia+1][ib]=addit(z,y) if ok and ((ib+1)<db) : if noisy : print(' put2 i,ia,ib'+str([i,ia,ib, i+1,ia,ib+1,x])) y=jb2*(jb2+2)-mb2*(mb2-2); y=[y,q]; y= sqroot(y); nput+=1 y=multip(y,x); z=cg[i+1][ia][ib+1]; cg[i+1][ia][ib+1]=addit(z,y) ma2 -= 2 m2 -=2; i+=1 if noisy : print('Downsweep done, nput='+str(nput)) def cgshow(j2,cg) : nj=len(cg); na=len(cg[0]); nb=len(cg[0][0]); ja2=na-1; jb2=nb-1; m2=j2 print('CG-Matrix Nj,j2,ja2,jb2='+str([nj,j2,ja2,jb2])); s='' for i in range(nj) : cgi=cg[i]; na=len(cgi); nb=len(cgi[0]) print('MatrixRow '+str(i)+' J2,M2='+str(j2)+','+str(m2)) for j in range(na) : for k in range(nb) : z=finetune(cgi[j][k]); s+=' '+displaynb(z) print(s); s=''; m2 -=2 def clebschgordan(ja2,jb2,noisy=True) : # ja2,jb2 are double of j-quanta. # ja=(ja2/2) may be called 'L' and jb=(jb2/2) called 'S' cgm= initclebsch(ja2,jb2,noisy) j2=ja2+jb2;k2=j2; jmax2=j2; jmin2=abs(ja2-jb2); i=0; z=[1] while j2>= jmin2 : if i>0 : x=gathervectors(cgm,jmax2,j2,ja2,jb2,noisy) if i>0 : z=orthoplement(x,noisy) downsweep(cgm, z, j2,ja2,jb2, noisy); j2-=2; i+=1 if noisy : for i in range(len(cgm)) : cgshow(k2,cgm[i]); k2-=2 return cgm def wikidump(cg) : # list nonzero entries in latex style def wnum(n2) : s='' if (n2>=0) else '-'; na=abs(n2); whole= (na%2)==0; n=div(na,2) if whole : s+=str(n) else : s+='\\tfrac{'+str(na)+'}2' return s def label(l2,l32,s2,s32,j2,j32,x) : s='\\langle'+wnum(l2)+','+wnum(s2)+';'+wnum(l32)+','+wnum(s32)+'|' s+=wnum(j2)+','+wnum(j32)+'\\rangle='+encodenb(reformat(x))+'\\\\\n' return s nj=len(cg); na=len(cg[0][0]); nb=len(cg[0][0][0]); s='\\begin{array}{l}\n' for ij in range(nj) : mj=len(cg[ij]) for im in range(mj) : for ia in range(na) : for ib in range(nb) : x=cg[ij][im][ia][ib]; ok= not iszero(x) j2=mj-1; l2=na-1; s2=nb-1; j32=j2-2*im; l32=l2-2*ia; s32=s2-2*ib if ok : s+= label(l2,l32,s2,s32,j2,j32,x) s+='\\end{array}\n'; return s def bugcheck(cgm,mxa,mxb) : def subcheck(mx,my) : na=len(mx); nb=len(my[0]); nc=len(my[0][0]); bugs=0 for i in range(na) : h=0 for j in range(nb) : for k in range(nc) : x=reformat(mx[i][h]); y=reformat(my[i][j][k]); z=lincomb(x,y,1,-1) h+=1; bugs+= (0 if iszero(z) else 1) return bugs n=len(cgm); ok=(n==2) na=len(cgm[0]); nb=len(cgm[1]); ma=cgm[0]; mb=cgm[1] ok= ok and (na==len(mxa)) and (nb==len(mxb)) if not ok : print('Dimensions wrong, n,na,nb='+str([n,na,nb])) wa=len(ma[0])*len(ma[0][0]); va=len(mxa[0]); ok= ok and (va==wa) wb=len(mb[0])*len(mb[0][0]); vb=len(mxb[0]); ok= ok and (vb==wb) if not ok : print('Also wrong, wa,wb,va,vb='+str([wa,wb,va,vb])) bugs= subcheck(mxa,ma)+subcheck(mxb,mb) print('Nb of mismatch: '+str(bugs)) def cgtest(n2=2) : # compare cgm to mxa,mxb, for spin-(1/2) 2nd factors. ok. # n2=2L added to 2S=1, n2=1,2,3... cgm= clebschgordan(n2,1) mxa,s,t =testcleb([n2+1,2],[n2,2]) mxb,s,t =testcleb([n2-1,2],[n2,2]) bugcheck(cgm,mxa,mxb) # special check for spin-(1/2), ok no mismatch. #clebschgordan(1,2); print('') def cg_getrow(cg,i) : # extract full row of CG matrix na=len(cg[0][0]); nb=len(cg[0][0][0]); n=na*nb; z=[[]]*n j=0; h=i; ia=0;ib=0 while h>=len(cg[j]) : h-=len(cg[j]);j+=1 # row is cg[j][h] for k in range(n) : z[k]=cg[j][h][ia][ib]; ib+=1; (ia,ib)= (ia+1,0) if (ib>=nb) else (ia,ib) return z def cg_getcol(cg,i) : na=len(cg[0][0]); nb=len(cg[0][0][0]); n=na*nb; z=[[]]*n ia= div(i,nb); ib=i-ia*nb; j=0;h=0 for k in range(n) : z[k]=cg[j][h][ia][ib]; h+=1; (j,h)= (j+1,0) if (h>=len(cg[j])) else (j,h) return z def orthotest(cg,mode,noisy) : # is cg orthogon.matrix? mode=0 rows =1 columns. na=len(cg[0][0]); nb=len(cg[0][0][0]); m=na*nb; x=[[]]*m hint=['Rows','Cols']; nerr=0;pairs=0 for i in range(m) : x[i]= cg_getrow(cg,i) if(mode==0) else cg_getcol(cg,i) for i in range(m) : for j in range(i,m) : sp=scalprod(x[i],x[j]); pairs+=1 if i==j : sp=lincomb(sp,reformat(1),1,-1) sok= iszero(sp); nerr= (nerr+1) if (not sok) else nerr if noisy : print('Orthonorm '+hint[mode]+' errs: '+str(nerr)+'/'+str(pairs)) return nerr def mirrortest(cg,noisy) : nonz=0 # test symmetry : <l,s;-l3,-s3|j,-j3>=(-1)^{l+s-j}<l,s;l3,s3|j,j3> nj=len(cg); na=len(cg[0][0]); nb=len(cg[0][0][0]); n=na*nb; nerr=0;pairs=0 for ij in range(nj) : mj=len(cg[ij]) for im in range(mj) : for ia in range(na) : for ib in range(nb) : x=cg[ij][im][ia][ib]; y=cg[ij][mj-1-im][na-1-ia][nb-1-ib] j2=mj-1; a2=na-1; b2=nb-1; z=a2+b2-j2; nonz+=(0 if iszero(x) else 1) if (z%2)==1 : print('Bug '+str([j2,a2,b2])); exit() z=div(abs(z),2); sign=1 if ((z%2)==0) else (-1); pairs+=1 w=lincomb(reformat(x),reformat(y),1,-sign); ok= iszero(w) if not ok : nerr +=1; print(' Err'+str([ij,im,ia,ib,x,y])) if noisy : print('Mirror errs/pairs/nonz '+str(nerr)+'/'+str(pairs)+'/'+str(nonz)) return nerr def racah(j2,l2,s2,l32) : # racah formula, alternative j sequence start def facto(n) : # n! y=1; i=n while i>1 : y*=i; i-=1 return y def fac2(n) : # (n/2)! if (n<=2) : return 1 if (n%2)!=0 : print('fac2 wrong argument '+str(n)+'. exit.'); exit() i=div(n,2); return facto(i) s32=j2-l32; pp=facto(j2+1)*fac2(l2+s2-j2) pq= fac2(l2+s2+j2+2)*fac2(j2+l2-s2)*fac2(j2+s2-l2) qp= fac2(l2+l32)*fac2(s2+s32); qq=fac2(l2-l32)*fac2(s2-s32) x=div(l2-l32,2); sign= 1 if (x%2)==0 else -1 z=sqroot([pp*qp,pq*qq]); z[0]=sign*z[0]; return z def racahtest(cg,noisy) : # test all cg[j][0] lines against the formula n=len(cg); j2=len(cg[0])-1; count=0; errs=0 for i in range(n) : row=cg[i][0]; nl=len(row); ns=len(row[0]); l2=nl-1; l32=l2 for k in range(nl) : s2= ns-1; s32=s2 for h in range(ns) : ok= (l32+s32)==j2; z=row[k][h] if ok : y=racah(j2,l2,s2,l32); count+=1; w=lincomb(y,reformat(z),1,-1) if ok : errs = errs if iszero(w) else (errs+1) s32 -=2 l32 -=2 j2 -=2 if noisy : print('Racahtest errors='+str(errs)+'/'+str(count)) return errs def cgexample(l2,s2,noisy) : cgm= clebschgordan(l2,s2,noisy); errs=0; s='' errs += orthotest(cgm,0,noisy) errs += orthotest(cgm,1,noisy) errs += mirrortest(cgm,noisy) errs += racahtest(cgm,noisy) if errs==0 : s=wikidump(cgm); print('<math>'+s+'</math>') return errs,s def fulltest() : tit=['','(1/2)','1','(3/2)','2'] for i in range(1,5) : for k in range(1,5) : title="\n\'\'\'Clebsch-Gordan L="+tit[i]+" S="+tit[k]+"\'\'\'<br>" print(title); cgexample(i,k,False) #errs=cgexample(6,6,True) fulltest()