csinpm
c234567
      program s3j
      implicit real*8 (a-h,o-z)
      call consfa
      call consea

      read(*,*)j1,j2,j3,m1,m2,m3
      print*,sex3j(j1,j2,j3,m1,m2,m3)
      end

     
*----------------------------------------------------------------------*
*  CLEBS-GORDAN COEF. <L1,L2,M1,M2,L,M>=CG(2*L1,2*L2,2*M1,2*M2,2*J,2*M)*
*----------------------------------------------------------------------*
      real*8 FUNCTION CGC(J1,M1,j2,M2,J,M)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      CNO=DSQRT(DFLOAT(J+1))
c      print*,j1,m1,j2,m2,j,m
      cno=cno*(-1)**(iabs(J1-J2+M)/2)
      CNO2=Sex3J(J1,J2,J,M1,M2,-M)
      CGC=CNO*CNO2
      RETURN
      END 
      real*8 FUNCTION CGC1(J1,M1,j2,M2,J,M)
      IMPLICIT real*8 (A-H,O-Z)
c      print*,j1,m1,j2,m2,j,m
      if((m1+m2-m).ne.0)then
      cgc1=0.d0
      return
      endif
      if((j.lt.iabs(j1-j2)).or.(j.gt.(j1+j2)))then
      cgc1=0.d0
      return
      endif
      if((iabs(m1).gt.j1).or.(iabs(m).gt.j))then
      cgc1=0.d0
      return
      endif 
      if((j-j1).eq.1)then
      if(m2.eq.-1)cgc1=dsqrt(dfloat(j1-m)/dfloat(j1+1)/2.d0)
      if(m2.eq.1)cgc1=dsqrt(dfloat(j1+m)/dfloat(j1+1)/2.d0)
c      print*,'caso 1',m2,cgc1
      return
      endif
      if((j-j1).eq.-1)then
      if(m2.eq.-1)cgc1=dsqrt(dfloat(j1+m)/dfloat(j1+1)/2.d0)
      if(m2.eq.1)cgc1=-dsqrt(dfloat(j1-m)/dfloat(j1+1)/2.d0)
c      print*,'caso2',m2,cgc1
      return
      endif
      print*,'algo falla'
      RETURN
      END 
      DOUBLE PRECISION FUNCTION sex3j(LA,LB,LC,NLF,NBT,NGA)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IF((NLF+NBT+NGA).NE.0.OR.IABS(LA-LB).GT.LC.OR.(LA+LB).LT.LC.
     .OR.IABS(NLF).GT.LA.OR.IABS(NBT).GT.LB.OR.IABS(NGA).GT.LC)THEN
      sex3j=0.0D00
      ELSE
      A1=ea((LA+NLF)/2)+ea((LB+NBT)/2)+ea((LB-NBT)/2)+ea((LC+NGA)/2)+
     +ea((LC-NGA)/2)+ea((LA-NLF)/2)
      D1=eelta(LA,LB,LC)/2.d0
      CNO=((-1)**(iabs(LA-LB-NGA)/2))
C
C    THIS IS DONE TO PREVENT OVERFLOW
C
      SUM=0.0D00
      N1=MAX0(LB-NLF-LC,LA+NBT-LC,0)/2
      N2=MIN0(LA+LB-LC,LA-NLF,LB+NBT)/2
      DO 100 N=N2,N1,-1
      B1=ea(N)+ea((LC-LB+NLF)/2+N)+ea((LC-LA-NBT)/2+N)+
     + ea((LA+LB-LC)/2-N)+ea((LA-NLF)/2-N)+ea((LB+NBT)/2-N)
      SUM=SUM+((-1)**n)*dexp(a1/2-b1)
100   CONTINUE
      sex3j=CNO*SUM*dexp(d1)
      ENDIF
      RETURN
      END
C
*************************** D E L T A ********************************
*--------------------------------------------------------------------*
*         AUXILIARY FUNCTION IN sex3j                                *
*--------------------------------------------------------------------*
C
      DOUBLE PRECISION FUNCTION eelta(LA,LB,LC)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IF((LA+LB+LC).LT.0)THEN
      eelta=0.0D00
      RETURN
      ELSE
      eelta=ea((LA+LB-LC)/2)+ea((LB+LC-LA)/2)+ea((LC+LA-LB)/2)-
     -ea((LA+LB+LC)/2+1)
      RETURN
      ENDIF
      END
C
C
C
************************** F A *****************************
C
      SUBROUTINE CONSea
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/eaCT/Aea(350)
      common/dbfc/dbf(350)
*------------------------------------------*
*     CALCULATES eaCTORIAL(N)              *
*------------------------------------------*
      FB=0.0D00
      DO 10 I=1,350
      FB=FB+dlog(Dfloat(I))
      Aea(I)=FB
 10   CONTINUE
      fd=1.d0
      do 20 j=1,250
      fd=1.d0
      do 30 k=1,j,2
      fd=fd*dfloat(k)
 30   continue
      dbf(j)=fd
 20   continue
      RETURN
      END
      double precision function dbeac(n)
      implicit double precision(a-h,o-z)
      common/dbfc/dbf(350)
      dbeac=dbf(n)
      return
      end
