c=======================================================================
c function  coef9j (it uses the function racah, i.e the 9j symbols
c are computing using Racah coeeficients).
c Numbering convention is used:
c the arguments of coef9j, when numbered sequentially 1 through 9,
c correspond to the array
c                               1  2  3
c                               4  5  6
c                               7  8  9
c
c=======================================================================
      real*8 function coef9j(j1,j2,j3,j4,j5,j6,j7,j8,j9)
      implicit real*8 (a-h,o-z)
      dimension lt(9)
      u9=0.0d0
      lt(1)=j1
      lt(2)=j2
      lt(3)=j4
      lt(4)=j5
      lt(5)=j3
      lt(6)=j6
      lt(7)=j7
      lt(8)=j8
      lt(9)=j9
      lmin=lt(1)
      imin=1
      do 20 i=2,9
      if(lt(i)-lmin) 15,20,20
   15 lmin=lt(i)
      imin=i
   20 continue
      kex=0
      go to (110,110,110,110,150,150,170,170,190),imin
  110 mm=(imin-1)/2+1
      m1=mm+mm-1
      m2=m1+1
      m3=mm+4
      l1=lt(7)
      lt(7)=lt(m1)
      lt(m1)=l1
      l1=lt(8)
      lt(8)=lt(m2)
      lt(m2)=l1
      l1=lt(9)
      lt(9)=lt(m3)
      lt(m3)=l1
      imin=imin+(7-m1)
      go to 175
  150 kex=1
      m1=7
      m2=8
      m3=imin+imin-9
      m4=m3+1
      go to 180
  170 kex=1
  175 m1=5
      m2=6
      m3=imin-6
      m4=m3+2
  180 l1=lt(m1)
      l1=lt(m1)
      lt(m1)=lt(m3)
      lt(m3)=l1
      l1=lt(m2)
      lt(m2)=lt(m4)
      lt(m4)=l1
      l1=lt(9)
      lt(9)=lt(imin)
      lt(imin)=l1
  190 if(lt(9)) 200,200,300
  200 if(lt(5)-lt(6)) 1000,210,1000
  210 if(lt(7)-lt(8)) 1000,220,1000
  220 rt=(lt(5)+1)*(lt(7)+1)
      k=(lt(5)+lt(7)-lt(1)-lt(4))/2
      rac= racah(lt(1),lt(2),lt(3),lt(4),lt(5),lt(7))
      ph=1.0d0
      if (2*(k/2) .ne. k) ph=-1.0d0
      u9=(rac/dsqrt(rt))*ph
      go to 370
  300 k1=iabs(lt(2)-lt(7))
      k2=iabs(lt(3)-lt(5))
      k3=iabs(lt(4)-lt(9))
      nmin=max0(k1,k2,k3)
      k1=lt(2)+lt(7)
      k2=lt(3)+lt(5)
      k3=lt(4)+lt(9)
      nmax=min0(k1,k2,k3)
      if (nmin-nmax) 320, 320, 1000
 320  do 350 n=nmin,nmax,2
        w1=n+1
        rac= racah(lt(2),lt(5),lt(7),lt(3),lt(1),n)
        if (rac) 321, 350, 321
 321    w1=w1*rac
        rac= racah(lt(2),lt(4),lt(7),lt(9),lt(8),n)
        if (rac) 322, 350, 322
 322    w1=w1*rac
        rac= racah(lt(3),lt(4),lt(5),lt(9),lt(6),n)
        if (rac) 323, 350, 323
 323    u9=u9+w1*rac
 350  continue
 370  if(kex) 400,1000,400
 400  kp=0
      do 410 i=1,9
 410  kp=kp+lt(i)
      k=kp/2
      ph=1.0d0
      if (2*(k/2) .ne. k) ph=-1.0d0
      u9=u9*ph
 1000 coef9j=u9
      return
      end


C     *********************************************************
c     calculates racah coefficients
c	Rayner and Kazuhito Dec 2008
c     *********************************************************
c     calculates racah coefficients
c     entries : racah , racahi , racahr
c     racah  : integer arguments = 2*j
c     racahi : arguments = true integer value
c     racahr : arguments = true real value
c     external : th_facinit , generates factorial table (ok)
C     *********************************************************             
      real*8 function racah(jad,jbd,jcd,jdd,jed,jfd)
      implicit real*8 (a-h,o-z)
      dimension i(16)
      real*8 g,s
      parameter (lfactc=2000)
      logical first
      common / logfac / first,lfact,g(lfactc)
      equivalence(i(1),i1),(i(2),i2),(i(3),i3),(i(4),i4),(i(5),i5),
     1     (i(6),i6),(i(7),i7),(i(8),i8),(i(9),i9),(i(10),i10),(i(11)
     $     ,i11),(i(12),i12),(i(13),i13),(i(14),i14),(i(15),i15),(i(16)
     $     ,i16)

