
      subroutine ffactor(idip,icharg,ireac,qnu2,ef1,ef2,eg1,fcoup)
      implicit real*8 (a-h,o-z)
      real*8  eff1(2),eff2(2),
     ,         ff1(2),ff2(2),f1(0:8),f2(0:8),g1(0:8),gg1(2)

 

      akp=1.79285d0

      akn=-1.91304d0





      alph=1.d0/137.04

      pi=3.1415926    





      

      sin2w=0.2330

      Gf=1.1664e-11



      sqr3=sqrt(3.0)

      amass=939.d0

      hc=197.3d0



      am2v=840.d0**2

      am2a=1030.d0**2





      f2s=-0.21

      g1s=-0.19

c      f2s=0.

c      g1s=0.





      f1s=0.

      f1(0)=1.+f1s

      f2(0)=f2s+(akp+akn)

      g1(0)=g1s+0.55/3.



      f1(3)=0.5

      f2(3)=0.5*(akp-akn)

      g1(3)=0.63



      f1(8)=0.5*sqr3

      f2(8)=0.5*sqr3*(akp+akn)

      g1(8)=0.275/sqr3



      z0v=-0.5

      z0a=0.5

      z3v=1.-2*sin2w

      z3a=-1.

      z8v=1./sqr3*(1.-2*sin2w)

      z8a=-1./sqr3

c      print*,z0v,z3v,z8v,f2(0),f2(3),f2(8)

c      write(10,*)'proton-nu   proton anu  neutron nu  neutron-anu'



      Q2=-qnu2*hc*hc





      dip=1.+q2/am2v

      dip=1./dip/dip



      dips=1.+Q2/am2a

      dips=1./dips/dips



      tau=Q2/4./amass/amass







c------------------------------------

c

c     factores de forma EM del neutron y proton

c





      if (icharg.eq.1)then

        anom=akp

        ge0=1.0

      else

        anom=akn

        ge0=0.0

      endif









      Gm0=ge0+anom









      if(idip.eq.0)then

        dip=1.d0

        dips=1.d0

      endif



      gm=gm0*dip

      ge=ge0*dip



      if(icharg.eq.2)then

      ge=-anom*tau/(1.+5.6*tau)*dip

      endif







      eff2(icharg)=(gm-ge)/(1.+tau)/anom

      eff1(icharg)=(ge+tau*gm)/(1.+tau)







      gmt=eff1(icharg)+anom*eff2(icharg)

      get=eff1(icharg)-tau*anom*eff2(icharg)

c      print*,'qnu2=',qnu2,ff1(1),ff2(1),gm/ge,gmt/get

c      if(ipl.eq.2)ff1=-anom*ff2

c      if(ipl.eq.3)ff2=0.0d0

c      if(ipl.eq.4)ff1=0.d0

c      falp=ff1

c      fbet=anom*ff2/2.d0/amassp



c      fff(1)=-anom*ff2/2.d0/amassp

c      fff(2)=ff1+anom*ff2

c      print*,fff(1),fff(2),'momentos'















      f20=f2(0)/(1.+tau)*dip

      f23=f2(3)/(1.+tau)*dip

      f28=f2(8)/(1.+tau)*dip



      f10=f1(0)*dip+tau*f20

      f13=f1(3)*dip+tau*f23 

      f18=f1(8)*dip+tau*f28

      

      g10=g1(0)*dips

      g13=g1(3)*dips

      g18=g1(8)*dips





      fac=gf*gf/2.

c     /3.1415926*(197.3**2)*2.*amass*1e16



c     cross section in fm^2/MeV x 1e-16=cm2/MeV x 1.e-42)





      ff1(1)=z0v*f10+z3v*f13+z8v*f18

      ff1(2)=z0v*f10-z3v*f13+z8v*f18



      ff2(1)=z0v*f20+z3v*f23+z8v*f28

      ff2(2)=z0v*f20-z3v*f23+z8v*f28



      gg1(1)=z0a*g10+z3a*g13+z8a*g18

      gg1(2)=z0a*g10-z3a*g13+z8a*g18



