program cfp
c     coefficients of fractional parentage
c      II: Racah, G. Phys. Rev. 62, 1942, 438-462. 
c     III: Racah, G. Phys. Rev. 63, 1943, 367-382.

      implicit none
      real*8 a,b,c,d,e,f,sw,bw,S,L,Sp,Lp,ll,f6,Spp,Lpp,lmax,k
      integer*4 ic,n,iwr,dstates
      character*5 oc(2)
      data oc/'large','small'/

      call InitFac
      iwr=1
      
1000  write(6,6001)oc(iwr)
6001  format(/,
     1'  1    w(abcd;ef)',/,
     1'  2    W(abcd;ef)',/,
     1'  3    (6):',/,
     1'       (s1l1s2l2(SpLp)s3l3SL|s1l1,s2l2s3l3(SppLpp),SL)',/,
     1'  4    (9) equation set for CFP:',/,
     1'       sum(SpLp)(l,ll(SppLpp),SL|l2(SpLp)lSL)',/,
     1'           x CFP(l2(SpLp)lSL||l3aSL)=0 Spp+Lpp odd',/,
     1'  5    (11) equation set for CFP, calculate',/,
     1'           CFP(l^(n-1)(SpLp) l SL]|l^n aSL)',/,
     1'  6    maximum L for n electrons in l-shell',/,
     1'  7    (SL||U(k)||SpLp) for n-electrons in l-shell',/,
     1'  8    output control: ',a5/,
     1'  9    (SL||V(k)||SpLp) for n-electrons in l-shell',/,
     1' 10    investigate states in a l^n configuration',/,
     1' 11    degeneracy of |S L> state of l^n configuration',/,
     1'  0    quit',/,
     1'  ---------',/)
      read(5,*)ic
      if(ic.eq.0)stop
      if(ic.eq.1)then
       write(6,501)
501    format(' a b c d e f: ',$)
       read(5,*)a,b,c,d,e,f,f
       write(6,*)sw(a,b,c,d,e,f)
      endif
      if(ic.eq.2)then
       write(6,501)
       read(5,*)a,b,c,d,e,f,f
       write(6,*)bw(a,b,c,d,e,f)
      endif
      if(ic.eq.3)then
       write(6,503)
503    format(' l Sp Lp S L Spp Lpp: ',$)
       read(5,*)ll,Sp,Lp,S,L,Spp,Lpp
       write(6,*)f6(ll,Sp,Lp,S,L,Spp,Lpp)
      endif
      if(ic.eq.4)then
       write(6,504)
504    format(' l S L: ',$)
       read(5,*)ll,S,L
       call e9(ll,S,L)
      endif
      if(ic.eq.5)then
       write(6,505)
505    format(' n l S L: ',$)
       read(5,*)n,ll,S,L
       call e11(n,ll,S,L)
      endif
      if(ic.eq.6)then
       write(6,506)
506    format(' n l: ',$)
       read(5,*)n,ll
       write(6,*)' Max L = ',NINT(lmax(ll,n))
      endif
      if(ic.eq.7)then
       write(6,507)
507    format(' n l S L Sp Lp k: ',$)
       read(5,*)n,ll,S,L,Sp,Lp,k
       call U44(n,ll,S,L,Sp,Lp,k,iwr)
      endif
      if(ic.eq.8)then
       iwr=iwr+1
       if(iwr.gt.2)iwr=1
      endif
      if(ic.eq.9)then
       write(6,507)
       read(5,*)n,ll,S,L,Sp,Lp,k
       call V44(n,ll,S,L,Sp,Lp,k,iwr)
      endif
      if(ic.eq.10)then
       write(6,510)
510    format(' n l: ',$)
       read(5,*)n,ll
       call states(n,ll,iwr)
      endif
      if(ic.eq.11)then
       write(6,511)
511    format(' n l S L: ',$)
       read(5,*)n,ll,S,L
       write(6,5111)dstates(n,ll,0,S,L)
5111   format(' Ndeg =',i2)
      endif
      goto 1000
      end

      function dstates(nn,ll,iwr,Slook,Llook)
c     all states in l-shell for n-electrons
      implicit none
      integer*4 nn,il,imax,i,j,n,m,NHBIG,iex,ic,orbit,spin,ML,MS,
     1smax,mm,smin,S,L,ns,is,ist,iwr,dstates,it,islook,illook
      real*8 ll,lmax,Slook,Llook
      character*1 ssy(17)
      data ssy/'s','p','d','f','g','h','i','k','l','m','n','o','q','r',
     1         't','u','v'/
      character*1 sy(17)
      data  sy/'S','P','D','F','G','H','I','K','L','M','N','O','Q','R',
     1         'T','U','V'/
      integer*4,allocatable::tab(:,:),STATE(:),si(:),lstab(:,:),
     1mtab(:,:)
      it=0
      islook=NINT(2.0d0*Slook)
      illook=NINT(Llook)
      smin=mod(nn,2)
      il=NINT(ll)
      imax=NINT(lmax(ll,nn))
      smax=min(nn,2*il+1)
      if(nn.gt.2*il+1)smax=4*il+2-nn
      if(iwr.eq.1)write(6,6001)ssy(il+1),nn
6001  format(/,a1,'^',i1,/)
c     -Lmax ...Lmax (ML) values:
      m=2*imax+1
c     -Smax ...Smax (MS) values (multiple of 1/2):
      n=smax+1
c     2*m ... number of spinorbitals
      if(iwr.eq.1)write(6,6009)nn,il,imax,smax
6009  format(i2,' electrons,',i2,' = l,',' Lmax =',i2,
     1' Smax =',i2,'/2')
      mm=2*(2*il+1)
      allocate (tab(m,n),STATE(mm),si(nn),lstab(m,n),mtab(m,n))
      do 1 i=1,m
      do 1 j=1,n
1     tab(i,j)=0
      NHBIG=0
c
c     states Nexc excited, Nexc=nn
c     distribute Nexc excitations upon centers 1...mm (N7..N)
c     initial indices - all to the first center:
      do 41 iex=1,nn
41    si(iex)=1
5     DO  1966 I=1,mm
1966  STATE(I)=0
      do 4 iex=1,nn
4     state(si(iex))=state(si(iex))+1
c     check that the state is allowed:
      do 42 I=1,mm
42    if(STATE(I).gt.1)goto 12

      MS=0
      ML=0
      do 11 I=1,mm
      spin=1-2*mod(I,2)
      orbit=-il+(I-1)/2
      MS=MS+spin *STATE(I)
11    ML=ML+orbit*STATE(I)
      if(iwr.eq.2)then
       if(mod(MS,2).eq.0)then
        write(6,6002)MS/2,ML,(STATE(I),I=1,mm)
6002    format(' MS =',i5,' ML =',i4,2x,30i1)
       else
        write(6,6012)MS,ML,(STATE(I),I=1,mm)
6012    format(' MS =',i3,'/2 ML =',i4,2x,30i1)
       endif
      endif
      i=ML+imax+1
      j=(MS+smax+2)/2
      tab(i,j)=tab(i,j)+1
      NHBIG=NHBIG+1
c
c     find index to be changed
12    do 8 ic=nn,1,-1
8     if(si(ic).lt.mm)goto 9
      goto 3
9     do 10 i=ic+1,nn
10    si(i)=si(ic)+1
      si(ic)=si(ic)+1
      goto 5