c     make usefull combinations
      k=jad+jbd-jed+2
      i1=k/2
      if((2*i1).ne.k) goto 300
      k=jcd+jdd-jed+2
      i4=k/2
      if((2*i4).ne.k) goto 300
      k=jad+jcd-jfd+2
      i7=k/2
      if((2*i7).ne.k) goto 300
      k=jbd+jdd-jfd+2
      i10=k/2
      if((2*i10).ne.k) goto 300
      i13=i1+jed
      i14=i4+jed
      i15=i7+jfd
      i16=i10+jfd
      i2=i13-jad
      i3=i13-jbd
      i5=i14-jcd
      i6=i14-jdd
      i8=i15-jad
      i9=i15-jcd
      i11=i16-jbd
      i12=i16-jdd
c     check triangular inequalities,find no. of terms in sum
      n=min0(i1,i2,i3,i4,i5,i6,i7,i8,i9,i10,i11,i12)-1
      if(n) 300,2,2
c     find minimum value of summation index
    2 il=max(i13,i14,i15,i16)
      if(min0(jad,jbd,jcd,jdd,jed,jfd)) 300,20,1
c     ..............
      entry racahi(ja1,jb1,jc1,jd1,je1,jf1)
c     make usefull combinations
      i13=ja1+jb1+je1+1
      i14=jc1+jd1+je1+1
      i15=ja1+jc1+jf1+1
      i16=jb1+jd1+jf1+1
      i1=i13-je1*2
      i2=i13-ja1*2
      i3=i13-jb1*2
      i4=i14-je1*2
      i5=i14-jc1*2
      i6=i14-jd1*2
      i7=i15-jf1*2
      i8=i15-ja1*2
      i9=i15-jc1*2
      i10=i16-jf1*2
      i11=i16-jb1*2
      i12=i16-jd1*2
c     check triangular inequalities,find no. of terms in sum
      n=min0(i1,i2,i3,i4,i5,i6,i7,i8,i9,i10,i11,i12)-1
      if(n) 300,4,4
c     find minimum value of summation index
    4 il=max(i13,i14,i15,i16)
      lmin=min0(ja1,jb1,jc1,jd1,je1,jf1)
      if(lmin)300,20,1
c     ............
      entry racahr(a,b,c,d,e,f)
c     convert arguments to integer 
      ja=dnint(2.0d0*a)
      jb=dnint(2.0d0*b)
      jc=dnint(2.0d0*c)
      jd=dnint(2.0d0*d)
      je=dnint(2.0d0*e)
      jf=dnint(2.0d0*f)
c     make usefull combinations
      k=ja+jb-je+2
      i1=k/2
      if((2*i1-k).ne.0) goto 300
      k=jc+jd-je+2
      i4=k/2
      if((2*i4-k).ne.0) goto 300
      k=ja+jc-jf+2
      i7=k/2
      if((2*i7-k).ne.0) goto 300
      k=jb+jd-jf+2
      i10=k/2
      if((2*i10-k).ne.0) goto 300
      i13=i1+je
      i14=i4+je
      i15=i7+jf
      i16=i10+jf
      i2=i13-ja
      i3=i13-jb
      i5=i14-jc
      i6=i14-jd
      i8=i15-ja
      i9=i15-jc
      i11=i16-jb
      i12=i16-jd
c     check triangular inequalities,find no. of terms in sum
      n=min0(i1,i2,i3,i4,i5,i6,i7,i8,i9,i10,i11,i12)-1
      if(n) 300,3,3
c     find minimum value of summation index
    3 il=max(i13,i14,i15,i16)
      lmin=min0(ja,jb,jc,jd,je,jf)
      if(lmin)300,20,1

    1 if(first) call th_facinit
      if(il.ge.lfact) stop 'racah: length factorial table insufficient'
      j1=il-i13+1 
      j2=il-i14+1 
      j3=il-i15+1 
      j4=il-i16+1
      j5=i13+i4-il 
      j6=i15+i5-il 
      j7=i16+i6-il
      ph=1.0d0
      if(2*(j5/2).eq.j5) ph=-1.0d0
      h=ph*dexp((g(i1)+g(i2)+g(i3)-g(i13+1)+g(i4)+g(i5)+g(i6)-
     1     g(i14+1)+g(i7)+g(i8)+g(i9)-g(i15+1)+g(i10)+g(i11)+g(i12)
     $     -g(i16+1))*.50d0+g(il+1)-g(j1)-g(j2)-g(j3)-g(j4)-g(j5)-g(j6)
     $     -g(j7))
      if(n)300,110,120