c      print*,'ff1(1),ff1(2),gg1(1),gg1(2)',ff1(1),ff1(2),gg1(1),gg1(2)

c      print*,'ff2 ',ff2(1),ff2(2)

      i=icharge



      write(13,'(i4,5e13.5)')i,q2,ff1(i),ff2(i),gg1(i)

c     i=1, proton

c     i=2, neutron



      if(ireac.ne.0)then



      ef1=ff1(icharg)

      ef2=ff2(icharg)

      eg1=gg1(icharg)



      fcoup=fac 



      else



      eg1=0.d0

      ef1=eff1(icharg)

      ef2=eff2(icharg)*anom

      fcoup=2.*(2.*pi)**2*alph*alph/qnu2/qnu2



      endif      

     

      write(57,'(4f10.5,f10.6,2f10.5)')qnu2,ef1,ef2,eg1,fcoup,
     ,                                 gm/ge,gmt/get



      end











      SUBROUTINE harmsq(spci,dpci)

      IMPLICIT real*8 (A-H,O-Z)

      include 'param.inc'

      common/harm/ysh(0:npw+1,-npw:npw)

*--------------------------------------------------*

*     CALCULATES spherical harmonics               *

*--------------------------------------------------*

c      print900,spci,dpci

c      write(13,900)spci,dpci

 900  format('sin in.=',f10.5,' cos in.=',f10.5,' sin fin.=',f10.5,

     ,' cos fin.=',f10.5)

      do 137 l=1,npw

      do 137 m=-npw,npw

      ysh(l,m)=0.d0

 137  continue

      PI=dacos(-1.d0)

      FOURPI=4.d0*PI

      spc=spci

      dpc=dpci

      if(dabs(spc).lt.1.d-7)then

        do 111 l=0,npw

        if=1

        if(dpc.lt.0.d0)if=(-1)**l

        ysh(l,0)=if*dsqrt((02.d0*l+1.d0)/fourpi)

 111  continue

      else



        DO 100 M=-npw,npw

*-------------------------------------------------------------------*

*     RECURENCE FORMULA FOR THE SPHER. HARM. ,KEEPTS IN YSH(L,M)    *

*     SUM IS THE TOTAL AMPLITUDE                                    *

*-------------------------------------------------------------------*

        A1=1.0D00

        A3=3.0D00

        IM=IABS(M)

        DM=Dfloat(M)

        DIM=Dfloat(IM)

C     DO 101 L=0,IM-1

C     YSH(L,M)=0.0D00

C101  CONTINUE

      IF(IM.EQ.0)THEN

      AUX1=1.0D0

      ELSE

      AUX1=SPC**IM

      if=1

      if(aux1.lt.0.d0)if=-1

      aux1=aux1*if**(im)

      ENDIF

      is=(-1)**im

      if(m.lt.0)is=1

c      print*,m,im,is,aux1,spc,2.d0**im,ea(im)

      YSH(IM,M)=is*DSQRT((2*DIM+A1)*dexp(eA(2*IM)))/DSQRT(FOURPI)/

     1((2.d0**IM)*dexp(eA(IM)))*AUX1

      YSH(IM+1,M)=DSQRT((2*DIM+A1)*(2*DIM+A3)/(DIM+DM+A1)/

     1(DIM-DM+A1))*DPC*YSH(IM,M)

c      print*,m,im,is,aux1,spc,2.d0**im,ea(im)

      DO 102 L=IM+1,npw

      DL=dfloat(L)

      CT=DSQRT((2*DL+A1)*(2*DL+A3)/(DL+A1+DM)/(DL+A1-DM))

      CI=DSQRT((DL+DM)*(DL-DM)/(2*DL+A1)/(2*DL-A1))

      YSH(L+1,M)=CT*(DPC*YSH(L,M)-CI*YSH(L-1,M))

 102  CONTINUE

 100  CONTINUE

      ENDIF

*-------------------------------------------------------------*