c
3     if(iwr.eq.2)then
       write(6,6003)NHBIG
6003   format(i6,' states',/,' Table:')
       do 13 i=1,m
13     write(6,6004)-imax+i-1,(tab(i,j),j=1,n)
6004   format(i3,1x,50i2)
      endif

c     extract S^L states:
      is=0
      ist=0
      do 14 S=smax,smin,-2
      do 14 L=imax,0,-1
c     fill SL table:
      do 15 i=1,m
      do 15 j=1,n
15    lstab(i,j)=0
      do 16 ML=-L,L
      do 16 MS=-S,S,2
      i=ML+imax+1
      j=(MS+smax+2)/2
16    lstab(i,j)=lstab(i,j)+1
      if(iwr.eq.2)then
       write(6,*)'LS Table:'
       do 20 i=1,m
20     write(6,6004)-imax+i-1,(lstab(i,j),j=1,n)
      endif

      ns=0
      ic=0
c     write total to a temporary one: 
172   do 18 i=1,m
      do 18 j=1,n
18    mtab(i,j)=tab(i,j)
c     try to subtract lstab:      
      do 17 i=1,m
      do 17 j=1,n
      if(lstab(i,j).gt.0)then
       if(mtab(i,j).gt.0)then
        mtab(i,j)=mtab(i,j)-1
       else
        ic=1
        goto 171
       endif
      endif
17    continue
171   if(ic.eq.0)then
       if(iwr.eq.2)write(6,*)'can'
c      state could be subtracted, update tab and try again:
       do 19 i=1,m
       do 19 j=1,n
19     tab(i,j)=mtab(i,j)
       if(iwr.eq.2)then
        write(6,*)'Updated Table:'
        do 21 i=1,m
21      write(6,6004)-imax+i-1,(tab(i,j),j=1,n)
       endif
       ns=ns+1
       ist=ist+1
       goto 172
      else
       if(iwr.eq.2)write(6,*)'cannot'
      endif
      if(ns.gt.0)then
       is=is+1
       if(illook.eq.L.and.islook.eq.S)it=ns
       if(iwr.eq.1)write(6,6007)is,S+1,sy(L+1),ns
6007   format(i5,i3,'^',a1,i4,' x')
      endif
14    continue
      if(iwr.eq.1)write(6,6008)is,ist,NHBIG
6008  format(/,i5,' SL labels',/,i5,' SL states',/,
     1i5,' configurations')
      dstates=it

      return
      end

      subroutine states(nn,ll,iwr)
c     all states in l-shell for n-electrons
      implicit none
      integer*4 nn,il,imax,i,j,n,m,NHBIG,iex,ic,orbit,spin,ML,MS,
     1smax,mm,smin,S,L,ns,is,ist,iwr
      real*8 ll,lmax
      character*1 ssy(17)
      data ssy/'s','p','d','f','g','h','i','k','l','m','n','o','q','r',
     1         't','u','v'/
      character*1 sy(17)
      data  sy/'S','P','D','F','G','H','I','K','L','M','N','O','Q','R',
     1         'T','U','V'/
      integer*4,allocatable::tab(:,:),STATE(:),si(:),lstab(:,:),
     1mtab(:,:)
      smin=mod(nn,2)
      il=NINT(ll)
      imax=NINT(lmax(ll,nn))
      smax=min(nn,2*il+1)
      if(nn.gt.2*il+1)smax=4*il+2-nn
      write(6,6001)ssy(il+1),nn
6001  format(/,a1,'^',i1,/)
c     -Lmax ...Lmax (ML) values:
      m=2*imax+1
c     -Smax ...Smax (MS) values (multiple of 1/2):
      n=smax+1
c     2*m ... number of spinorbitals
      if(iwr.eq.1)write(6,6009)nn,il,imax,smax
6009  format(i2,' electrons,',i2,' = l,',' Lmax =',i2,
     1' Smax =',i2,'/2')
      mm=2*(2*il+1)
      allocate (tab(m,n),STATE(mm),si(nn),lstab(m,n),mtab(m,n))
      do 1 i=1,m
      do 1 j=1,n
1     tab(i,j)=0
      NHBIG=0
c
c     states Nexc excited, Nexc=nn
c     distribute Nexc excitations upon centers 1...mm (N7..N)
c     initial indices - all to the first center:
      do 41 iex=1,nn
41    si(iex)=1
5     DO  1966 I=1,mm
1966  STATE(I)=0
      do 4 iex=1,nn
4     state(si(iex))=state(si(iex))+1
c     check that the state is allowed:
      do 42 I=1,mm
42    if(STATE(I).gt.1)goto 12

      MS=0
      ML=0
      do 11 I=1,mm
      spin=1-2*mod(I,2)
      orbit=-il+(I-1)/2
      MS=MS+spin *STATE(I)
11    ML=ML+orbit*STATE(I)
      if(iwr.eq.2)then
       if(mod(MS,2).eq.0)then
        write(6,6002)MS/2,ML,(STATE(I),I=1,mm)
6002    format(' MS =',i5,' ML =',i4,2x,30i1)
       else
        write(6,6012)MS,ML,(STATE(I),I=1,mm)
6012    format(' MS =',i3,'/2 ML =',i4,2x,30i1)
       endif
      endif
      i=ML+imax+1
      j=(MS+smax+2)/2
      tab(i,j)=tab(i,j)+1
      NHBIG=NHBIG+1
c
c     find index to be changed
12    do 8 ic=nn,1,-1
8     if(si(ic).lt.mm)goto 9
      goto 3
9     do 10 i=ic+1,nn
10    si(i)=si(ic)+1
      si(ic)=si(ic)+1
      goto 5
c
3     if(iwr.eq.2)then
       write(6,6003)NHBIG
6003   format(i6,' states',/,' Table:')
       do 13 i=1,m
13     write(6,6004)-imax+i-1,(tab(i,j),j=1,n)
6004   format(i3,1x,50i2)
      endif

c     extract S^L states:
      is=0
      ist=0
      do 14 S=smax,smin,-2
      do 14 L=imax,0,-1
c     fill SL table:
      do 15 i=1,m
      do 15 j=1,n
15    lstab(i,j)=0
      do 16 ML=-L,L
      do 16 MS=-S,S,2
      i=ML+imax+1
      j=(MS+smax+2)/2
16    lstab(i,j)=lstab(i,j)+1
      if(iwr.eq.2)then
       write(6,*)'LS Table:'
       do 20 i=1,m
20     write(6,6004)-imax+i-1,(lstab(i,j),j=1,n)
      endif

      ns=0
      ic=0
c     write total to a temporary one: 
172   do 18 i=1,m
      do 18 j=1,n
18    mtab(i,j)=tab(i,j)
c     try to subtract lstab:      
      do 17 i=1,m
      do 17 j=1,n
      if(lstab(i,j).gt.0)then
       if(mtab(i,j).gt.0)then
        mtab(i,j)=mtab(i,j)-1
       else
        ic=1
        goto 171
       endif
      endif
17    continue
171   if(ic.eq.0)then
      if(iwr.eq.2)write(6,*)'can'
c      state could be subtracted, update tab and try again:
       do 19 i=1,m
       do 19 j=1,n
19     tab(i,j)=mtab(i,j)
       if(iwr.eq.2)then
        write(6,*)'Updated Table:'
        do 21 i=1,m