c
 110  racah=h 
      return
c
 120  s=1.0d0
      k=n-1
      kl=il+1
      j5=j5-1
      j6=j6-1
      j7=j7-1
      do 130 j=1,n              ! k=n-j
        s=1.0d0
     &       -((kl+k)*(j5-k)*(j6-k)*(j7-k))*s/((j1+k)*(j2+k)*(j3+k)*(j4
     $       +k))
        k=k-1
 130  continue  
      racah=h*s
      return
c     one of the arguments =0
 20   iad=il
      ibd=il
      do 21 j=13,16
        if(iad.lt.i(j)) goto 22
        if(iad.lt.ibd) ibd=iad
        iad=i(j)
        goto 21
 22     if(ibd.gt.i(j)) ibd=i(j)
 21   continue
      j5=i13+i4-il 
      ph=1.0d0
      if(2*(j5/2).eq.j5) ph=-1.0d0
      racah=ph/dsqrt(dfloat(iad*ibd))
      return
c     impossible combination of arguments
 300  racah=0.0d0
      return
      end      


c      *********************************************************
c     th_facinit , generates factorial table
c	Rayner and Kazuhito Dec 2008
C     ********************************************************* 
      subroutine th_facinit
      implicit real*8 (a-h,o-z)
      parameter (lfactc=2000)
      logical first
      real*8 faclog,fn
      common / logfac / first,lfact,faclog(lfactc)
      data first/.true./,lfact/lfactc/

      first=.false.
      faclog(1)=0.0d0
      faclog(2)=0.0d0
      fn=1.0d0
      do 10 i=3,lfactc
        fn=fn+1.0d0
        faclog(i)=faclog(i-1)+dlog(fn)
 10   continue
      return
      end      


c ********************************
c sixj with integer input (2*i)
c given in terms of a racah coeff
c	Rayner and Kazuhito Dec 2008
c ********************************
      real*8 function sj2i(j1,j2,k2,k1,j3,k3)
      implicit real*8 (a-h,o-z)
      sj2i = par(dfloat(j1+j2+j3+k1)/2.0d0)*
     >        racah(j1,j2,j3,k1,k2,k3)
      return
      end
cc ********************************
c parity, i.e. (-1.)**a
c	Rayner and Kazuhito Dec 2008
cc ********************************
      real*8 function par(a)
      implicit real*8 (a-h,o-z)
c***********************************
      i = a
      i = i/2
      p = a/2.0d0
      q = i
      if(p-q) 1,2,1
 1    par = -1.0d0
      return
 2    par = 1.0d0
      return
      end            
cc***********************************
c three-j symbol
c given in terms of a clebsh
c	Rayner and Kazuhito Dec 2008
cc***********************************
      real*8 function tj2i(j21,j22,j23,m21,m22,m23)
      implicit real*8 (a-h,o-z)
      a1=dfloat(j21)/2.0d0
      b1=dfloat(m21)/2.0d0
      a2=dfloat(j22)/2.0d0
      b2=dfloat(m22)/2.0d0
      a3=dfloat(j23)/2.0d0
      b3=dfloat(m23)/2.0d0
      goto 9
      entry tji(j1,j2,j3,m1,m2,m3)
      a1=dfloat(j1)
      b1=dfloat(m1)
      a2=dfloat(j2)
      b2=dfloat(m2)
      a3=dfloat(j3)
      b3=dfloat(m3)
      entry tj(aa1,aa2,aa3,bb1,bb2,bb3)
      a1=aa1
      b1=bb1
      a2=aa2
      b2=bb2
      a3=aa3
      b3=bb3
 9    tj = par(a1-a2-b3)*clebr(a1,b1,a2,b2,a3,-b3)/
     &     dsqrt(2.0d0*a3+1.0d0)
      return
      end
      
c*************************************************************
c  cleb, clebi, cleb2i
c	Rayner and Kazuhito Dec 2008
c*************************************************************
      real*8 function clebr(a,b,c,d,e,f)