c
c
c
      DOUBLE PRECISION FUNCTION ea(N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/eaCT/Aea(350)
C     IF(N.EQ.-1)THEN
C     ea=1.E100
C     ELSE
      IF(N.EQ.0)THEN
      ea=0.0D0
      ELSE
      ea=Aea(N)
      ENDIF
C     ENDIF
      RETURN
      END

C 
C     CLEBS-GORDAN  <L1,L,M1,M2,L,M>=CG(2*L1,2*L2,2*M1,2*M2,2*J,2*M)   
C                                                                        
      DOUBLE PRECISION FUNCTION CG(J1,m1,j2,m2,J,M)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      CNO=DSQRT(DBLE(J+1))*DBLE((-1)**((J1-J2+M)/2))
      CG=CNO*Sex3J(J1,J2,J,M1,M2,-M)
      RETURN
      END 
C
C     3 J   S I M B O L S, L A , L B , L C = 2 * ( J 1 , J 2 , J 3 )
C
C                    N L F , N B T , N G A = 2 * ( M 1 , M 2 , M 3 )
C
      DOUBLE PRECISION FUNCTION SIM3J(LA,LB,LC,NLF,NBT,NGA) 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IF((NLF+NBT+NGA).NE.0.OR.IABS(LA-LB).GT.LC.OR.(LA+LB).LT.LC.
     .OR.IABS(NLF).GT.LA.OR.IABS(NBT).GT.LB.OR.IABS(NGA).GT.LC)THEN
      SIM3J=0.0D00
      ELSE
      A1=FA((LA+NLF)/2)
      A1=DSQRT(A1)
      A2=FA((LB+NBT)/2)
      A2=DSQRT(A2)
      A3=FA((LB-NBT)/2)
      A3=DSQRT(A3)
      A4=FA((LC+NGA)/2)
      A4=DSQRT(A4)
      A5=FA((LC-NGA)/2)
      A5=DSQRT(A5)
      A6=FA((LA-NLF)/2)
      A6=DSQRT(A6)
      D1=DELTA(LA,LB,LC)
      isg=1
      if(mod((la-lb-nga)/2,2).ne.0)isg=-1
      CNO=isg*DSQRT(D1)
      SUM=0.0D00
      N1=MAX0(LB-NLF-LC,LA+NBT-LC,0)/2
      N2=MIN0(LA+LB-LC,LA-NLF,LB+NBT)/2 
      DO 100 N=N2,N1,-1
      B1=FA(N)/A1
      B2=FA((LC-LB+NLF)/2+N)/A2
      B3=FA((LC-LA-NBT)/2+N)/A3
      B4=FA((LA+LB-LC)/2-N)/A4
      B5=FA((LA-NLF)/2-N)/A5
      B6=FA((LB+NBT)/2-N)/A6
      CNTER=DBLE((-1)**N)*B1*B2*B3*B4*B5*B6 
      IF(CNTER.NE.0.0D00)SUM=SUM+1.D0/CNTER
100   CONTINUE
      SIM3J=CNO*SUM 
      ENDIF
      RETURN
      END 
C
C     F A C T O R I A L ( N )  
C

      SUBROUTINE CONSFA
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/FACT/AFA(350)
      FB=1.0D00
      DO 10 I=1,250 
      FB=FB*DBLE(I) 
      AFA(I)=FB
 10   CONTINUE
      RETURN
      END 
C
      DOUBLE PRECISION FUNCTION FA(N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/FACT/AFA(350)
      IF(N.EQ.0)THEN
      FA=1.0D0
      ELSE
      FA=AFA(N)
      ENDIF
      RETURN
      END 
      DOUBLE PRECISION FUNCTION SIM6J(J1,J2,J3,J4,J5,J6)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      SUM=0.D0
      DO 100 N=MAX0(J1+J2+J3,J1+J5+J6,J4+J2+J6,J3+J4+J5),
     ,MIN0(J1+J2+J4+J5,J2+J3+J5+J6,J3+J1+J6+J4),2
      C1=dble((-1)**(N/2))*FA((N+2)/2)/FA((N-J1-J2-J3)/2)
      C2=FA((N-J1-J5-J6)/2)*FA((N-J4-J2-J6)/2)
      C3=FA((N-J4-J5-J3)/2)*FA((J1+J2+J4+J5-N)/2)
      C4=FA((J2+J3+J5+J6-N)/2)*FA((J3+J1+J6+J4-N)/2)
      SUM=SUM+C1/C2/C3/C4
  100 CONTINUE
      IF(SUM.NE.0.D0)THEN
      SIM6J=SUM*DSQRT(delta(j1,j2,j3)*delta(j1,j5,j6)*delta(j4,j2,j6)*
     *delta(j4,j5,j3))
      else
      sim6j=0.0d0
      endif
      return
      end
C 
C     R a c a h   W - c o e f f i c i e n t
C
C         a s   a   f u n c t i o n   
C
C        o f   6 - j   S y m b o l s
C
      function w(j1,j2,l2,l1,j3,l3)
      implicit double precision(a-h,o-z)
      fase=dble((-1)**((j1+j2+l1+l2)/2))
      w=fase*sim6j(j1,j2,j3,l1,l2,l3)
      return
      end
      DOUBLE PRECISION FUNCTION DELTA(LA,LB,LC)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IF((LA+LB+LC).LT.0)THEN 
      DELTA=0.0D00
      RETURN
      ELSE
      DELTA=FA((LA+LB-LC)/2)*FA((LB+LC-LA)/2)*FA((LC+LA-LB)/2)/
     1FA((LA+LB+LC)/2+1)
      RETURN
      ENDIF
      END