21      write(6,6004)-imax+i-1,(tab(i,j),j=1,n)
       endif
       ns=ns+1
       ist=ist+1
       goto 172
      else
       if(iwr.eq.2)write(6,*)'cannot'
      endif
      if(ns.gt.0)then
       is=is+1
       write(6,6007)is,S+1,sy(L+1),ns
6007   format(i5,i3,'^',a1,i4,' x')
      endif
14    continue
      write(6,6008)is,ist,NHBIG
6008  format(/,i5,' SL labels',/,i5,' SL states',/,
     1i5,' configurations')

      return
      end
C
      subroutine U44(nn,ll,S,L,Sp,Lp,k,iwr)
c     (j1 j2 j||U(k)||j1 j2p jp)
c     =(-1)^(j1+k-j2-jp)(j2||U(k)||j2p)sqrt((2j+1)*(2jp+1)))
c     W(j2 j j2p jp;j1 k)
c     (Racah II 44b)
c
c     j1=j2=j2p=ll, j=L,jp=Lp:
c     (ll ll L||U(k)||ll  ll  Lp)
c     =(-1)^(   k   -Lp)(j2||U(k)||j2p)sqrt((2L+1)*(2Lp+1)))
c     W(j2 j j2p jp;j1 k)
      implicit none
      integer*4 nn,N,Np,Nm,Npm,i,j,isol,isolp,iwr
      real*8 S,Sp,L,ll,Lp,u,k,u44e,l1,s1,l1p,s1p,u44f
      real*8,allocatable::P(:,:),ls(:),ss(:),Pp(:,:),lsp(:),ssp(:)
      real*8 z0,tol
      data z0,tol/0.0d0,1.0d-10/
      if(nn.eq.1)then
       write(6,*)' n=1 not implemented'
       return
      endif
      if(S.ne.Sp)then
       u=z0
      else
       if(nn.eq.2)then
        u=dble(nn)*u44e(ll,L,Lp,k)
       else
c       do CFPs (l^n) for SL and SpLp:
        call e11(nn,ll,S,L)
        call e11(nn,ll,Sp,Lp)
c       read in CFPs for SL:
        N=0
        Nm=0
        call open_cfp(ll,nn,S,L)
        read(11,*,end=992,err=992)N,Nm
992     close(11)
        allocate(P(N,N),ls(N),ss(N))
        call read_CFP(N,Nm,ll,L,S,P,ls,ss,nn)
c       read in CFPs for SpLp:
        Np=0
        Npm=0
        call open_cfp(ll,nn,Sp,Lp)
        read(11,*,end=991,err=991)Np,Npm
991     close(11)
        allocate(Pp(Np,Np),lsp(Np),ssp(Np))
        call read_CFP(Np,Npm,ll,Lp,Sp,Pp,lsp,ssp,nn)

        do 1 isol=1,Nm
        if(iwr.eq.1)write(6,601)isol
 601    format(' SL  state ',i4)
        do 1 isolp=1,Npm
        if(iwr.eq.1)write(6,602)isolp
 602    format('              SpLp  state ',i4)
        u=z0
        do 2 i=1,N
        s1=ss(i)
        l1=ls(i)
        do 2 j=1,Np
        s1p=ssp(j)
        l1p=lsp(j)
        if(iwr.eq.1)
     1  write(6,6005)i, s1, l1, j, s1p, l1p,P(isol,i),Pp(isolp,j),
     1  u44f(l1,l1p,ll,L,Lp,k)
6005    format(i2,2f4.1,i2,2f4.1,3f8.4)
2       if(dabs(s1-s1p).lt.tol.and.dabs(l1-l1p).lt.tol)
     1  u=u+P(isol,i)*Pp(isolp,j)*u44f(l1,l1p,ll,L,Lp,k)
        u=dble(nn)*u
1       write(6,600)u
        return
       endif
      endif
      write(6,600)u
 600  format(' u = ',f12.6)
      return
      end

      subroutine V44(nn,ll,S,L,Sp,Lp,k,iwr)
c     (j1 j2 j||U(k)||j1 j2p jp)
c     =(-1)^(j1+k-j2-jp)(j2||U(k)||j2p)sqrt((2j+1)*(2jp+1)))
c     W(j2 j j2p jp;j1 k)
c     (Racah II 44b)
c
c     j1=j2=j2p=ll, j=L,jp=Lp:
c     (ll ll L||U(k)||ll  ll  Lp)
c     =(-1)^(   k   -Lp)(j2||U(k)||j2p)sqrt((2L+1)*(2Lp+1)))
c     W(j2 j j2p jp;j1 k)
      implicit none
      integer*4 nn,N,Np,Nm,Npm,i,j,isol,isolp,iwr
      real*8 S,Sp,L,ll,Lp,u,k,u44e,l1,s1,l1p,s1p,u44f
      real*8,allocatable::P(:,:),ls(:),ss(:),Pp(:,:),lsp(:),ssp(:)
      real*8 z0,half,one,tol
      data z0,half,one,tol/0.0d0,0.5d0,1.0d0,1.0d-10/
      if(nn.eq.1)then
       write(6,*)' n=1 not implemented'
       return
      endif
      if(nn.eq.2)then
       u=dble(nn)*u44e(ll,L,Lp,k)*u44e(half,S,Sp,one)
      else
c      do CFPs (l^n) for SL and SpLp:
       call e11(nn,ll,S,L)
       call e11(nn,ll,Sp,Lp)
c      read in CFPs for SL:
       N=0
       Nm=0
       call open_cfp(ll,nn,S,L)
       read(11,*,end=992,err=992)N,Nm
992    close(11)
       allocate(P(N,N),ls(N),ss(N))
       call read_CFP(N,Nm,ll,L,S,P,ls,ss,nn)
c      read in CFPs for SpLp:
       Np=0
       Npm=0
       call open_cfp(ll,nn,Sp,Lp)
       read(11,*,end=991,err=991)Np,Npm
991    close(11)
       allocate(Pp(Np,Np),lsp(Np),ssp(Np))
       call read_CFP(Np,Npm,ll,Lp,Sp,Pp,lsp,ssp,nn)

       do 1 isol=1,Nm
       if(iwr.eq.1)write(6,601)isol
 601   format(' SL  state ',i4)
       do 1 isolp=1,Npm
       if(iwr.eq.1)write(6,602)isolp
 602   format('              SpLp  state ',i4)
       u=z0
       do 2 i=1,N
       s1=ss(i)
       l1=ls(i)
       do 2 j=1,Np
       s1p=ssp(j)
       l1p=lsp(j)
       if(iwr.eq.1)
     1 write(6,6005)i, s1, l1, j, s1p, l1p,P(isol,i),Pp(isolp,j),
     1 u44f(l1,l1p,ll,L,Lp,k)
6005   format(i2,2f4.1,i2,2f4.1,3f8.4)
2      if(dabs(s1-s1p).lt.tol.and.dabs(l1-l1p).lt.tol)
     1 u=u+P(isol,i)*Pp(isolp,j)*u44f(l1,l1p,ll,L,Lp,k)
c                          (L1 ll L||U(k)||L1p  ll  Lp)
       u=dble(nn)*u
1      write(6,600)u
       return
      endif
      write(6,600)u
 600  format(' v = ',f12.6)
      return
      end

      function u44e(ll,L,Lp,k)
