refine Subroutine

subroutine refine ( kd, m, q, a, iflag )

  !*****************************************************************************80
  !
  !! REFINE refines an estimate of the characteristic value of Mathieu functions.
  !
  !  Discussion:
  !
  !    This procedure calculates the accurate characteristic value
  !    by the secant method.
  !
  !  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:
  !
  !    20 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, integer ( kind = 4 ) KD, the case code:
  !    1, for cem(x,q)  ( m = 0,2,4,...)
  !    2, for cem(x,q)  ( m = 1,3,5,...)
  !    3, for sem(x,q)  ( m = 1,3,5,...)
  !    4, for sem(x,q)  ( m = 2,4,6,...)
  !
  !    Input, integer ( kind = 4 ) M, the order of the Mathieu functions.
  !
  !    Input, real ( kind = 8 ) Q, the parameter of the Mathieu functions.
  !
  !    Input/output, real ( kind = 8 ) A, the characteristic value, which
  !    should have been refined on output.
  !
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) ca
  real ( kind = 8 ) delta
  real ( kind = 8 ) eps
  real ( kind = 8 ) f
  real ( kind = 8 ) f0
  real ( kind = 8 ) f1
  integer ( kind = 4 ) it
  integer ( kind = 4 ) iflag
  integer ( kind = 4 ) kd
  integer ( kind = 4 ) m
  integer ( kind = 4 ) mj
  real ( kind = 8 ) q
  real ( kind = 8 ) x
  real ( kind = 8 ) x0
  real ( kind = 8 ) x1

  eps = 1.0D-14
  mj = 10 + m
  ca = a
  delta = 0.0D+00
  x0 = a
  call cvf ( kd, m, q, x0, mj, f0 )
  x1 = 1.002D+00 * a
  call cvf ( kd, m, q, x1, mj, f1 )

  do

     do it = 1, 100
        mj = mj + 1
        x = x1 - ( x1 - x0 ) / ( 1.0D+00 - f0 / f1 )
        call cvf ( kd, m, q, x, mj, f )
        if ( abs ( 1.0D+00 - x1 / x ) < eps .or. f == 0.0D+00 ) then
           exit
        end if
        x0 = x1
        f0 = f1
        x1 = x
        f1 = f
     end do

     a = x

     if ( 0.05D+00 < delta ) then
        a = ca
        if ( iflag < 0 ) then
           iflag = -10
        end if
        return
     end if

     if ( abs ( ( a - ca ) / ca )  <= 0.05D+00 ) then
        exit
     end if

     x0 = ca
     delta = delta + 0.005D+00
     call cvf ( kd, m, q, x0, mj, f0 )
     x1 = ( 1.0D+00 + delta ) * ca
     call cvf ( kd, m, q, x1, mj, f1 )

  end do

  return
end subroutine refine