jelp Subroutine

subroutine jelp ( u, hk, esn, ecn, edn, eph )

  !*****************************************************************************80
  !
  !! JELP computes Jacobian elliptic functions SN(u), CN(u), DN(u).
  !
  !  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:
  !
  !    08 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 ) U, the argument.
  !
  !    Input, real ( kind = 8 ) HK, the modulus, between 0 and 1.
  !
  !    Output, real ( kind = 8 ) ESN, ECN, EDN, EPH, the values of
  !    sn(u), cn(u), dn(u), and phi (in degrees).
  !
  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 ) d
  real ( kind = 8 ) dn
  real ( kind = 8 ) ecn
  real ( kind = 8 ) edn
  real ( kind = 8 ) eph
  real ( kind = 8 ) esn
  real ( kind = 8 ) hk
  integer ( kind = 4 ) j
  integer ( kind = 4 ) n 
  real ( kind = 8 ) pi
  real ( kind = 8 ) r(40)
  real ( kind = 8 ) sa
  real ( kind = 8 ) t
  real ( kind = 8 ) u

  pi = 3.14159265358979D+00
  a0 = 1.0D+00
  b0 = sqrt ( 1.0D+00 - hk * hk )

  do n = 1, 40

     a = ( a0 + b0 ) / 2.0D+00
     b = sqrt ( a0 * b0 )
     c = ( a0 - b0 ) / 2.0D+00
     r(n) = c / a

     if ( c < 1.0D-07 ) then
        exit
     end if

     a0 = a
     b0 = b

  end do

  dn = 2.0D+00 ** n * a * u

  do j = n, 1, -1
     t = r(j) * sin ( dn )
     sa = atan ( t / sqrt ( abs ( 1.0D+00 - t * t )))
     d = 0.5D+00 * ( dn + sa )
     dn = d
  end do

  eph = d * 180.0D+00 / pi
  esn = sin ( d )
  ecn = cos ( d )
  edn = sqrt ( 1.0D+00 - hk * hk * esn * esn )

  return
end subroutine jelp