c     (Racah II 44b)
c     (j1 j2 j||U(k)||j1 j2p jp)
c     (ll ll L||U(k)||ll  ll  Lp)
      implicit none
      real*8 u44e,si,ll,L,Lp,k,bw
      real*8 one,two
      data one,two/1.0d0,2.0d0/
      if(mod(NINT(k)-NINT(Lp),2).eq.0)then
       si=one
      else
       si=-one
      endif
c     (j2||U(k)||j2p)=delta(j2,j2p)=delta(ll,ll)=1="one" in the next line
      u44e=si*one*dsqrt((two*L+one)*(two*Lp+one))*bw(ll,L,ll,Lp,ll,k)
      return
      end

      function u44f(L1,L1p,ll,L,Lp,k)
c     (Racah II 44b)
c     (j1 j2 j||U(k)||j1  j2p jp)
c     (L1 ll L||U(k)||L1p  ll  Lp)
      implicit none
      real*8 u44f,si,ll,L,Lp,k,L1,L1p,bw
      real*8 z0,one,two,tol
      data z0,one,two,tol/0.0d0,1.0d0,2.0d0,1.0d-10/
      if(dabs(L1-L1p).gt.tol)then
       u44f=z0
      else
       if(mod(NINT(L1+k-ll-Lp),2).eq.0)then
        si=one
       else
        si=-one
       endif
       u44f=si*dsqrt((two*L+one)*(two*Lp+one))*bw(ll,L,ll,Lp,L1,k)
      endif
      return
      end

      subroutine e11(nn,ll,S,L)
      implicit none
c     calculate coefficients of fractional parentage 
c     (ll^(nn-1)(apSpLp)ll SL]|ll^nn aSL) based on
c     III: Racah, G. Phys. Rev. 63, 1943, 367-382. Eq. (11)
      integer*4 nn,tn,N,NP,dstates
      real*8 ll,S,L,Sp,Lp,spmin,spmax,lpmin,lpmax,lmax
      real*8,allocatable::ls(:),ss(:),P(:,:)
      real*8 z0,half,one,two,tol
      data z0,half,one,two,tol/0.0d0,0.5d0,1.0d0,2.0d0,1.0d-10/

c     nn=1   ... CFP not relevant
c     nn=2   ... CFP = 1, but some state not allowed (S+L odd)
c     write all possible L S:
      Sp=-one
101   Sp=Sp+one
      Lp=-one
102   Lp=Lp+one
      if(mod(NINT(Lp+Sp),2).eq.0)then
       N=1
       NP=1
       allocate(ss(1),ls(1),P(1,1))
       ss(1)=half
       ls(1)=ll
       P(1,1)=one
       call write_CFP(1,N,N,NP,ll,Lp,Sp,P,ls,ss,2,0)
       deallocate(ss,ls,P)
      endif
      if(Lp.lt.two*ll-tol)goto 102
      if(Sp.lt.one-tol)goto 101

      if(nn.lt.3)then
       return
      endif

c     for (ll^n-1 (LpSp) ll (LS)]| ll^n SL)
c     need to know all under
c     for (ll^tn-1 (LppSpp) ll (LpSp)]| ll^tn SpLp)
      do 1 tn=3,nn-1
      if(mod(tn,2).ne.0)then
       spmin=half
      else
       spmin=z0
      endif
      lpmin=z0
      spmax=dble(tn)*half
      lpmax=lmax(ll,tn)
      Sp=spmin-one
103   Sp=Sp+one
      Lp=lpmin-one
104   Lp=Lp+one
      if(dstates(tn,ll,0,Sp,Lp).gt.0)call e11_a(tn,ll,Sp,Lp)
      if(Lp.lt.lpmax-tol)goto 104
      if(Sp.lt.spmax-tol)goto 103
1     continue

c     finally, get the required ones:
      call e11_a(tn,ll,S,L)
   
      return
      end

      subroutine e11_a(tn,ll,S,L)
c     III: Racah, G. Phys. Rev. 63, 1943, 367-382. Eq. (11)
c     calculate and save (ll^nn-1 (apLpSp) ll (LS)]| ll^nn SL)
      implicit none 
      integer*4 N,M,IER,i,j,IP,JP,NP,Nm,NPm,nn,tn,nmax,izl,ii,k,M0,
     1npos,dstates,code,tcd,nd
      real*8 ll,S,L,sstart,send,lstart,lend,Sp,Lp,Spp,Lpp,f6g,anorm,
     1getdet,SIJ,tem,sppp,lppp,suml,lmax,lij,pnorm
      real*8,allocatable::a(:,:),P(:,:),e(:,:),ai(:,:),ls(:),ss(:),
     1C(:),Pm(:,:),lsm(:),ssm(:),at(:,:),ati(:,:),et(:,:),v(:),
     2Ptem(:)
      logical lt
      integer*4,allocatable::cd(:)
      real*8 z0,half,one,two,tol
      data z0,half,one,two,tol/0.0d0,0.5d0,1.0d0,2.0d0,1.0d-10/
      character*1 sy(17)
      data  sy/'S','P','D','F','G','H','I','K','L','M','N','O','Q','R',
     1         'T','U','V'/

      nn=tn
      write(6,609)nint(ll),nn-1,nint(ll),nn,S,NINT(L),NINT(two*S+one),
     1sy(NINT(L)+1)
609   format(/,'** setting equations for',
     1i2,'^',i1,'|]',i1,'^',i1,' (',f3.1,' ',i2,' = ',i1,'^',
     2a1,') **',/)

      if(L.gt.lmax(ll,nn))then
       write(6,*)' L > Lmax, no solution'
       return
      endif

c     if on disk, do not do anything:
      call open_cfp(ll,nn,S,L)
      N=0
      Nm=0
      read(11,*,end=992,err=992)N,Nm
992   close(11)
      if(N.gt.0)then
       write(6,*)'Found on disk'
       return
      endif

c     find how many allowed:
      npos=dstates(nn,ll,0,S,L)
      if(npos.lt.1)then
       write(6,*)'Impossible combination'
       return
      endif

      sstart=dabs(S-half)
      send=S+half
      lstart=dabs(L-ll)
      lend=min(L+ll,lmax(ll,nn-1))

c     how many coefficients will we have (alphap,Sp,Lp):
      N=0
      Sp=Sstart-one
101   Sp=Sp+one
      Lp=lstart-one
102   Lp=Lp+one
c     read  (ll^nn-2 (appLppSpp) ll (LpSp)]| ll^nn-1 apSpLp)
      call open_cfp(ll,nn-1,Sp,Lp)
      Nm=0
      Npm=0
      read(11,*,end=991,err=991)Nm,Npm
991   close(11)
      if(Nm.gt.0)then
       N=N+Npm
       write(6,601)Npm,Sp,Lp,nn-1, NINT(two*Sp+one),sy(NINT(Lp)+1)
601    format(i10,' CFPs for Sp Lp nn: ',2f4.1,i2,i3,'^',a1)
      endif
      if(lp.lt.lend-tol)goto 102
      if(sp.lt.send-tol)goto 101

c     set maximum M:
      M=0
      Sppp=-one
107   Sppp=sppp+one
      Lppp=-one
