      implicit real*8(a-h,o-z)
      character*1 it
      logical half
      write(6,*) ' calculates Clebsch-Gordon coefficients'
      write(6,*) '  for <j1m1j2m2|jm>'
       loop: do
      write(6,*)'enter "I"/"H" - whole or half integral; other to stop'
      read(5,'(a1)') it
      if (it.ne.'I'.and.it.ne.'H') exit
      half=.false.
      if(it.eq.'H') half=.true.
      write(6,*) 'input j1,m1,j2,m2,j,m (x2 if "H")'
      read(5,*)j1,m1,j2,m2,j,m
      xx=clebsch(half,j1,m1,j2,m2,j,m)
      write(6,*) xx
      enddo loop
      stop
      end
       function clebsch(half,j1,m1,j2,m2,j,m)
      implicit real*8 (a-h,o-z)
      logical half
c  returns clebsch-gordon coeffient
c condon & shortly convention -- using racah equation quoted by rotenberg
      clebsch=0.d0
      if ((-m1-m2+m).ne.0) return
      if(iabs(m1).gt.j1) return
      if(iabs(m2).gt.j2) return
      if(iabs(m).gt.j) return
      jj=j1+j2+j
      if(half.and.(jj/2)*2.ne.jj) return
      jj1=jj-2*j1
      jj2=jj-2*j2
      jj3=jj-2*j
      if(jj1.lt.0) return
      if(jj2.lt.0) return
      if(jj3.lt.0) return
      if (half) then
       xx=j+1
       jj1=jj1/2
       jj2=jj2/2
       jj3=jj3/2
       jj=jj/2
       xx=(xx*factrl(jj1)*factrl(jj2)/(factrl(jj3)*factrl(jj+1)))
        if(m1.ge.0) xx=xx*rfctrl((j1+m1)/2,m1)
        if(m1.lt.0) xx=xx/rfctrl((j1-m1)/2,-m1)
        if(m2.ge.0) xx=xx/rfctrl((j2+m2)/2,m2)
        if(m2.lt.0) xx=xx*rfctrl((j2-m2)/2,-m2)
        if(m.ge.0) xx=xx*rfctrl((j+m)/2,m)
        if(m.lt.0) xx=xx/rfctrl((j-m)/2,-m)
        xx=dsqrt(xx)
       else
       xx=2*j+1
       xx=(xx*factrl(jj1)*factrl(jj2)/(factrl(jj3)*factrl(jj+1)))
        if(m1.ge.0) xx=xx*rfctrl(j1+m1,2*m1)
        if(m1.lt.0) xx=xx/rfctrl(j1-m1,-2*m1)
        if(m2.ge.0) xx=xx/rfctrl(j2+m2,2*m2)
        if(m2.lt.0) xx=xx*rfctrl(j2-m2,-2*m2)
        if(m.ge.0) xx=xx*rfctrl(j+m,2*m)
        if(m.lt.0) xx=xx/rfctrl(j-m,-2*m)
        xx=dsqrt(xx)
      endif
       clebsch=xx
      ii1=j-j2+m1
      ii2=j-j1-m2
      nn1=j1-m1
      nn2=j2+m2
      jm=j-m
      sum=0.d0
      sgn=1.d0
       if(half) then
        ii1=ii1/2
        ii2=ii2/2
        nn1=nn1/2
        nn2=nn2/2
        jm=jm/2
       endif
      k1=-min0(0,ii1,ii2)
      k2=min0(jj3,nn1,nn2)
      if(k1.gt.0) sgn=(-1)**k1
      do k=k1,k2
      arg=sgn*(rfctrl(nn1,k)/factrl(k))*factrl(jm)
      arg=arg*(rfctrl(nn2,k)/factrl(ii1+k))
      arg=arg*(rfctrl(jj3,k)/factrl(ii2+k))
      sum=sum+arg
      sgn=-sgn
      enddo
      clebsch=clebsch*sum
      return
      end
      function factrl(k)
      implicit real*8 (a-h,o-z)
c   function to return factorial of k
c     k must be ge 0
      if(k.lt.0) then
      write(6,'("***error in factrl**** k = ",i10)') k
       stop
       endif
      factrl=1.d0
      if(k.eq.0) return
      do  i=1,k
      factrl=factrl*i
      enddo
      return
      end
      function rfctrl(k,i)
      implicit real*8 (a-h,o-z)
c  calculate (k!)/((k-i)!)
c     k, k-i must be ge 0
      if(k.lt.0.or.k-i.lt.0) then
      write(6,'("***error in rfctrl**** k , i = ",2i10)') k,i
       stop
       endif
      rfctrl=1.d0
      if(i.eq.0) return
       in=k-i+1
       it=k
       do  l=in,it
       rfctrl=rfctrl*l
       enddo
       return
       end
