program racahtest
c
c     tests racah-coeff
c
      IMPLICIT REAL*8 (A-H,O-Z)
      implicit integer(i-n)
      COMMON/REDFAK/FAKM(100)

      CALL FACTOR

      Do i=1,10
       Print *,i,'     fakm(i)',FAKM(i)
      enddo

      A1=1.
      A2=1.
      A3=1.
      B1=1.
      B2=1.
      B3=0.

      rac=racah(A1,A2,A3,B1,B2,B3)

      print *,'rac ',rac

      Stop
      End
C
C=====================================================================
C      
      Double Precision FUNCTION RACAH(A1,A2,A3,B1,B2,B3)
C
C     6-J SYMBOLS, RACAH FORMULA,  A.MESSIAH II (1962), APENDIX C
C     FAKM(N) = GAMMA(N)*10**(-N+1)
C
      IMPLICIT REAL*8 (A-H,O-Z)
      implicit integer(i-n)
      COMMON/REDFAK/FAKM(100)
      external delta
C
      RACAH=0.0D0
      DEL1=DELTA(A1,A2,A3)
      IF(DEL1.EQ.0.0D0) RETURN
      DEL2=DELTA(A1,B2,B3)
      IF(DEL2.EQ.0.0D0) RETURN
      DEL3=DELTA(B1,A2,B3)
      IF(DEL3.EQ.0.0D0) RETURN
      DEL4=DELTA(B1,B2,A3)
      IF(DEL4.EQ.0.0D0) RETURN
      K1=A1+A2+A3+0.001D0
      K2=A1+B2+B3+0.001D0
      K3=B1+A2+B3+0.001D0
      K4=B1+B2+A3+0.001D0
      K5=A1+A2+B1+B2+2.001D0
      K6=A2+A3+B2+B3+2.001D0
      K7=A3+A1+B3+B1+2.001D0
      MU=MAX0(K1,K2,K3,K4)+1
      MO=MIN0(K5,K6,K7)-1
      IF(MO.LT.MU) RETURN
      SU=0.0D0
      DO 1 K=MU,MO
      J1=K-K1
      J2=K-K2
      J3=K-K3
      J4=K-K4
      J5=K5-K
      J6=K6-K
      J7=K7-K
1     SU=-SU+FAKM(K+1)/(FAKM(J1)*FAKM(J2)*FAKM(J3)*FAKM(J4)*FAKM(J5)
     1*FAKM(J6)*FAKM(J7))
      SS=-1.0D0
      IF(MOD(MO,2).EQ.1) SS=1.0D0
      RACAH=SS*SU*DSQRT(DEL1*DEL2*DEL3*DEL4)*10.0D0
      IF(DABS(RACAH).LT.1.0D-6)RACAH=0.0D0
      RETURN
      END
C
C=====================================================================
C      
      Double Precision FUNCTION DELTA(A,B,C)
C
C     SUBPROGRAM OF RACAH-ROUTINE
C     FAKM(N) = GAMMA(N)*10**(-N+1)
C
      IMPLICIT REAL*8 (A-H,O-Z)
      implicit integer(i-n)
      COMMON/REDFAK/FAKM(100)
      DELTA=0.0D0
      TEST=A+B+C-DINT(A+B+C+0.001d0)+0.001D0 ! Integer relation
      IF(TEST.GT.0.1D0) RETURN
      J1=A+B-C+100.01D0          ! triangular relation
      J2=B+C-A+100.01D0
      J3=C+A-B+100.01D0
      J4=A+B+C+2.01D0
      J1=J1-99
      J2=J2-99
      J3=J3-99
      IF(MIN0(J1,J2,J3 ).LE.0) RETURN
      DELTA=0.1d0*FAKM(J1)*FAKM(J2)/FAKM(J4)*FAKM(J3)
      RETURN
      END
C
C=====================================================================
C      
      SUBROUTINE FACTOR
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON/REDFAK/FAKM(100)
      FAKM(1)=1.0D0
      DO 1 I=1,40
        AI=dfloat(I)
        FAKM(I+1)=FAKM(I)*0.1D0*AI
1     CONTINUE
      RETURN
      END