108   Lppp=Lppp+one
c     block odd (Lppp+Sppp) states:
      if(mod(nint(Sppp+Lppp),2).ne.0)M=M+1
      if(lppp.lt.two*ll-tol)goto 108
      if(sppp.lt.one-tol)goto 107

      write(6,*)'N M: ',N,M

      nmax=max(M,N)
      allocate(a(nmax,nmax),P(nmax,nmax),ai(nmax,nmax),
     1e(nmax,2*nmax),ls(nmax),ss(nmax),C(nmax),cd(nmax))
      do 21 i=1,nmax
      do 21 j=1,nmax
21    A(i,j)=0.0d0

      write(6,604)
604   format('     f6g      Pm M N  i j l  S   L Sp Lp',
     1' Spp Lpp Sppp Lppp')

      N=0
c     (Lp Sp):
      Sp=Sstart-one
1     Sp=Sp+one
      Lp=lstart-one
2     Lp=Lp+one
c     read  (ll^nn-2 (appLppSpp) ll (LpSp)]| ll^nn-1 apSpLp)
      call open_cfp(ll,nn-1,Sp,Lp)
      Nm=0
      Npm=0
      read(11,*,end=99,err=99)Nm,Npm,tcd
99    close(11)
      if(Nm.gt.0)then
       allocate(Pm(Nm,Nm),lsm(Nm),ssm(Nm))
       call read_CFP(Nm,NPm,ll,Lp,Sp,Pm,lsm,ssm,nn-1)

c      (alpfap):       
       do 111 j=1,Npm
       N=N+1
       ss(N)=Sp
       ls(N)=Lp
       cd(N)=tcd

       M=0
       sppp=-one
7      sppp=sppp+one
       lppp=-one
8      lppp=lppp+one
       if(mod(nint(Sppp+Lppp),2).ne.0)then
        M=M+1
c       sum (alphapp):
c       sum (Spp Lpp):
        do 20 i=1,Nm
        Spp=ssm(i)
        Lpp=lsm(i)
        if(lt(Sppp,Spp,S).and.lt(Spp,half,Sp).and.lt(Sp,half,S).and.
     1     lt(Lppp,Lpp,L).and.lt(Lpp,ll  ,Lp).and.lt(Lp,ll  ,L))then
         write(6,603)f6g(Spp,half,half,Lpp,ll,ll,Sp,Lp,S,L,Sppp,Lppp),
     1   Pm(j,i),M,N,i,j,NINT(ll),S,NINT(L),Sp,NINT(Lp),Spp,NINT(Lpp),
     2   Sppp,NINT(Lppp)
603      format(2f8.4,2i2,1x,3i2,1x,4(f4.1,i2))
         a(M,N)=a(M,N)
     1   +f6g(Spp,half,half,Lpp,ll,ll,Sp,Lp,S,L,Sppp,Lppp)*Pm(j,i)
        endif
20      continue
       endif
       if(lppp.lt.two*ll-tol)goto 8
       if(sppp.lt.one-tol)goto 7

111    continue

       deallocate(Pm,lsm,ssm)
      endif
      if(lp.lt.lend-tol)goto 2
      if(sp.lt.send-tol)goto 1

      call wmtrio(6,a,nmax,M,N)
      
c     delete zero lines
      izl=0
      i=0
22    i=i+1
      suml=z0
      do 23 j=1,N
23    suml=suml+dabs(A(i,j))
      if(suml.lt.tol)then
       izl=izl+1
       do 24 ii=i,M-1
       do 24 j=1,N
24     A(ii,j)=A(ii+1,j)
       M=M-1
       i=i-1
      endif 
      if(i.lt.M)goto 22
      if(izl.gt.0) write(6,*)izl,' zero lines deleted'

      write(6,*)'m-matrix:'
      call wmtrio(6,a,nmax,M,N)

      if(M.gt.N.and.N.gt.0)then
       write(6,*)' M > N'
       M0=M
       M=N
201    tem=getdet(a,nmax,M)
       write(6,*)M,' x ',M,' determinant: ',tem
       if(dabs(tem).gt.tol)then
        allocate(at(M,M),ati(M,M),et(M,2*M),v(M))
        do 202 i=1,M
        do 202 j=1,M
202     at(i,j)=a(j,i)
        call INV(M,At,AtI,M,tol,et,IER)
c       check that remaining lines can be obtained as lin .
c       combination:
        do 203 i=M+1,M0
c       line i:
        do 205 k=1,M
        v(k)=z0
        do 205 j=1,M
205     v(k)=v(k)+Ati(k,j)*a(i,j)
        write(6,*)'i:',i
        write(6,*)'v:'
        write(6,*)(v(k),k=1,M)
c       is it for all N?:
        do 203 j=1,N
        lij=z0
        do 204 k=1,M
204     lij=lij+a(k,j)*v(k)
        if(dabs(lij-a(i,j)).gt.tol)then
         write(6,*)i,j,lij,a(i,j),' no lin dep, no solution'
         return
        endif
203     continue
        deallocate(at,ati,et,v)
        write(6,*)' additional lines linear dependent and forgotten'
       else
        if(M.gt.1)then
         M=M-1
         goto 201
        else
         write(6,*)' M > N, no solution'
         return
        endif
       endif
      else
c      complete to square matrix:
       do 19 i=M+1,N
       do 19 j=1,N
19     a(i,j)=a(M,j)
       write(6,*)' m completed'
      endif
     
      tem=getdet(a,nmax,N)
      write(6,*)'determinant: ',tem

      if(dabs(tem).gt.tol)then
       write(6,*)'No solution'
       return
      endif

      M=N
      NP=0
c     NP ... number of possible solutions

12    M=M-1
      NP=NP+1
      do 6 IP=1,NP
      DO 6 i =N-NP+1,N
6     P(IP,i)=z0
      do 13 IP=1,NP
13    P(IP,N-IP+1)=one
      if(dabs(getdet(a,nmax,M)).lt.tol.and.M.gt.0)goto 12

      do 25 j=1,N
25    write(6,*)'N:',j,' Sp:',ss(j),'LP:',NINT(ls(j)),' code: ',cd(j)
  
      if(M.gt.0)then
       call INV(nmax,A,AI,M,tol,e,IER)
       if(IER.eq.0)then
        if(NP.eq.1)then
         write(6,*)NP,' solution possible'
        else
         write(6,*)NP,' solutions possible'
        endif
        if(NP.gt.npos)then
         write(6,*)'warning: only ',npos,' allowed'
c        write(6,*)NP-npos,' forgotten'
c        NP=npos
         code=npos
        else
         code=0
        endif
        do 3 IP=1,NP
        DO 15 i=1,M
        C(i)=z0
        do 15 j=M+1,N
15      C(i)=C(i)-A(i,j)*P(IP,j)
        do 3 i=1,M
        P(IP,i)=z0
        do 3 j=1,M
3       P(IP,i)=P(IP,i)+AI(i,j)*C(j)
       else
        write(6,*)'no inversion'
        return
       endif
      endif
    
c     orthonormalize:
      do 11 IP=1,NP
      SIJ=z0
      do 18 JP=1,IP-1
      do 17 i=1,N
17    SIJ=SIJ+P(IP,i)*P(JP,i)
      do 18 i=1,N
18    P(IP,i)=P(IP,i)-SIJ*P(JP,i)
      anorm=z0
      do 10 i=1,N
10    anorm=anorm+P(IP,i)**2
      do 11 i=1,N