*     THIS PRINTS A TEST, Y10=DSQRT(3/4PI)*DCOS(THET)...       *

*-------------------------------------------------------------*

c      WRITE(*,99)YSHi(2,0)/DSQRT(5.D0/FOURPI),

c     ,YSHi(10,0)/DSQRT(21.D0/FOURPI),

c     ,YSHi(2,0),DSQRT(1.25D0/FOURPI)*(3.*DPC*DPC-1.0),

c     ,YSHi(2,-1),YSHi(2,1)

 99   FORMAT(2E12.5)

      return

      end

C

c

c

*----------------------------------------------------------------------*

*----------------------------------------------------------------------*

*  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





      real*8 function djl(l,xpm)

      implicit real*8 (a-h,o-z)





      if(l.eq.0) then

      djl=dsin(xpm)/xpm

      return

      endif



      if(l.eq.1) then

      djl=dsin(xpm)/xpm/xpm-dcos(xpm)/xpm

      return

      endif



      if(l.eq.2) then

      djl=(3.d0/xpm/xpm/xpm-1.d0/xpm)*dsin(xpm)-3.d0/xpm/xpm*dcos(xpm)

      return

      endif





      if(l.eq.3) then

      djl=dcos(xpm)-15.d0*dcos(xpm)/xpm/xpm+15.d0*dsin(xpm)/xpm**3-

     -    6.d0*dsin(xpm)/xpm

      djl=djl/xpm

      return

      endif



      if(l.eq.4) then

      djl=-105d0*dcos(xpm)/xpm**3+10.d0*dcos(xpm)/xpm+dsin(xpm)+

     +     105.d0*dsin(xpm)/xpm**4-45.d0*dsin(xpm)/xpm/xpm

      djl=djl/xpm

      return

      endif





      if(l.ge.5)then

      print*,' preparar version con hl(pm), esta ya no es eficiente'

      stop

      return

      endif



      end 







      SUBROUTINE HL(X)

      IMPLICIT  double precision(A-H,O-Z)

      include 'param.inc'

      COMMON/UNO/Hp(0:100)

      dimension ahl(0:120)

      if(x.lt.0.025d0)then

      hp(0)=1.d0

      do 10 l=1,npw+3

c      print*,dbeac(l+l+1)

      hp(l)=x**l/(dbeac(l+l+1))

 10   continue

      return

      endif

      IMAX=50+idint(x)



c      imax=60

      aHL(IMAX)=1.D-280 

      aHL(IMAX-1)=0.0D0

      DO 100 I=IMAX-1,1,-1

      aHL(I-1)=(2.D0*Dfloat(I)+1.D0)/X*aHL(I)-aHL(I+1) 

 100  CONTINUE

      REF=DSIN(X)/X 

      ESCALA=REF/aHL(0)

      DO 200 I=0,npw+3

      Hp(I)=ESCALA*aHL(I)

 200  CONTINUE

      RETURN

      END 

      SUBROUTINE azh(L,X,ZH)

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      COMPLEX*16 ZI,ZH,ZEROZ,coef

      ZI=dCMPLX(0.D0,1.D0)

      ZEROZ=dCMPLX(0.D0,0.D0)

      ZH=ZEROZ

      DO 100 N=1,l+1

      COEF=FA(L+N-1)/FA(N-1)/FA(1+L-N)*(2.D0**(1-N))

      ZH=ZH+COEF/X**N*(ZI**(2+N-L))

 100  CONTINUE

      ZH=ZH*CdEXP(ZI*X)

      RETURN

      END 

C

      SUBROUTINE SPLINE(X,Y,N,YP1,YPN,Y2)

