ADAS Subroutine r8fdip0

      function r8fdip0(e1,l1,e2,l2,eps)
C-----------------------------------------------------------------------
C
C  ****************** fortran77 program: r8fdip0.for *******************
C
C  purpose: calculates the function i0(k1,l1,k2,l2,1) defined in
C           Phil. Trans. Roy. Soc. a266,255,1970, where
C           e1=k1*k1, e2=k2*k2, and the relative accuracy is
C           approximately eps.
C
C           It is suitable for use in equations (13) etc. of
C           J.Phys.B. 7,l364,1974
C           (original by A. Burgess, DAMTP, University of Cambridge)
C
C
C  subroutine:
C
C  input : (r*8)  e1       = initial electron energy (Ryd)
C  input : (r*8)  l1       = inital orbital angular momentum
C  input : (r*8)  e2       = final electron energy (Ryd)
C  input : (r*8)  l2       = final orbital angular momentum
C  input : (r*8)  eps      = accuracy setting
C
C  output: (r*8)  r8fdip0  = dipole matrix element for neutral atom
C
C
C  Routines:
C          routine    source    brief description
C          -------------------------------------------------------------
C          f21        adas      special quadrature for Burgess codes
C          i4unit     adas      fetch unit number for output of messages
C
C
C  Author:  H. P. Summers, University of Strathclyde
C           ja7.08
C           tel. 0141-548-4196
C
C  Date:   24/02/03
C
C  Update: HP Summers    24/05/04
C          restructure and addded standard warning
C  Update: AD Whiteford  16/03/05
C          renamed to r8fdip
C          Now calls r8f21 instead of just f21, this routine was renamed
C  Update: AD Whiteford  17/05/07
C          Updated comments as part of subroutine documentation
C          procedure.
C
C-----------------------------------------------------------------------
      INTEGER             L1,          L2
      REAL*8              E1,          E2,          EPS