11    P(IP,i)=P(IP,i)/dsqrt(anorm)

      do 26 j=1,N
      Sp=ss(j)
      Lp=ls(j)
      if(cd(j).gt.0)then
       nd=1
       do 27 jp=j+1,N
27     if(ss(jp).eq.Sp.and.ls(jp).eq.Lp)nd=nd+1
       if(nd.eq.2.and.code.eq.0.and.NP.eq.1)then
        write(6,*)'Repairing substate:'
        write(6,*)'N:',j,' Sp:',Sp,'Lp:',NINT(Lp),' code: ',cd(j)
        call open_cfp(ll,nn-1,Sp,Lp)
        Nm=0
        Npm=0
        read(11,*,end=991,err=991)Nm,Npm
        close(11)
        allocate(Pm(Nm,Nm),lsm(Nm),ssm(Nm),Ptem(Nm))
        call read_CFP(Nm,NPm,ll,Lp,Sp,Pm,lsm,ssm,nn-1)
        pnorm=z0
        do 29 i=1,nd
        write(6,*)'jp P(1,jp)',j+i-1,P(1,j+i-1)
29      pnorm=pnorm+P(1,j+i-1)**2
        pnorm=dsqrt(pnorm)
        write(6,*)'pnorm: ',pnorm
        do 28 jp=1,Nm
        Ptem(jp)=z0
        do 28 i=1,nd
28      Ptem(jp)=Ptem(jp)+P(1,j+i-1)*Pm(i,jp)/pnorm
        do 30 jp=1,Nm
30      Pm(1,jp)=Ptem(jp)
        Npm=1
        call write_CFP(1,Nm,Nm,Npm,ll,Lp,Sp,Pm,lsm,ssm,nn-1,0)
        write(6,*)'Repairing state:'
        if(P(1,j).lt.z0)then
         P(1,j)=-pnorm
        else
         P(1,j)= pnorm
        endif
        do 31 jp=j+1,N-1
        P(1,jp)=P(1,jp+1)
        ss(jp)=ss(jp+1)
31      ls(jp)=ls(jp+1)
        deallocate(Pm,lsm,ssm,Ptem)
        N=N-1
       endif
      endif
26    continue

      call write_CFP(1,nmax,N,NP,ll,L,S,P,ls,ss,nn,code)
      close(11)

      return

      end

      subroutine e9(ll,S,L)
      implicit none
c     III: Racah, G. Phys. Rev. 63, 1943, 367-382. Eq. (9)
      integer*4 N,M,ic,Mmax,IER,i,j,IP,JP,NP,Ndim
      real*8 ll,S,L,sstart,send,lstart,lend,Sp,Lp,Spp,Lpp,f6,anorm,
     1getdet,SIJ,tem
      real*8,allocatable::a(:,:),P(:,:),e(:,:),ai(:,:),ls(:),ss(:),
     1C(:)
      real*8 z0,half,one,two,tol
      data z0,half,one,two,tol/0.0d0,0.5d0,1.0d0,2.0d0,1.0d-10/
      character*1 sy(17)
      data  sy/'S','P','D','F','G','H','I','K','L','M','N','O','Q','R',
     1         'T','U','V'/

      write(6,609)nint(ll),3,nint(ll),2
609   format('** setting equations for',i2,'^',i1,'|]',i1,'^',i1,' **')
      sstart=dabs(S-half)
      send=S+half
      lstart=dabs(L-ll)
      lend=L+ll
      
      do 1000 ic=1,2
      Mmax=0
      N=0
      if(ic.eq.2)write(6,*)' Sp  Lp:'
      Sp=Sstart-one
1     Sp=Sp+one
      Lp=lstart-one
2     Lp=Lp+1
      if(mod(nint(Sp+Lp),2).eq.0)then
      if(Sp.le.one+tol.and.Lp.le.two*ll+tol)then
c      this is allowed (aantisymmetric) state, and we want its CFPs:
       N=N+1
       if(ic.eq.2)then
        write(6,600)Sp,NINT(Lp),NINT(two*Sp+one),sy(NINT(Lp)+1)
600     format(f3.1,i3,' (',i1,'^',a1,')',/,12x,'Spp  Lpp:')
        ss(N)=Sp
        ls(N)=Lp
       endif
       M=0
       spp=sstart-one
7      spp=spp+one
       lpp=lstart-one
8      lpp=lpp+1
       if(mod(nint(Spp+Lpp),2).ne.0)then
       if(Spp.le.one+tol.and.Lpp.le.two*ll+tol)then
c       block this state:
        M=M+1
        if(ic.eq.2)then
         a(M,N)=f6(ll,Sp,Lp,S,L,Spp,Lpp)
         write(6,602)Spp,NINT(Lpp),NINT(two*Spp)+1,
     1   sy(NINT(Lpp)+1),f6(ll,Sp,Lp,S,L,Spp,Lpp)
602      format(17x,f3.1,i3,' (',i1,'^',a1,') ',f10.4)
        endif
       endif
       endif
       if(lpp.lt.lend-tol)goto 8
       if(spp.lt.send-tol)goto 7
       if(M.gt.Mmax)Mmax=M
      endif
      endif
      if(lp.lt.lend-tol)goto 2
      if(sp.lt.send-tol)goto 1
      if(ic.eq.1)then
       write(6,*)' number of CFPs ',N
       write(6,*)' number of equations ',Mmax
       Ndim=max(Mmax,N)
       N=Ndim
       allocate(a(N,N),P(N,N),ai(N,N),e(N,2*N),ls(N),ss(N),C(N))
      else
       if(Mmax.gt.N)then
        write(6,*)' M > N, no solution'
        return
       else
c       complete to square matrix:
        do 19 i=Mmax+1,Ndim
        do 19 j=1,N
19      a(i,j)=a(Mmax,j)
       endif
      endif
1000  continue
      
      write(6,*)'m-matrix:'
      call wmtrio(6,a,N,N,N)

      tem=getdet(a,N,N)
      write(6,*)'determinant: ',tem

      if(dabs(tem).gt.tol)then
       write(6,*)'No solution'
       return
      endif

      M=N
      NP=0
c     NP ... number of possible solutions

12    M=M-1
      NP=NP+1
      do 6 IP=1,NP
      DO 6 i =N-NP+1,N
6     P(IP,i)=z0
      do 13 IP=1,NP
13    P(IP,N-IP+1)=one
      if(dabs(getdet(a,N,M)).lt.tol.and.M.gt.0)goto 12

      if(M.gt.0)then
       call INV(N,A,AI,M,tol,e,IER)
       if(IER.eq.0)then
        write(6,*)NP,' solutions possible'
        do 3 IP=1,NP
        DO 15 i=1,M
        C(i)=z0
        do 15 j=M+1,N
15      C(i)=C(i)-A(i,j)*P(IP,j)
        do 3 i=1,M
        P(IP,i)=z0
        do 3 j=1,M
3       P(IP,i)=P(IP,i)+AI(i,j)*C(j)
       else
        write(6,*)'no inversion'
        return
       endif
      endif
    
c     orthonormalize:
      do 11 IP=1,NP
      SIJ=z0
      do 18 JP=1,IP-1
      do 17 i=1,N
17    SIJ=SIJ+P(IP,i)*P(JP,i)
      do 18 i=1,N
18    P(IP,i)=P(IP,i)-SIJ*P(JP,i)
      anorm=z0
      do 10 i=1,N
