      FUNCTION RF(X,Y,Z)
C
C     * NUMERICAL RECIPES ROUTINE USED TO EVALUATE INCOMPLETE
C     * ELLIPTIC INTEGRALS. CALLED BY GWLOOKU THROUGH FUNCTIONS
C     * EF,KF
C     * "PAUSE" STATEMENT IS REPLACED WITH WITH A CALL TO THE 
C     * "XIT" SUBROUTINE.
C     (C) COPR. 1986-92 NUMERICAL RECIPES SOFTWARE #=!E=#,)]UBCJ
C
C     * MAR 24/99 J SCINOCCA
C
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)
      REAL RF,X,Y,Z,ERRTOL,TINY,BIG,THIRD,C1,C2,C3,C4

      PARAMETER (ERRTOL=0.08E0,TINY=1.5E-76,BIG=3.E74,THIRD=1.E0/3.E0,
     1           C1=1.E0/24.E0,C2=0.1E0,C3=3.E0/44.E0,C4=1.E0/14.E0)

      REAL ALAMB,AVE,DELX,DELY,DELZ,E2,E3,SQRTX,SQRTY,SQRTZ,
     1     XT,YT,ZT
C
      IF(MIN(X,Y,Z).LT.0.E0.OR.MIN(X+Y,X+Z,Y+Z).LT.TINY
     1                   .OR.MAX(X,Y,Z).GT.BIG) THEN
        WRITE(6,'(A)') ' INVALID ARGUMENTS IN RF'
        CALL                                       XIT ('RF',-1)
      ENDIF

      XT=X
      YT=Y
      ZT=Z

   1  CONTINUE
        SQRTX=SQRT(XT)
        SQRTY=SQRT(YT)
        SQRTZ=SQRT(ZT)
        ALAMB=SQRTX*(SQRTY+SQRTZ)+SQRTY*SQRTZ
        XT=0.25E0*(XT+ALAMB)
        YT=0.25E0*(YT+ALAMB)
        ZT=0.25E0*(ZT+ALAMB)
        AVE=THIRD*(XT+YT+ZT)
        DELX=(AVE-XT)/AVE
        DELY=(AVE-YT)/AVE
        DELZ=(AVE-ZT)/AVE
      IF(MAX(ABS(DELX),ABS(DELY),ABS(DELZ)).GT.ERRTOL)GOTO 1

      E2=DELX*DELY-DELZ**2
      E3=DELX*DELY*DELZ
      RF=(1.E0+(C1*E2-C2-C3*E3)*E2+C4*E3)/SQRT(AVE)

      RETURN
      END