c      PARAMETER (NMAX=npunt1)

      implicit  real*8 (a-h,o-z)

      include 'param.inc'

      DIMENSION X(npunt2),Y(npunt2),Y2(npunt2),U(npunt2)

      IF (YP1.GT..99E30) THEN

        Y2(1)=0.d0

        U(1)=0.d0

      ELSE

        Y2(1)=-0.5d0

        U(1)=(3.d0/(X(2)-X(1)))*((Y(2)-Y(1))/(X(2)-X(1))-YP1)

      ENDIF

      DO 11 I=2,N-1

        SIG=(X(I)-X(I-1))/(X(I+1)-X(I-1))

        P=SIG*Y2(I-1)+2.d0

        Y2(I)=(SIG-1.d0)/P

        U(I)=(6.d0*((Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1))

     *      /(X(I)-X(I-1)))/(X(I+1)-X(I-1))-SIG*U(I-1))/P

11    CONTINUE

      IF (YPN.GT..99d30) THEN

        QN=0.d0

        UN=0.d0

      ELSE

        QN=0.5d0

        UN=(3.d0/(X(N)-X(N-1)))*(YPN-(Y(N)-Y(N-1))/(X(N)-X(N-1)))

      ENDIF

      Y2(N)=(UN-QN*U(N-1))/(QN*Y2(N-1)+1.d0)

      DO 12 K=N-1,1,-1

        Y2(K)=Y2(K)*Y2(K+1)+U(K)

12    CONTINUE

      RETURN

      END

 

      SUBROUTINE SPLINT(XA,YA,Y2A,N,X,Y)

      implicit real*8 (a-h,o-z)

      include 'param.inc'

      DIMENSION XA(npunt2),YA(npunt2),Y2A(npunt2)

      KLO=1

      KHI=N

1     IF (KHI-KLO.GT.1) THEN

        K=(KHI+KLO)/2

        IF(XA(K).GT.X)THEN

          KHI=K

        ELSE

          KLO=K

        ENDIF

      GOTO 1

      ENDIF

      H=XA(KHI)-XA(KLO)

      IF (H.EQ.0.d0) PAUSE 'Bad XA input.'

      A=(XA(KHI)-X)/H

      B=(X-XA(KLO))/H

      Y=A*YA(KLO)+B*YA(KHI)+

     *      ((A**3-A)*Y2A(KLO)+(B**3-B)*Y2A(KHI))*(H**2)/6.d0

      RETURN

      END

      subroutine zeul(zs,zt,j,zw,jm)     

      implicit real*8 (a-h,o-z)

      include 'param.inc'

      dimension rzw(0:npw),dizw(0:npw),zw(0:j)

      complex*16 zs,zt,zw

      do 100 k=0,j

c        print*,k,j

        rzw(k)=dble(zw(k))

        dizw(k)=dimag(zw(k))

 100  continue

      rterm=dble(zt)

      diter=dimag(zt)

      call eulsum(rsum,rterm,j,rzw)

      call eulsum(disum,diterm,j,dizw)

      zs=dcmplx(rsum,disum)                  

      do 101 k=0,j

        print*,j,k,npw,'    estamos aqui'

        zw(k)=dcmplx(rzw(k),dizw(k))

        print*,j,k,npw,zw(k),'    estamos aqui'

 101  continue

      return

      end

c

c

c

      SUBROUTINE EULSUM(SUM,TERM,JTERM,WKSP)

      implicit real*8 (a-h,o-z)

      include 'param.inc'

      DIMENSION WKSP(0:jterm+1)

c      print*,'jterm=',jterm

      IF(JTERM.EQ.0)THEN

        NTERM=1

        WKSP(0)=TERM

        SUM=0.5d0*TERM

      ELSE

        TMP=WKSP(0)

        WKSP(0)=TERM

        DO 11 J=0,NTERM

c

c          print*,jterm,j+1

c

          DUM=WKSP(J+1)

          WKSP(J+1)=0.5d0*(WKSP(J)+TMP)

          TMP=DUM

11      CONTINUE

        IF(dABS(WKSP(NTERM+1)).LE.dABS(WKSP(NTERM)))THEN

c        print*,nterm+1

          SUM=SUM+0.5d0*WKSP(NTERM+1)

          NTERM=NTERM+1

        ELSE

          SUM=SUM+WKSP(NTERM+1)

        ENDIF

      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