10    anorm=anorm+P(IP,i)**2
      do 11 i=1,N
11    P(IP,i)=P(IP,i)/dsqrt(anorm)

      call write_CFP(1,N,N,NP,ll,L,S,P,ls,ss,3,0)

      return

      end

      SUBROUTINE INV(nr,a,ai,n,TOL,e,IERR)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      dimension ai(nr,nr),a(nr,nr),e(nr,2*nr)
C
10000 IERR=0
      DO 1  ii=1,n
      DO 1  jj=1,n
      e(ii,jj)=a(ii,jj)
1     e(ii,jj+n)=0.0D0
      do 13 ii=1,n
13    e(ii,ii+n)=1.0D0
c
      DO 2 ii=1,n-1
      if (ABS(e(ii,ii)).LE.TOL) then
       DO 3 io=ii+1,n
3      if (ABS(e(io,ii)).GT.TOL) goto 11
       IERR=1
       write(6,*)ii
       write(6,*)'tol = ',tol
       tol=tol*0.50d0
       if(tol.gt.1.0d-20)goto 10000
       RETURN
c
11     CONTINUE
       DO 4 kk=1,2*n
       w=e(ii,kk)
       e(ii,kk)=e(io,kk)
4      e(io,kk)=w
      ENDIF
      eii=e(ii,ii)
      DO 5 jj=ii+1,n
      e1=e(jj,ii)/eii
      DO 6 kk=ii+1, 2*n
6     e(jj,kk)=e(jj,kk)-e(ii,kk)*e1
5     e(jj,ii)=0.0D0
2     CONTINUE
c
      DO 7  i2=n,2,-1
      eii=e(i2,i2)
      DO 7 j2=i2-1,1,-1
      e1=e(j2,i2)/eii
      DO 9 kk=1, n
9     e(j2,kk+n)=e(j2,kk+n)-e(i2,kk+n)*e1
7     e(j2,i2)=0.0d0
c
      DO 10 ii=1,n
      ei=1.0d0/e(ii,ii)
      DO 12  jj=1,n
12    ai(ii,jj)=e(ii,jj+n)*ei
10    CONTINUE
c
      RETURN
      END             

      subroutine wmtrio(io,A,n0,n,m)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      dimension A(n0,n0)
      N1=1
1     N3=MIN(N1+4,M)
      WRITE(io,18)(I,I=N1,N3)
18    FORMAT(4X,5I14)
      DO 130 LN=1,N
130   WRITE(io,17)LN,(A(LN,J),J=N1,N3)
      N1=N1+5
      IF(N3.LT.M)GOTO 1
17    FORMAT(I4,5D14.6)
      return
      end

      function f6(ll,Sp,Lp,S,L,Spp,Lpp)
c     (s1l1s2l2(SpLp)s3l3SL|s1l1,s2l2s3l3(SppLpp),SL)
c     Racah III (6)
c     si=0.5,li=ll
      implicit none
      real*8 ll,Sp,Lp,S,L,Spp,Lpp,f6,bw
      real*8 half,one,two
      data half,one,two/0.5d0,1.0d0,2.0d0/
      f6=dsqrt((two*Sp+one)*(two*Spp+one)*(two*Lp+one)*(two*Lpp+one))
     1*bw(half,half,S,half,Sp,Spp)*bw(ll,ll,L,ll,Lp,Lpp)
      return
      end

      function f6g(s1,s2,s3,l1,l2,l3,Sp,Lp,S,L,Spp,Lpp)
c     general (s1l1s2l2(SpLp)s3l3SL|s1l1,s2l2s3l3(SppLpp),SL)
c     Racah III (6)
      implicit none
      real*8 Sp,Lp,S,L,Spp,Lpp,f6g,bw,s1,s2,s3,l1,l2,l3
      real*8 one,two
      data one,two/1.0d0,2.0d0/
      f6g=dsqrt((two*Sp+one)*(two*Spp+one)*(two*Lp+one)*(two*Lpp+one))
     1*bw(s1,s2,S,s3,Sp,Spp)*bw(l1,l2,L,l3,Lp,Lpp)
600   format(2f8.4,6f4.1,1x,6f4.1)
      return
      end

      function sw(a,b,c,d,e,ff)
c     Racah II (36)
      implicit none
      integer*4 a1,a2,a3,a4,a5,a6,a7,z
      real*8 a,b,c,d,e,ff,sw,si,s
      integer*4 nfac0
      parameter (nfac0=100)
      real*8 f0,f
      common/npar/f0,f(nfac0)
      s=0.0d0
      si=-1.0d0
      do 1 z=0,nint(a+b+c+d)+1
      si=-si
      a1=nint(a + b+ c+d)+1-z
      a2=nint(a + b- e)-z
      a3=nint(c + d- e)-z
      a4=nint(a + c-ff)-z
      a5=nint(b + d-ff)-z
      a6=nint(e +ff- a-d)+z
      a7=nint(e +ff- b-c)+z
1     if(a1.ge.0.and.a2.ge.0.and.a3.ge.0.and.a4.ge.0.and.a5.ge.0.and.
     1   a6.ge.0.and.a7.ge.0)s=
     2   s+si*f(a1)/f(a2)/f(a3)/f(a4)/f(a5)/f(z)/f(a6)/f(a7)
      sw=s
      return
      end

      function bw(a,b,c,d,e,ff)
c     Racah II, eq. 36'
      implicit none
      integer*4 a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16
      real*8 a,b,c,d,e,ff,sw,bw
      integer*4 nfac0
      parameter (nfac0=100)
      real*8 f0,f
      common/npar/f0,f(nfac0)
      a1 =nint(a+b-e)
      a2 =nint(a+e-b)
      a3 =nint(b+e-a)
      a4 =nint(c+d-e)
      a5 =nint(c+e-d)
      a6 =nint(d+e-c)

      a7 =nint(a+c-ff)
      a8 =nint(a+ff-c)
      a9 =nint(c+ff-a)
      a10=nint(b+d-ff)
      a11=nint(b+ff-d)
      a12=nint(d+ff-b)

      a13=nint(a+b+ e)+1
      a14=nint(c+d+ e)+1
      a15=nint(a+c+ff)+1
      a16=nint(b+d+ff)+1
      if(a1 .ge.0.and.a2 .ge.0.and.a3 .ge.0.and.a4 .ge.0.and.
     1   a5 .ge.0.and.a6 .ge.0.and.a7 .ge.0.and.a8 .ge.0.and.
     1   a9 .ge.0.and.a10.ge.0.and.a11.ge.0.and.a12.ge.0.and.
     1   a13.ge.0.and.a14.ge.0.and.a15.ge.0.and.a16.ge.0)then
       bw=dsqrt(f(a1)*f(a2)*f(a3)*f(a4)*f(a5)*f(a6)*f(a7)*f(a8)*f(a9)
     1 *f(a10)*f(a11)*f(a12)/f(a13)/f(a14)/f(a15)/f(a16))
     2 *sw(a,b,c,d,e,ff)
      else
       bw=0.0d0
      endif
      return
      end

      Function Fac(n)
      implicit none
      integer*4 n,ii
      real*8 Fac,bb
      bb=1.0d0
      do 1 ii=2 , n
1     bb=bb*dble(ii)
      fac=bb
      return
      end

      Function DFac(n)
