elit Subroutine

subroutine elit ( hk, phi, fe, ee )

  !*****************************************************************************80
  !
  !! ELIT: complete and incomplete elliptic integrals F(k,phi) and E(k,phi).
  !
  !  Licensing:
  !
  !    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However, 
  !    they give permission to incorporate this routine into a user program 
  !    provided that the copyright is acknowledged.
  !
  !  Modified:
  !
  !    12 July 2012
  !
  !  Author:
  !
  !    Shanjie Zhang, Jianming Jin
  !
  !  Reference:
  !
  !    Shanjie Zhang, Jianming Jin,
  !    Computation of Special Functions,
  !    Wiley, 1996,
  !    ISBN: 0-471-11963-6,
  !    LC: QA351.C45.
  !
  !  Parameters:
  !
  !    Input, real ( kind = 8 ) HK, the modulus, between 0 and 1.
  !
  !    Input, real ( kind = 8 ) PHI, the argument in degrees.
  !
  !    Output, real ( kind = 8 ) FE, EE, the values of F(k,phi) and E(k,phi).
  !
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) a0
  real ( kind = 8 ) b
  real ( kind = 8 ) b0
  real ( kind = 8 ) c
  real ( kind = 8 ) ce
  real ( kind = 8 ) ck
  real ( kind = 8 ) d
  real ( kind = 8 ) d0
  real ( kind = 8 ) ee
  real ( kind = 8 ) fac
  real ( kind = 8 ) fe
  real ( kind = 8 ) g
  real ( kind = 8 ) hk
  integer ( kind = 4 ) n
  real ( kind = 8 ) phi
  real ( kind = 8 ) pi
  real ( kind = 8 ) r

  g = 0.0D+00
  pi = 3.14159265358979D+00
  a0 = 1.0D+00
  b0 = sqrt ( 1.0D+00 - hk * hk )
  d0 = ( pi / 180.0D+00 ) * phi
  r = hk * hk

  if ( hk == 1.0D+00 .and. phi == 90.0D+00 ) then

     fe = 1.0D+300
     ee = 1.0D+00

  else if ( hk == 1.0D+00 ) then

     fe = log ( ( 1.0D+00 + sin ( d0 ) ) / cos ( d0 ) )
     ee = sin ( d0 )

  else

     fac = 1.0D+00
     do n = 1, 40
        a = ( a0 + b0 ) /2.0D+00
        b = sqrt ( a0 * b0 )
        c = ( a0 - b0 ) / 2.0D+00
        fac = 2.0D+00 * fac
        r = r + fac * c * c
        if ( phi /= 90.0D+00 ) then
           d = d0 + atan ( ( b0 / a0 ) * tan ( d0 ) )
           g = g + c * sin( d )
           d0 = d + pi * int ( d / pi + 0.5D+00 )
        end if
        a0 = a
        b0 = b
        if ( c < 1.0D-07 ) then
           exit
        end if
     end do

     ck = pi / ( 2.0D+00 * a )
     ce = pi * ( 2.0D+00 - r ) / ( 4.0D+00 * a )
     if ( phi == 90.0D+00 ) then
        fe = ck
        ee = ce
     else
        fe = d / ( fac * a )
        ee = fe * ce / ck + g
     end if

  end if

  return
end subroutine elit