c     arguments are real and of true value; j1,m1,j2,m2,j3,m3
      implicit real*8 (a-h,o-z)
      real*8 faclog
      parameter (lfactc=2000)
      logical first
      common / logfac / first,lfact,faclog(lfactc)
      ia=dnint(2.0d0*a)
      ib=dnint(2.0d0*c)
      ic=dnint(2.0d0*e)
      id=dnint(2.0d0*b)
      ie=dnint(2.0d0*d)
      if=dnint(2.0d0*f)
      goto 7000
c ...............  clebi  ...........................
c      arguments are integer and of true value
      entry clebi2(ll1,lm1,ll2,lm2,ll3,lm3)
      ia=2*ll1
      ib=2*ll2
      ic=2*ll3
      id=2*lm1
      ie=2*lm2
      if=2*lm3
      goto 7000
c ..............  cleb  ..............................
c      arguments are integer and represent twice the real value
      entry cleb(i2j1,i2m1,i2j2,i2m2,i2j3,i2m3)
      ia=i2j1
      ib=i2j2
      ic=i2j3
      id=i2m1
      ie=i2m2
      if=i2m3
 7000 if (first) call th_facinit
      rac=0.0d0
      if(id+ie-if) 1000,105,1000
 105  k1=ia+ib+ic
      if((-1)**k1) 1000,110,110
 110  k1=ia+ib-ic
      k2=ic+ia-ib
      k3=ib+ic-ia
      k4=ia-iabs (ib-ic)
      k5=ib-iabs (ic-ia)
      k6=ic-iabs (ia-ib)
      k7= min0 (k1,k2,k3,k4,k5,k6)
      if(k7) 1000,120,120
 120  if((-1)**(ia+id)) 1000,1000,130
 130  if((-1)**(ib+ie)) 1000,1000,140
 140  if((-1)**(ic+if)) 1000,1000,150
 150  if(ia-iabs (id)) 1000,152,152
 152  if(ib-iabs (ie)) 1000,154,154
 154  if(ic-iabs (if)) 1000,160,160
 160  signfc=1.0d0
      iam=ia
      ibm=ib
      icm=ic
      idm=id
      iem=ie
      ifm=if
      if(ia-ib) 210,220,220
 210  if(ia-ic) 215,225,225
 215  it=ia
      ia=ib
      ib=it
      it=id
      id=ie
      ie=it
      signfc=(-1.0d0)**((ia+ib-ic)/2)
      go to 235
 220  if(ic-ib) 225,235,235
 225  it=ic
      ic=ib
      ib=it
      it=if
      if=-ie
      ie=-it
      fibm=ibm+1
      ficm=icm+1
      signfc=(-1.d0)**((iam-idm)/2)*dsqrt (ficm/fibm)
 235  if(ib) 237,236,237
 236  rac=signfc
      go to 900
 237  if(ie) 250,250,240
 240  signfc=signfc*((-1.0d0)**((ia+ib-ic)/2))
      id=-id
      ie=-ie
      if=-if
 250  fc2=ic+1
      iabcp=(ia+ib+ic)/2+1
      iabc=iabcp-ic
      icab=iabcp-ib
      ibca=iabcp-ia
      iapd=(ia+id)/2+1
      iamd=iapd-id
      ibpe=(ib+ie)/2+1
      ibme=ibpe-ie
      icpf=(ic+if)/2+1
      icmf=icpf-if
      sqfclg=0.5d0*(dlog(fc2)-faclog(iabcp+1)
     1     +faclog(iabc)+faclog(icab)+faclog(ibca)
     2     +faclog(iapd)+faclog(iamd)+faclog(ibpe)
     3     +faclog(ibme)+faclog(icpf)+faclog(icmf))
      nzmic2=(ib-ic-id)/2
      nzmic3=(ia-ic+ie)/2
      nzmi= max0 (0,nzmic2,nzmic3)+1
      nzmx= min0 (iabc,iamd,ibpe)
      if(nzmi-nzmx) 310,310,900
 310  ss=0.0d0
      s1=(-1.0d0)**(nzmi-1)
      do 400 nz=nzmi,nzmx
        nzm1=(nz-1)
        nzt1=iabc-nzm1
        nzt2=iamd-nzm1
        nzt3=ibpe-nzm1
        nzt4=nz-nzmic2
        nzt5=nz-nzmic3
        termlg=sqfclg-faclog(nz)-faclog(nzt1)-faclog(nzt2)
     1       -faclog(nzt3)-faclog(nzt4)-faclog(nzt5)
        ssterm=s1*dexp(termlg)
        ss=ss+ssterm
 400  s1=-s1
      rac=signfc*ss
 900  ia=iam
      ib=ibm
      ic=icm
      id=idm
      ie=iem
      if=ifm
 1000 cleb=rac
      return
      end