c     double factorial = n!!
      implicit none
      integer*4 n,ii
      real*8 DFac,bb
      if(n.eq.-1.or.n.eq.0)then
       bb=1.0d0
      else
       bb=dble(n)
       if(mod(n,2).ne.0)then
        do 1 ii=n-2,1,-2
1       bb=bb*dble(ii)
       else
        do 2 ii=n-2,2,-2
2       bb=bb*dble(ii)
       endif
      endif
      dfac=bb
      return
      end

      subroutine Initfac
      implicit none
      integer*4 nfac0,i
      parameter (nfac0=100)
      real*8 f0,factorial,fac
      common/npar/f0,factorial(nfac0)
      real*8 dm,d0,df,dfac
      common/dpar/dm,d0,df(nfac0)
      do  1 i=0 , nfac0
1     factorial(i)=fac(i)
      do  2 i=-1, nfac0
2     df(i)=dfac(i)
      return
      end

      FUNCTION getdet(A,N0,N)
      IMPLICIT none
      INTEGER*4 N, N0, I, J, K
      REAL*8 A(N0,N0),T, TEMP,L,getdet
      LOGICAL DETEXISTS
      REAL*8,allocatable:: ELEM(:,:)

      allocate(ELEM(N,N))
      DO 1 I=1,N
      DO 1 J=1,N
1     ELEM(I,J)=A(I,J)
        
      DETEXISTS = .TRUE.
      L = 1.0d0
c     CONVERT TO UPPER TRIANGULAR FORM
      DO 3 K = 1, N-1
      IF (DABS(ELEM(K,K)).LE.1.0D-10) THEN
       DETEXISTS = .FALSE.
       DO 4 I = K+1, N
       IF (dabs(ELEM(I,K)).gt.1.0d-10) THEN
        DO 2 J = 1, N
        TEMP = ELEM(I,J)
        ELEM(I,J)= ELEM(K,J)
2       ELEM(K,J) = TEMP
        DETEXISTS = .TRUE.
        L=-L
        EXIT
       ENDIF
4      continue
       IF (DETEXISTS .EQV. .FALSE.) THEN
        getdet = 0.0d0
        RETURN
       ENDIF
      ENDIF
      DO 3 J = K+1, N
      T = ELEM(J,K)/ELEM(K,K)
      DO 3 I = K, N
3     ELEM(J,I) = ELEM(J,I) - T*ELEM(K,I)
     
      !CALCULATE DETERMINANT BY FINDING PRODUCT OF DIAGONAL ELEMENTS
      getdet = L
      DO 5 I = 1, N
5     getdet = getdet * ELEM(I,I)
      END

      subroutine open_cfp(ll,n,S,L)
      implicit none
      real*8 ll,S,L
      integer*4 il,jl,n
      character*10 ss
      il=NINT(ll)
      jl=NINT(L)
      write(ss(1:1),607)il
607   format(i1)
      ss(2:2)='^'
      write(ss(3:3),607)n
      ss(4:4)='_'
      write(ss(5:7),608)S
608   format(f3.1)
      ss(8:8)='_'
      write(ss(9:10),609)jl
609   format(i2)
      if(ss(9:9).eq.' ')ss(9:9)='_'
      open(11,file=ss)
      return
      end

      subroutine write_CFP(ic,n0,N,NP,ll,L,S,P,ls,ss,nn,code)
      implicit none
      integer*4 N,NP,ic,il,jl,IP,i,nn,n0,code
      real*8 ll,L,S,ls(*),ss(*)
      real*8 one,two,P(n0,n0)
      data one,two/1.0d0,2.0d0/
      character*1 sy(17),ssy(17)
c     index     1   2   3   4   5   6   7   8   9  10  11  12  13  14
c              15  16  17
c     L:        0   1   2   3   4   5   6   7   8   9  10  11  12  13
c              14  15  16
      data  sy/'S','P','D','F','G','H','I','K','L','M','N','O','Q','R',
     1         'T','U','V'/
      data ssy/'s','p','d','f','g','h','i','k','l','m','n','o','q','r',
     1         't','u','v'/
c
      call open_cfp(ll,nn,S,L)
      il=NINT(ll)
      jl=NINT(L)
      write(11,1100)N,NP,code
1100  format(3i4)
c     N number of nn-1 basis functions
c     NP number of solutions (additional quantum number) for ll^nn SL
c     code 0 if ok, 1 if solution is incomplete

      if(ic.eq.1)then
      write(6,603)il,nn,S,NINT(L),ssy(il+1),nn,NINT(two*S+one),sy(jl+1)
      write(6,604)
      do 1 IP=1,NP
      write( 6,606)IP
      do 1 i=1,N
1     write( 6,605)i,ssy(il+1),nn-1,ss(i),NINT(ls(i)),
     1NINT(two*ss(i)+one),sy(NINT(ls(i))+1),P(IP,i)
      endif

      write(11,603)il,nn,S,NINT(L),ssy(il+1),nn,NINT(two*S+one),sy(jl+1)
 603  format(/,' Term |l^n,S,L> = |',i1,'^',i1,',',f3.1,',',i1,
     1'> = |l^n M^L> = |',a1,'^',i1,' ',i1,'^',a1,'>:',/)
      write(11,604)
604   format(
     1' |l^n-1 (Sp Lp)=(M^L) l>  x  (l^n SL|[l^n-1 (Sp Lp) l SL)')
      do 9 IP=1,NP
      write(11,606)IP
606   format(/,' Solution ',i2,':',/)
      do 9 i=1,N
9     write(11,605)i,ssy(il+1),nn-1,ss(i),NINT(ls(i)),
     1NINT(two*ss(i)+one),sy(NINT(ls(i))+1),P(IP,i)
605   format(i3,'    |',a1,'^',i1,
     1' (',f3.1,i2,')=(',i1,'^',a1,')',f22.15)

      close(11)
      return
      end

      subroutine read_CFP(N,NP,ll,L,S,P,ls,ss,nn)
      implicit none
      integer*4 N,NP,IP,i,nn,ils
      real*8 ll,L,S,ls(*),ss(*)
      real*8 P(N,N)
      call open_cfp(ll,nn,S,L)
      do 1 i=1,5
1     read(11,*)
      do 9 IP=1,NP
      do 2 i=1,3
2     read(11,*)
      do 9 i=1,N
      read(11,605)ss(i),ils,P(IP,i)
9     ls(i)=dble(ils)
605   format(3x,10x,f3.1,i2,3x,1x,1x,1x,1x,f22.15)
      close(11)
      return
      end

      function lmax(ll,nn)
c     maximum L in ll-shell with nn electrons
      implicit none
      real*8 ll,lc,lmax
      integer*4 nn,m,ie,ic
      ie=0
      lc=0.0d0
      do 1 m=NINT(ll),-NINT(ll),-1
      do 1 ic=1,2
      if(ie.lt.nn)then
       ie=ie+1
       lc=lc+dble(m)
      else
       goto 2
      endif
1     continue
2     lmax=lc
      return
      end

      function lt(a,b,c)
c     checks triangle conditions
      implicit none
      real*8 a,b,c,tol
      logical lt
      data tol/1.0d-10/
      lt=a.le.b+c+tol.and.a.ge.dabs(b-c)-tol.and.
     1   b.le.c+a+tol.and.b.ge.dabs(c-a)-tol.and.
     1   c.le.a+b+tol.and.c.ge.dabs(a-b)-tol
      return
      end