c--------------------------> by Lorenzo Fortunato Apr. 2002 function RACAH3J(aj1,aj2,aj3,bm1,bm2,bm3,err3) c---------------------------------------------------------| c This function calculates the Racah V coefficient (also| c called Wigner 3J) as stated in closed formula 15.36 | c in the book "Nuclear Shell Theory" by deShalit-Talmi | c (Academic press 1963). Uses function fuckt. | c-> The first 6 entries are real number 1st row and 2nd | c consecutively. | c err9=FALSE turns off the triangle rule check | c Written by Lorenzo Fortunato in April 2002. | c Improved by Antonio Mason in April 2006. | c Updated in Jan. 2007 | c---------------------------------------------------------| implicit real(a-h,o-z) logical err3 racah3j=0. som=bm1+bm2+bm3 if (som.ne.0.) then if(err3.eqv..TRUE.)write(*,*)" m1+m2+m3 is not zero!" goto 120 else if (aj3.gt.(aj1+aj2).or.aj3.lt.abs(aj1-aj2)) then if(err3.eqv..TRUE.)write(*,*)"Triangle rule not satisfied!" goto 120 endif sq1=sqrt(fuckt(aj1+aj2-aj3)*fuckt(aj2+aj3-aj1)* . fuckt(aj3+aj1-aj2)/fuckt(aj1+aj2+aj3+1.)) sq2=sqrt(fuckt(aj1+bm1)*fuckt(aj1-bm1)* . fuckt(aj2+bm2)*fuckt(aj2-bm2)* . fuckt(aj3+bm3)*fuckt(aj3-bm3)) smt=0. do 100 t=max((aj2 -aj3-bm1),(-aj3+aj1+bm2)), & min((aj1+aj2-aj3),(aj1-bm1),(aj2+bm2)),1. term=1./(fuckt(t)*fuckt(aj1+aj2-aj3-t)* . fuckt(aj3-aj2+bm1+t)*fuckt(aj3-aj1-bm2+t)* . fuckt(aj1-bm1-t)*fuckt(aj2+bm2-t)) 50 chksgn=aj1-aj2-bm3+t segno=(-1)**(mod(abs(chksgn),2.)) smt=smt+segno*term oldsmt=smt 100 enddo 110 racah3j=sq1*sq2*smt 120 return end function RACAH6J(aj1,aj2,aj3,bl1,bl2,bl3,err6) c---------------------------------------------------------| c This routine calculates the Racah W coefficient (also | c called Wigner 6J) as stated in closed formula 15.38 | c in the book "Nuclear Shell Theory" by deShalit-Talmi | c (Academic press 1963). Uses function fuckt and Delta. | c-> Entries are real number 1st row and 2nd consecutively.| c err9=FALSE turns off the triangle rule check | c Written by Lorenzo Fortunato in April 2002. | c Improved by Antonio Mason in April 2006. | c Error option included by L.F. in January 2006 | c---------------------------------------------------------| implicit real(a-h,o-z) logical err6 racah6j=0. if (aj3.lt.abs(aj1-aj2).or.aj3.gt.(aj1+aj2)) then if (err6.eqv..TRUE.) write(*,*)'1,2,3 don`t obey trgl.rule' goto 110 elseif (bl3.lt.abs(aj1-bl2).or.bl3.gt.(aj1+bl2)) then if (err6.eqv..TRUE.) write(*,*)'1,5,6 don`t obey trgl. rule' goto 110 elseif (bl3.lt.abs(bl1-aj2).or.bl3.gt.(bl1+aj2)) then if (err6.eqv..TRUE.) write(*,*)'2,4,6 don`t obey trgl. rule' goto 110 elseif (aj3.lt.abs(bl1-bl2).or.aj3.gt.(bl1+bl2)) then if (err6.eqv..TRUE.) write(*,*)'3,4,5 don`t obey trgl. rule' goto 110 endif coef=Delta(aj1,aj2,aj3)*Delta(aj1,bl2,bl3)* . Delta(bl1,aj2,bl3)*Delta(bl1,bl2,aj3) somt=0. bo1=max((aj1+aj2+aj3),(aj1+bl2+bl3),(bl1+aj2+bl3),(bl1+bl2+aj3)) bo2=min((aj1+aj2+bl1+bl2),(aj2+aj3+bl2+bl3),(aj3+aj1+bl3+bl1)) if (bo1.ne.aint(bo1)) bo1=aint(bo1)+1. if (bo2.ne.aint(bo2)) bo2=aint(bo2) do t=bo1,bo2,1. segno=(-1)**(mod(t,2)) term2=fuckt(1+t)/(fuckt(t-aj1-aj2-aj3)* . fuckt(t-aj1-bl2-bl3)*fuckt(t-bl1-aj2-bl3)* . fuckt(t-bl1-bl2-aj3)*fuckt(aj1+aj2+bl1+bl2-t)* . fuckt(aj2+aj3+bl2+bl3-t)*fuckt(aj3+aj1+bl3+bl1-t)) if (abs(fuckt(t+1)/term2).le.1.e-30) goto 100 somt=somt+segno*term2 c print *,'t,somt,segno*term2=',t,somt,segno*term2 90 enddo 100 racah6j=coef*somt 110 return end function fuckt(a) c-------------------------------------------------------------| c Very simple-minded factorial (OK for not too big entries) | c It can handle half integer (to a limited precision!). | c-------------------------------------------------------------| implicit real(a-h,o-z) pi=3.14159265359 chk=a-aint(a) rst=1. if (chk.eq.0.5) then do cc=a,0.5,-1. rst=rst*cc enddo fuckt=rst*pi**.5 goto 500 else do cc=a,1.,-1. rst=rst*cc enddo fuckt=rst goto 500 endif 500 return end function Delta(a,b,c) c-------------------------------------------------------------| c Delta function as in "Nuclear Shell theory" deShalit-Talmi | c 15.39 (Academic Press 1963). Essential for 6J oefficients. | c-------------------------------------------------------------| implicit real(a-h,o-z) Delta=sqrt(fuckt(a+b-c)*fuckt(b+c-a)*fuckt(c+a-b)/ . fuckt(a+b+c+1.)) return end function cleb(sj1,sj2,sm1,sm2,rJ,rM) c-------------------------------------------------------------| c This routine calculates the Clebsch-Gordan vector coupling | c coefficients by calling the racah3J routine written above. | c It is an inversion of formula 13.52 in the cited book. | c-> It uses also fuckt function. The s-entries are the left | c ones in the CG while the r-entries are the right ones. | c Here m1+m2=M and in racah3j we have m1+m2+m3=0 !! | c Be carefull to the minus!!If one exchanges j1 and j2 -> | c (-1)**(J-j1-j2)Written by Lorenzo Fortunato in April 2002. | c-------------------------------------------------------------| implicit real(a-h,o-z) segno=-1. c print *,'segno CG ',sj1-sj2+rM if (mod(abs(sj1-sj2+rM),2.).eq.0.) then segno=1. endif rrr=racah3J(sj1,sj2,rJ,sm1,sm2,-rM,.FALSE.) cleb=segno*sqrt(2.*rJ+1.)*rrr return end function RACAH9J(aj1,aj2,aj3, . bj1,bj2,bj3, . cj1,cj2,cj3,err9) c-------------------------------------------------------------| c This function calculates the Racah X coefficient (also | c called Wigner 9J)as it appears in Jahn,Hope, PRC93(1954) | c The triangle rule is checked for every row and column. | c err9=FALSE turns off the triangle rule check | c Written by Lorenzo Fortunato in Apr. 2002.Update Jan.2007 | c-------------------------------------------------------------| implicit real(a-h,o-z) logical err9 external racah6j w9j=0. if (aj3.lt.abs(aj1-aj2).or.aj3.gt.(aj1+aj2)) then if (err9.eqv..TRUE.)write(*,*)'1st row don`t obey trgle rule' goto 100 elseif (bj3.lt.abs(bj1-bj2).or.bj3.gt.(bj1+bj2)) then if (err9.eqv..TRUE.)write(*,*)'2nd row don`t obey trgle rule' goto 100 elseif (cj3.lt.abs(cj1-cj2).or.cj3.gt.(cj1+cj2)) then if (err9.eqv..TRUE.)write(*,*)'3rd row don`t obey trgle rule' goto 100 elseif (cj1.lt.abs(aj1-bj1).or.cj1.gt.(aj1+bj1)) then if (err9.eqv..TRUE.)write(*,*)'1st col don`t obey trgle rule' goto 100 elseif (cj2.lt.abs(aj2-bj2).or.cj2.gt.(aj2+bj2)) then if (err9.eqv..TRUE.)write(*,*)'2nd col don`t obey trgle rule' goto 100 elseif (cj3.lt.abs(aj3-bj3).or.cj3.gt.(aj3+bj3)) then if (err9.eqv..TRUE.)write(*,*)'3rd col don`t obey trgle rule' goto 100 endif xmi=min(abs(cj2-cj3),abs(bj1-bj3),abs(aj1-aj2)) xma=max((cj2+cj3),(bj1+bj3),(aj1+aj2)) do x=xmi,xma,.5 rac1=racah6j(aj1,bj1,cj1,cj2,cj3,x,.FALSE.) rac2=racah6j(aj2,bj2,cj2,bj1,x,bj3,.FALSE.) rac3=racah6j(aj3,bj3,cj3,x,aj1,aj2,.FALSE.) term=rac1*rac2*rac3*(2*x+1)*(-1.)**(2*x) w9j=w9j+term enddo c 10 continue 100 racah9j=w9j return end function elsphar(aIf,aMf,alam,amu,aIi,aMi) c-------------------------------------------------------------| c This function calculates the matrix element of lambda-mu | c spherical harmonic between an initial (right) state and | c a final (left) state. The entries are real numbers. | c Needs the functions racah3J and fuckt. | c Written by Lorenzo Fortunato in May 2002. | c-------------------------------------------------------------| implicit real(a-h,o-z) pi=3.14159265359 co1= RACAH3J(aIf,alam,aIi,-aMf,amu,aMi,.FALSE.) co2= RACAH3J(aIf,alam,aIi,0.,0.,0.,.FALSE.) if (mod((-aMf),2.).eq.0.) then segn=1. else segn=-1. endif elsphar=segn*co1*co2* . sqrt((2.*aIi+1.)*(2.*alam+1.)*(2.*aIf+1.)/(4.*pi)) return end function relsph(aIf,alam,aIi) c-------------------------------------------------------------| c This function calculates the reduced matrix element of | c lambda spherical harmonic between an initial (right) state | c and a final (left) state. The entries are real numbers. | c Taken from Rotemberg's tables.Care to the sign convention. | c Needs the functions racah3J and fuckt. | c Written by Lorenzo Fortunato in June 2002. | c-------------------------------------------------------------| implicit real(a-h,o-z) pi=3.14159265359 co2= RACAH3J(aIf,alam,aIi,0.,0.,0.,.FALSE.) if (mod(aIf,2.).eq.0.) then segn=1. else segn=-1. endif relsph=segn*co2* . sqrt((2.*aIi+1.)*(2.*alam+1.)*(2.*aIf+1.)/(4.*pi)) return end