quadpack.f90 Source File


Contents

Source Code


Source Code

! Downloaded from http://people.sc.fsu.edu/~jburkardt/f_src/quadpack/quadpack.f90
! Written by Robert Piessens, Elise deDoncker-Kapenga, Christian Ueberhuber, David Kahaner.
! Reference:
! Robert Piessens, Elise deDoncker-Kapenga, Christian Ueberhuber, David Kahaner,
! QUADPACK: A Subroutine Package for Automatic Integration,
! Springer, 1983, ISBN: 3540125531.

module quadpack
use warning_helpers, only: is_zero, is_not_zero, exactly_equal, not_exactly_equal
implicit none

private :: is_zero, is_not_zero, exactly_equal, not_exactly_equal

interface
  real function scalar_func(x_)
    real, intent(in) :: x_
  end function scalar_func
end interface

contains
subroutine aaaa

!*****************************************************************************80
!
!! AAAA is a dummy subroutine with QUADPACK documentation in its comments.
!
! 1. introduction
!
!    quadpack is a fortran subroutine package for the numerical
!    computation of definite one-dimensional integrals. it originated
!    from a joint project of r. piessens and e. de doncker (appl.
!    math. and progr. div.- k.u.leuven, belgium), c. ueberhuber (inst.
!    fuer math.- techn.u.wien, austria), and d. kahaner (nation. bur.
!    of standards- washington d.c., u.s.a.).
!
! 2. survey
!
!    - qags : is an integrator based on globally adaptive interval
!             subdivision in connection with extrapolation (de doncker,
!             1978) by the epsilon algorithm (wynn, 1956).
!
!    - qagp : serves the same purposes as qags, but also allows
!             for eventual user-supplied information, i.e. the
!             abscissae of internal singularities, discontinuities
!             and other difficulties of the integrand function.
!             the algorithm is a modification of that in qags.
!
!    - qagi : handles integration over infinite intervals. the
!             infinite range is mapped onto a finite interval and
!             then the same strategy as in qags is applied.
!
!    - qawo : is a routine for the integration of cos(omega*x)*f(x)
!             or sin(omega*x)*f(x) over a finite interval (a,b).
!             omega is is specified by the user
!             the rule evaluation component is based on the
!             modified clenshaw-curtis technique.
!             an adaptive subdivision scheme is used connected with
!             an extrapolation procedure, which is a modification
!             of that in qags and provides the possibility to deal
!             even with singularities in f.
!
!    - qawf : calculates the fourier cosine or fourier sine
!             transform of f(x), for user-supplied interval (a,
!             infinity), omega, and f. the procedure of qawo is
!             used on successive finite intervals, and convergence
!             acceleration by means of the epsilon algorithm (wynn,
!             1956) is applied to the series of the integral
!             contributions.
!
!    - qaws : integrates w(x)*f(x) over (a,b) with a < b finite,
!             and   w(x) = ((x-a)**alfa)*((b-x)**beta)*v(x)
!             where v(x) = 1 or log(x-a) or log(b-x)
!                            or log(x-a)*log(b-x)
!             and   -1 < alfa, -1 < beta.
!             the user specifies a, b, alfa, beta and the type of
!             the function v.
!             a globally adaptive subdivision strategy is applied,
!             with modified clenshaw-curtis integration on the
!             subintervals which contain a or b.
!
!    - qawc : computes the cauchy principal value of f(x)/(x-c)
!             over a finite interval (a,b) and for
!             user-determined c.
!             the strategy is globally adaptive, and modified
!             clenshaw-curtis integration is used on the subranges
!             which contain the point x = c.
!
!  each of the routines above also has a "more detailed" version
!    with a name ending in e, as qage.  these provide more
!    information and control than the easier versions.
!
!
!   the preceeding routines are all automatic.  that is, the user
!      inputs his problem and an error tolerance.  the routine
!      attempts to perform the integration to within the requested
!      absolute or relative error.
!   there are, in addition, a number of non-automatic integrators.
!      these are most useful when the problem is such that the
!      user knows that a fixed rule will provide the accuracy
!      required.  typically they return an error estimate but make
!      no attempt to satisfy any particular input error request.
!
!      qk15
!      qk21
!      qk31
!      qk41
!      qk51
!      qk61
!           estimate the integral on [a,b] using 15, 21,..., 61
!           point rule and return an error estimate.
!      qk15i 15 point rule for (semi)infinite interval.
!      qk15w 15 point rule for special singular weight functions.
!      qc25c 25 point rule for cauchy principal values
!      qc25o 25 point rule for sin/cos integrand.
!      qmomo integrates k-th degree chebychev polynomial times
!            function with various explicit singularities.
!
! 3. guidelines for the use of quadpack
!
!    here it is not our purpose to investigate the question when
!    automatic quadrature should be used. we shall rather attempt
!    to help the user who already made the decision to use quadpack,
!    with selecting an appropriate routine or a combination of
!    several routines for handling his problem.
!
!    for both quadrature over finite and over infinite intervals,
!    one of the first questions to be answered by the user is
!    related to the amount of computer time he wants to spend,
!    versus his -own- time which would be needed, for example, for
!    manual subdivision of the interval or other analytic
!    manipulations.
!
!    (1) the user may not care about computer time, or not be
!        willing to do any analysis of the problem. especially when
!        only one or a few integrals must be calculated, this attitude
!        can be perfectly reasonable. in this case it is clear that
!        either the most sophisticated of the routines for finite
!        intervals, qags, must be used, or its analogue for infinite
!        intervals, qagi. these routines are able to cope with
!        rather difficult, even with improper integrals.
!        this way of proceeding may be expensive. but the integrator
!        is supposed to give you an answer in return, with additional
!        information in the case of a failure, through its error
!        estimate and flag. yet it must be stressed that the programs
!        cannot be totally reliable.
!
!    (2) the user may want to examine the integrand function.
!        if bad local difficulties occur, such as a discontinuity, a
!        singularity, derivative singularity or high peak at one or
!        more points within the interval, the first advice is to
!        split up the interval at these points. the integrand must
!        then be examinated over each of the subintervals separately,
!        so that a suitable integrator can be selected for each of
!        them. if this yields problems involving relative accuracies
!        to be imposed on -finite- subintervals, one can make use of
!        qagp, which must be provided with the positions of the local
!        difficulties. however, if strong singularities are present
!        and a high accuracy is requested, application of qags on the
!        subintervals may yield a better result.
!
!        for quadrature over finite intervals we thus dispose of qags
!        and
!        - qng for well-behaved integrands,
!        - qag for functions with an oscillating behavior of a non
!          specific type,
!        - qawo for functions, eventually singular, containing a
!          factor cos(omega*x) or sin(omega*x) where omega is known,
!        - qaws for integrands with algebraico-logarithmic end point
!          singularities of known type,
!        - qawc for cauchy principal values.
!
!        remark
!
!        on return, the work arrays in the argument lists of the
!        adaptive integrators contain information about the interval
!        subdivision process and hence about the integrand behavior:
!        the end points of the subintervals, the local integral
!        contributions and error estimates, and eventually other
!        characteristics. for this reason, and because of its simple
!        globally adaptive nature, the routine qag in particular is
!        well-suited for integrand examination. difficult spots can
!        be located by investigating the error estimates on the
!        subintervals.
!
!        for infinite intervals we provide only one general-purpose
!        routine, qagi. it is based on the qags algorithm applied
!        after a transformation of the original interval into (0,1).
!        yet it may eventuate that another type of transformation is
!        more appropriate, or one might prefer to break up the
!        original interval and use qagi only on the infinite part
!        and so on. these kinds of actions suggest a combined use of
!        different quadpack integrators. note that, when the only
!        difficulty is an integrand singularity at the finite
!        integration limit, it will in general not be necessary to
!        break up the interval, as qagi deals with several types of
!        singularity at the boundary point of the integration range.
!        it also handles slowly convergent improper integrals, on
!        the condition that the integrand does not oscillate over
!        the entire infinite interval. if it does we would advise
!        to sum succeeding positive and negative contributions to
!        the integral -e.g. integrate between the zeros- with one
!        or more of the finite-range integrators, and apply
!        convergence acceleration eventually by means of quadpack
!        subroutine qelg which implements the epsilon algorithm.
!        such quadrature problems include the fourier transform as
!        a special case. yet for the latter we have an automatic
!        integrator available, qawf.
!
  return
end subroutine
subroutine qag ( f, a, b, epsabs, epsrel, key, result, abserr, neval, ier )

!*****************************************************************************80
!
!! QAG approximates an integral over a finite interval.
!
!  Discussion:
!
!    The routine calculates an approximation RESULT to a definite integral   
!      I = integral of F over (A,B),
!    hopefully satisfying
!      || I - RESULT || <= max ( EPSABS, EPSREL * ||I|| ).
!
!    QAG is a simple globally adaptive integrator using the strategy of 
!    Aind (Piessens, 1973).  It is possible to choose between 6 pairs of
!    Gauss-Kronrod quadrature formulae for the rule evaluation component. 
!    The pairs of high degree of precision are suitable for handling
!    integration difficulties due to a strongly oscillating integrand.
!
!  Author:
!
!    Robert Piessens, Elise de Doncker-Kapenger, 
!    Christian Ueberhuber, David Kahaner
!
!  Reference:
!
!    Robert Piessens, Elise de Doncker-Kapenger, 
!    Christian Ueberhuber, David Kahaner,
!    QUADPACK, a Subroutine Package for Automatic Integration,
!    Springer Verlag, 1983
!
!  Parameters:
!
!    Input, external real F, the name of the function routine, of the form
!      function f ( x )
!      real f
!      real x
!    which evaluates the integrand function.
!
!    Input, real A, B, the limits of integration.
!
!    Input, real EPSABS, EPSREL, the absolute and relative accuracy requested.
!
!    Input, integer KEY, chooses the order of the local integration rule:
!    1,  7 Gauss points, 15 Gauss-Kronrod points,
!    2, 10 Gauss points, 21 Gauss-Kronrod points,
!    3, 15 Gauss points, 31 Gauss-Kronrod points,
!    4, 20 Gauss points, 41 Gauss-Kronrod points,
!    5, 25 Gauss points, 51 Gauss-Kronrod points,
!    6, 30 Gauss points, 61 Gauss-Kronrod points.
!
!    Output, real RESULT, the estimated value of the integral.
!
!    Output, real ABSERR, an estimate of || I - RESULT ||.
!
!    Output, integer NEVAL, the number of times the integral was evaluated.
!
!    Output, integer IER, return code.
!    0, normal and reliable termination of the routine.  It is assumed that the 
!      requested accuracy has been achieved.
!    1, maximum number of subdivisions allowed has been achieved.  One can 
!      allow more subdivisions by increasing the value of LIMIT in QAG. 
!      However, if this yields no improvement it is advised to analyze the
!      integrand to determine the integration difficulties.  If the position
!      of a local difficulty can be determined, such as a singularity or
!      discontinuity within the interval) one will probably gain from 
!      splitting up the interval at this point and calling the integrator 
!      on the subranges.  If possible, an appropriate special-purpose 
!      integrator should be used which is designed for handling the type 
!      of difficulty involved.
!    2, the occurrence of roundoff error is detected, which prevents the
!      requested tolerance from being achieved.
!    3, extremely bad integrand behavior occurs at some points of the
!      integration interval.
!    6, the input is invalid, because EPSABS < 0 and EPSREL < 0.
!
!  Local parameters:
!
!    LIMIT is the maximum number of subintervals allowed in
!    the subdivision process of QAGE.
!
  implicit none

  integer, parameter :: limit = 500

  real a
  real abserr
  real alist(limit)
  real b
  real blist(limit)
  real elist(limit)
  real epsabs
  real epsrel
  procedure(scalar_func) f
  integer ier
  integer iord(limit)
  integer key
  integer last
  integer neval
  real result
  real rlist(limit)

  call qage ( f, a, b, epsabs, epsrel, key, limit, result, abserr, neval, &
    ier, alist, blist, rlist, elist, iord, last )

  return
end subroutine
subroutine qage ( f, a, b, epsabs, epsrel, key, limit, result, abserr, neval, &
  ier, alist, blist, rlist, elist, iord, last )

!*****************************************************************************80
!
!! QAGE estimates a definite integral.
!
!  Discussion:
!
!    The routine calculates an approximation RESULT to a definite integral   
!      I = integral of F over (A,B),
!    hopefully satisfying
!      || I - RESULT || <= max ( EPSABS, EPSREL * ||I|| ).
!
!  Author:
!
!    Robert Piessens, Elise de Doncker-Kapenger, 
!    Christian Ueberhuber, David Kahaner
!
!  Reference:
!
!    Robert Piessens, Elise de Doncker-Kapenger, 
!    Christian Ueberhuber, David Kahaner,
!    QUADPACK, a Subroutine Package for Automatic Integration,
!    Springer Verlag, 1983
!
!  Parameters:
!
!    Input, external real F, the name of the function routine, of the form
!      function f ( x )
!      real f
!      real x
!    which evaluates the integrand function.
!
!    Input, real A, B, the limits of integration.
!
!    Input, real EPSABS, EPSREL, the absolute and relative accuracy requested.
!
!    Input, integer KEY, chooses the order of the local integration rule:
!    1,  7 Gauss points, 15 Gauss-Kronrod points,
!    2, 10 Gauss points, 21 Gauss-Kronrod points,
!    3, 15 Gauss points, 31 Gauss-Kronrod points,
!    4, 20 Gauss points, 41 Gauss-Kronrod points,
!    5, 25 Gauss points, 51 Gauss-Kronrod points,
!    6, 30 Gauss points, 61 Gauss-Kronrod points.
!
!    Input, integer LIMIT, the maximum number of subintervals that
!    can be used.
!
!    Output, real RESULT, the estimated value of the integral.
!
!    Output, real ABSERR, an estimate of || I - RESULT ||.
!
!    Output, integer NEVAL, the number of times the integral was evaluated.
!
!    Output, integer IER, return code.
!    0, normal and reliable termination of the routine.  It is assumed that the 
!      requested accuracy has been achieved.
!    1, maximum number of subdivisions allowed has been achieved.  One can 
!      allow more subdivisions by increasing the value of LIMIT in QAG. 
!      However, if this yields no improvement it is advised to analyze the
!      integrand to determine the integration difficulties.  If the position
!      of a local difficulty can be determined, such as a singularity or
!      discontinuity within the interval) one will probably gain from 
!      splitting up the interval at this point and calling the integrator 
!      on the subranges.  If possible, an appropriate special-purpose 
!      integrator should be used which is designed for handling the type 
!      of difficulty involved.
!    2, the occurrence of roundoff error is detected, which prevents the
!      requested tolerance from being achieved.
!    3, extremely bad integrand behavior occurs at some points of the
!      integration interval.
!    6, the input is invalid, because EPSABS < 0 and EPSREL < 0.
!
!    Workspace, real ALIST(LIMIT), BLIST(LIMIT), contains in entries 1 
!    through LAST the left and right ends of the partition subintervals.
!
!    Workspace, real RLIST(LIMIT), contains in entries 1 through LAST
!    the integral approximations on the subintervals.
!
!    Workspace, real ELIST(LIMIT), contains in entries 1 through LAST
!    the absolute error estimates on the subintervals.
!
!    Output, integer IORD(LIMIT), the first K elements of which are pointers 
!    to the error estimates over the subintervals, such that
!    elist(iord(1)), ..., elist(iord(k)) form a decreasing sequence, with
!    k = last if last <= (limit/2+2), and k = limit+1-last otherwise.
!
!    Output, integer LAST, the number of subintervals actually produced 
!    in the subdivision process.
!
!  Local parameters:
!
!    alist     - list of left end points of all subintervals
!                       considered up to now
!    blist     - list of right end points of all subintervals
!                       considered up to now
!    elist(i)  - error estimate applying to rlist(i)
!    maxerr    - pointer to the interval with largest error estimate
!    errmax    - elist(maxerr)
!    area      - sum of the integrals over the subintervals
!    errsum    - sum of the errors over the subintervals
!    errbnd    - requested accuracy max(epsabs,epsrel*abs(result))
!    *****1    - variable for the left subinterval
!    *****2    - variable for the right subinterval
!    last      - index for subdivision
!
  implicit none

  integer limit

  real a
  real abserr
  real alist(limit)
  real area
  real area1
  real area12
  real area2
  real a1
  real a2
  real b
  real blist(limit)
  real b1
  real b2
  real c
  real defabs
  real defab1
  real defab2
  real elist(limit)
  real epsabs
  real epsrel
  real errbnd
  real errmax
  real error1
  real error2
  real erro12
  real errsum
  procedure(scalar_func) :: f
  integer ier
  integer iord(limit)
  integer iroff1
  integer iroff2
  integer key
  integer keyf
  integer last
  integer maxerr
  integer neval
  integer nrmax
  real resabs
  real result
  real rlist(limit)
!
!  Test on validity of parameters.
!
  ier = 0
  neval = 0
  last = 0
  result = 0.0e+00
  abserr = 0.0e+00
  alist(1) = a
  blist(1) = b
  rlist(1) = 0.0e+00
  elist(1) = 0.0e+00
  iord(1) = 0

  if ( epsabs < 0.0e+00 .and. epsrel < 0.0e+00 ) then
    ier = 6
    return
  end if
!
!  First approximation to the integral.
!
  keyf = key
  keyf = max ( keyf, 1 )
  keyf = min ( keyf, 6 )

  c = keyf
  neval = 0

  if ( keyf == 1 ) then
    call qk15 ( f, a, b, result, abserr, defabs, resabs )
  else if ( keyf == 2 ) then
    call qk21 ( f, a, b, result, abserr, defabs, resabs )
  else if ( keyf == 3 ) then
    call qk31 ( f, a, b, result, abserr, defabs, resabs )
  else if ( keyf == 4 ) then
    call qk41 ( f, a, b, result, abserr, defabs, resabs )
  else if ( keyf == 5 ) then
    call qk51 ( f, a, b, result, abserr, defabs, resabs )
  else if ( keyf == 6 ) then
    call qk61 ( f, a, b, result, abserr, defabs, resabs )
  end if

  last = 1
  rlist(1) = result
  elist(1) = abserr
  iord(1) = 1
!
!  Test on accuracy.
!
  errbnd = max ( epsabs, epsrel * abs ( result ) )

  if ( abserr <= 5.0e+01 * epsilon ( defabs ) * defabs .and. &
    errbnd < abserr ) then
    ier = 2
  end if

  if ( limit == 1 ) then
    ier = 1
  end if

  if ( ier /= 0 .or. &
    ( abserr <= errbnd .and. not_exactly_equal(abserr, resabs)) .or. &
    is_zero(abserr)) then

    if ( keyf /= 1 ) then
      neval = (10*keyf+1) * (2*neval+1)
    else
      neval = 30 * neval + 15
    end if

    return

  end if
!
!  Initialization.
!
  errmax = abserr
  maxerr = 1
  area = result
  errsum = abserr
  nrmax = 1
  iroff1 = 0
  iroff2 = 0

  do last = 2, limit
!
!  Bisect the subinterval with the largest error estimate.
!
    a1 = alist(maxerr)
    b1 = 0.5E+00 * ( alist(maxerr) + blist(maxerr) )
    a2 = b1
    b2 = blist(maxerr)

    if ( keyf == 1 ) then
      call qk15 ( f, a1, b1, area1, error1, resabs, defab1 )
    else if ( keyf == 2 ) then
      call qk21 ( f, a1, b1, area1, error1, resabs, defab1 )
    else if ( keyf == 3 ) then
      call qk31 ( f, a1, b1, area1, error1, resabs, defab1 )
    else if ( keyf == 4 ) then
      call qk41 ( f, a1, b1, area1, error1, resabs, defab1)
    else if ( keyf == 5 ) then
      call qk51 ( f, a1, b1, area1, error1, resabs, defab1 )
    else if ( keyf == 6 ) then
      call qk61 ( f, a1, b1, area1, error1, resabs, defab1 )
    end if

    if ( keyf == 1 ) then
      call qk15 ( f, a2, b2, area2, error2, resabs, defab2 )
    else if ( keyf == 2 ) then
      call qk21 ( f, a2, b2, area2, error2, resabs, defab2 )
    else if ( keyf == 3 ) then
      call qk31 ( f, a2, b2, area2, error2, resabs, defab2 )
    else if ( keyf == 4 ) then
      call qk41 ( f, a2, b2, area2, error2, resabs, defab2 )
    else if ( keyf == 5 ) then
      call qk51 ( f, a2, b2, area2, error2, resabs, defab2 )
    else if ( keyf == 6 ) then
      call qk61 ( f, a2, b2, area2, error2, resabs, defab2 )
    end if
!
!  Improve previous approximations to integral and error and
!  test for accuracy.
!
    neval = neval + 1
    area12 = area1 + area2
    erro12 = error1 + error2
    errsum = errsum + erro12 - errmax
    area = area + area12 - rlist(maxerr)

    if (not_exactly_equal(defab1, error1) .and. not_exactly_equal(defab2, error2)) then

      if ( abs ( rlist(maxerr) - area12 ) <= 1.0e-05 * abs ( area12 ) &
        .and. 9.9e-01 * errmax <= erro12 ) then
        iroff1 = iroff1 + 1
      end if

      if ( 10 < last .and. errmax < erro12 ) then
        iroff2 = iroff2 + 1
      end if

    end if

    rlist(maxerr) = area1
    rlist(last) = area2
    errbnd = max ( epsabs, epsrel * abs ( area ) )
!
!  Test for roundoff error and eventually set error flag.
!
    if ( errbnd < errsum ) then

      if ( 6 <= iroff1 .or. 20 <= iroff2 ) then
        ier = 2
      end if
!
!  Set error flag in the case that the number of subintervals
!  equals limit.
!
      if ( last == limit ) then
        ier = 1
      end if
!
!  Set error flag in the case of bad integrand behavior
!  at a point of the integration range.
!
      if ( max ( abs ( a1 ), abs ( b2 ) ) <= ( 1.0e+00 + c * 1.0e+03 * &
        epsilon ( a1 ) ) * ( abs ( a2 ) + 1.0e+04 * tiny ( a2 ) ) ) then
        ier = 3
      end if

    end if
!
!  Append the newly-created intervals to the list.
!
    if ( error2 <= error1 ) then
      alist(last) = a2
      blist(maxerr) = b1
      blist(last) = b2
      elist(maxerr) = error1
      elist(last) = error2
    else
      alist(maxerr) = a2
      alist(last) = a1
      blist(last) = b1
      rlist(maxerr) = area2
      rlist(last) = area1
      elist(maxerr) = error2
      elist(last) = error1
    end if
!
!  Call QSORT to maintain the descending ordering
!  in the list of error estimates and select the subinterval
!  with the largest error estimate (to be bisected next).
!
    call qsort ( limit, last, maxerr, errmax, elist, iord, nrmax )
 
    if ( ier /= 0 .or. errsum <= errbnd ) then
      exit
    end if

  end do
!
!  Compute final result.
!
  result = sum ( rlist(1:last) )

  abserr = errsum

  if ( keyf /= 1 ) then
    neval = ( 10 * keyf + 1 ) * ( 2 * neval + 1 )
  else
    neval = 30 * neval + 15
  end if

  return
end subroutine
subroutine qagi ( f, bound, inf, epsabs, epsrel, result, abserr, neval, ier )

!*****************************************************************************80
!
!! QAGI estimates an integral over a semi-infinite or infinite interval.
!
!  Discussion:
!
!    The routine calculates an approximation RESULT to a definite integral   
!      I = integral of F over (A, +Infinity), 
!    or 
!      I = integral of F over (-Infinity,A)
!    or 
!      I = integral of F over (-Infinity,+Infinity),
!    hopefully satisfying
!      || I - RESULT || <= max ( EPSABS, EPSREL * ||I|| ).
!
!  Author:
!
!    Robert Piessens, Elise de Doncker-Kapenger, 
!    Christian Ueberhuber, David Kahaner
!
!  Reference:
!
!    Robert Piessens, Elise de Doncker-Kapenger, 
!    Christian Ueberhuber, David Kahaner,
!    QUADPACK, a Subroutine Package for Automatic Integration,
!    Springer Verlag, 1983
!
!  Parameters:
!
!    Input, external real F, the name of the function routine, of the form
!      function f ( x )
!      real f
!      real x
!    which evaluates the integrand function.
!
!    Input, real BOUND, the value of the finite endpoint of the integration
!    range, if any, that is, if INF is 1 or -1.
!
!    Input, integer INF, indicates the type of integration range.
!    1:  (  BOUND,    +Infinity),
!    -1: ( -Infinity,  BOUND),
!    2:  ( -Infinity, +Infinity).
!
!    Input, real EPSABS, EPSREL, the absolute and relative accuracy requested.
!
!    Output, real RESULT, the estimated value of the integral.
!
!    Output, real ABSERR, an estimate of || I - RESULT ||.
!
!    Output, integer NEVAL, the number of times the integral was evaluated.
!
!    Output, integer IER, error indicator.
!    0, normal and reliable termination of the routine.  It is assumed that 
!      the requested accuracy has been achieved.
!    > 0,  abnormal termination of the routine.  The estimates for result
!      and error are less reliable.  It is assumed that the requested
!      accuracy has not been achieved.
!    1, maximum number of subdivisions allowed has been achieved.  One can 
!      allow more subdivisions by increasing the data value of LIMIT in QAGI
!      (and taking the according dimension adjustments into account).
!      However, if this yields no improvement it is advised to analyze the
!      integrand in order to determine the integration difficulties.  If the
!      position of a local difficulty can be determined (e.g. singularity,
!      discontinuity within the interval) one will probably gain from
!      splitting up the interval at this point and calling the integrator 
!      on the subranges.  If possible, an appropriate special-purpose 
!      integrator should be used, which is designed for handling the type
!      of difficulty involved.
!    2, the occurrence of roundoff error is detected, which prevents the
!      requested tolerance from being achieved.  The error may be
!      under-estimated.
!    3, extremely bad integrand behavior occurs at some points of the
!      integration interval.
!    4, the algorithm does not converge.  Roundoff error is detected in the
!      extrapolation table.  It is assumed that the requested tolerance
!      cannot be achieved, and that the returned result is the best which 
!      can be obtained.
!    5, the integral is probably divergent, or slowly convergent.  It must 
!      be noted that divergence can occur with any other value of IER.
!    6, the input is invalid, because INF /= 1 and INF /= -1 and INF /= 2, or
!      epsabs < 0 and epsrel < 0.  result, abserr, neval are set to zero.
!
!  Local parameters:
!
!            the dimension of rlist2 is determined by the value of
!            limexp in QEXTR.
!
!           alist     - list of left end points of all subintervals
!                       considered up to now
!           blist     - list of right end points of all subintervals
!                       considered up to now
!           rlist(i)  - approximation to the integral over
!                       (alist(i),blist(i))
!           rlist2    - array of dimension at least (limexp+2),
!                       containing the part of the epsilon table
!                       which is still needed for further computations
!           elist(i)  - error estimate applying to rlist(i)
!           maxerr    - pointer to the interval with largest error
!                       estimate
!           errmax    - elist(maxerr)
!           erlast    - error on the interval currently subdivided
!                       (before that subdivision has taken place)
!           area      - sum of the integrals over the subintervals
!           errsum    - sum of the errors over the subintervals
!           errbnd    - requested accuracy max(epsabs,epsrel*
!                       abs(result))
!           *****1    - variable for the left subinterval
!           *****2    - variable for the right subinterval
!           last      - index for subdivision
!           nres      - number of calls to the extrapolation routine
!           numrl2    - number of elements currently in rlist2. if an
!                       appropriate approximation to the compounded
!                       integral has been obtained, it is put in
!                       rlist2(numrl2) after numrl2 has been increased
!                       by one.
!           small     - length of the smallest interval considered up
!                       to now, multiplied by 1.5
!           erlarg    - sum of the errors over the intervals larger
!                       than the smallest interval considered up to now
!           extrap    - logical variable denoting that the routine
!                       is attempting to perform extrapolation. i.e.
!                       before subdividing the smallest interval we
!                       try to decrease the value of erlarg.
!           noext     - logical variable denoting that extrapolation
!                       is no longer allowed (true-value)
!
  implicit none

  integer, parameter :: limit = 500

  real abseps
  real abserr
  real alist(limit)
  real area
  real area1
  real area12
  real area2
  real a1
  real a2
  real blist(limit)
  real boun
  real bound
  real b1
  real b2
  real correc
  real defabs
  real defab1
  real defab2
  real dres
  real elist(limit)
  real epsabs
  real epsrel
  real erlarg
  real erlast
  real errbnd
  real errmax
  real error1
  real error2
  real erro12
  real errsum
  real ertest
  logical extrap
  procedure(scalar_func) :: f
  integer id
  integer ier
  integer ierro
  integer inf
  integer iord(limit)
  integer iroff1
  integer iroff2
  integer iroff3
  integer jupbnd
  integer k
  integer ksgn
  integer ktmin
  integer last
  integer maxerr
  integer neval
  logical noext
  integer nres
  integer nrmax
  integer numrl2
  real resabs
  real reseps
  real result
  real res3la(3)
  real rlist(limit)
  real rlist2(52)
  real small
!
!  Test on validity of parameters.
!
  ier = 0
  neval = 0
  last = 0
  result = 0.0e+00
  abserr = 0.0e+00
  alist(1) = 0.0e+00
  blist(1) = 1.0e+00
  rlist(1) = 0.0e+00
  elist(1) = 0.0e+00
  iord(1) = 0

  if ( epsabs < 0.0e+00 .and. epsrel < 0.0e+00 ) then
    ier = 6
    return
  end if
!
!  First approximation to the integral.
!
!  Determine the interval to be mapped onto (0,1).
!  If INF = 2 the integral is computed as i = i1+i2, where
!  i1 = integral of f over (-infinity,0),
!  i2 = integral of f over (0,+infinity).
!
  if ( inf == 2 ) then
    boun = 0.0e+00
  else
    boun = bound
  end if

  call qk15i ( f, boun, inf, 0.0e+00, 1.0e+00, result, abserr, defabs, resabs )
!
!  Test on accuracy.
!
  last = 1
  rlist(1) = result
  elist(1) = abserr
  iord(1) = 1
  dres = abs ( result )
  errbnd = max ( epsabs, epsrel * dres )

  if ( abserr <= 100.0E+00 * epsilon ( defabs ) * defabs .and. &
    errbnd < abserr ) then
    ier = 2
  end if

  if ( limit == 1 ) then
    ier = 1
  end if

  if (ier /= 0 .or. (abserr <= errbnd .and. not_exactly_equal(abserr, resabs)) &
       .or. is_zero(abserr)) go to 130
!
!  Initialization.
!
  rlist2(1) = result
  errmax = abserr
  maxerr = 1
  area = result
  errsum = abserr
  abserr = huge ( abserr )
  nrmax = 1
  nres = 0
  ktmin = 0
  numrl2 = 2
  extrap = .false.
  noext = .false.
  ierro = 0
  iroff1 = 0
  iroff2 = 0
  iroff3 = 0

  if ( ( 1.0e+00 - 5.0e+01 * epsilon ( defabs ) ) * defabs <= dres ) then
    ksgn = 1
  else
    ksgn = -1
  end if

  do last = 2, limit
!
!  Bisect the subinterval with nrmax-th largest error estimate.
!
    a1 = alist(maxerr)
    b1 = 5.0e-01 * ( alist(maxerr) + blist(maxerr) )
    a2 = b1
    b2 = blist(maxerr)
    erlast = errmax
    call qk15i ( f, boun, inf, a1, b1, area1, error1, resabs, defab1 )
    call qk15i ( f, boun, inf, a2, b2, area2, error2, resabs, defab2 )
!
!  Improve previous approximations to integral and error
!  and test for accuracy.
!
    area12 = area1 + area2
    erro12 = error1 + error2
    errsum = errsum + erro12 - errmax
    area = area + area12 - rlist(maxerr)

    if (not_exactly_equal(defab1, error1) .and. not_exactly_equal(defab2, error2)) then

      if ( abs ( rlist(maxerr) - area12 ) <= 1.0e-05 * abs ( area12 ) &
        .and. 9.9e-01 * errmax <= erro12 ) then

        if ( extrap ) then
          iroff2 = iroff2 + 1
        end if

        if ( .not. extrap ) then
          iroff1 = iroff1 + 1
        end if

      end if

      if ( 10 < last .and. errmax < erro12 ) then
        iroff3 = iroff3 + 1
      end if

    end if

    rlist(maxerr) = area1
    rlist(last) = area2
    errbnd = max ( epsabs, epsrel * abs ( area ) )
!
!  Test for roundoff error and eventually set error flag.
!
    if ( 10 <= iroff1 + iroff2 .or. 20 <= iroff3 ) then
      ier = 2
    end if

    if ( 5 <= iroff2 ) then
      ierro = 3
    end if
!
!  Set error flag in the case that the number of subintervals equals LIMIT.
!
    if ( last == limit ) then
      ier = 1
    end if
!
!  Set error flag in the case of bad integrand behavior
!  at some points of the integration range.
!
    if ( max ( abs(a1), abs(b2) ) <= (1.0e+00 + 1.0e+03 * epsilon ( a1 ) ) * &
    ( abs(a2) + 1.0e+03 * tiny ( a2 ) )) then
      ier = 4
    end if
!
!  Append the newly-created intervals to the list.
!
    if ( error2 <= error1 ) then
      alist(last) = a2
      blist(maxerr) = b1
      blist(last) = b2
      elist(maxerr) = error1
      elist(last) = error2
    else
      alist(maxerr) = a2
      alist(last) = a1
      blist(last) = b1
      rlist(maxerr) = area2
      rlist(last) = area1
      elist(maxerr) = error2
      elist(last) = error1
    end if
!
!  Call QSORT to maintain the descending ordering
!  in the list of error estimates and select the subinterval
!  with NRMAX-th largest error estimate (to be bisected next).
!
    call qsort ( limit, last, maxerr, errmax, elist, iord, nrmax )

    if ( errsum <= errbnd ) go to 115

    if ( ier /= 0 ) then
      exit
    end if

    if ( last == 2 ) then
      small = 3.75e-01
      erlarg = errsum
      ertest = errbnd
      rlist2(2) = area
      cycle
    end if

    if ( noext ) then
      cycle
    end if

    erlarg = erlarg - erlast

    if ( small < abs ( b1 - a1 ) ) then
      erlarg = erlarg + erro12
    end if
!
!  Test whether the interval to be bisected next is the
!  smallest interval.
!
    if ( .not. extrap ) then

      if ( small < abs ( blist(maxerr) - alist(maxerr) ) ) then
        cycle
      end if

      extrap = .true.
      nrmax = 2

    end if

    if ( ierro == 3 .or. erlarg <= ertest ) then
      go to 60
    end if
!
!  The smallest interval has the largest error.
!  before bisecting decrease the sum of the errors over the
!  larger intervals (erlarg) and perform extrapolation.
!
    id = nrmax
    jupbnd = last

    if ( (2+limit/2) < last ) then
      jupbnd = limit + 3 - last
    end if

    do k = id, jupbnd
      maxerr = iord(nrmax)
      errmax = elist(maxerr)
      if ( small < abs ( blist(maxerr) - alist(maxerr) ) ) then
        go to 90
      end if
      nrmax = nrmax + 1
    end do
!
!  Extrapolate.
!
60  continue

    numrl2 = numrl2 + 1
    rlist2(numrl2) = area
    call qextr ( numrl2, rlist2, reseps, abseps, res3la, nres ) 
    ktmin = ktmin+1

    if ( 5 < ktmin .and. abserr < 1.0e-03 * errsum ) then
      ier = 5
    end if

    if ( abseps < abserr ) then

      ktmin = 0
      abserr = abseps
      result = reseps
      correc = erlarg
      ertest = max ( epsabs, epsrel * abs(reseps) )

      if ( abserr <= ertest ) then
        exit
      end if

    end if
!
!  Prepare bisection of the smallest interval.
!
    if ( numrl2 == 1 ) then
      noext = .true.
    end if

    if ( ier == 5 ) then
      exit
    end if

    maxerr = iord(1)
    errmax = elist(maxerr)
    nrmax = 1
    extrap = .false.
    small = small * 5.0e-01
    erlarg = errsum

90  continue

  end do
!
!  Set final result and error estimate.
!
  if (exactly_equal(abserr, huge(abserr))) go to 115

  if ( ( ier + ierro ) == 0 ) then
    go to 110
  end if

  if ( ierro == 3 ) then
    abserr = abserr + correc
  end if

  if ( ier == 0 ) then
    ier = 3
  end if

  if (is_not_zero(result) .and. is_not_zero(area)) then
    go to 105
  end if

  if ( errsum < abserr ) then
    go to 115
  end if

  if (is_zero(area)) go to 130

  go to 110

105 continue

  if ( errsum / abs ( area ) < abserr / abs ( result )  ) then
    go to 115
  end if
!
!  Test on divergence
!
110 continue

  if ( ksgn == (-1) .and. &
  max ( abs(result), abs(area) ) <=  defabs * 1.0e-02) go to 130

  if ( 1.0e-02 > (result/area) .or. &
    (result/area) > 1.0e+02 .or. &
    errsum > abs(area)) then
    ier = 6
  end if

  go to 130
!
!  Compute global integral sum.
!
  115 continue

  result = sum ( rlist(1:last) )

  abserr = errsum
  130 continue

  neval = 30 * last - 15
  if ( inf == 2 ) then
    neval = 2 * neval
  end if

  if ( 2 < ier ) then
    ier = ier - 1
  end if

  return
end subroutine
subroutine qagp ( f, a, b, npts2, points, epsabs, epsrel, result, abserr, &
  neval, ier )

!*****************************************************************************80
!
!! QAGP computes a definite integral.
!
!  Discussion:
!
!    The routine calculates an approximation RESULT to a definite integral   
!      I = integral of F over (A,B),
!    hopefully satisfying
!      || I - RESULT || <= max ( EPSABS, EPSREL * ||I|| ).
!
!    Interior break points of the integration interval,
!    where local difficulties of the integrand may occur, such as
!    singularities or discontinuities, are provided by the user.
!
!  Author:
!
!    Robert Piessens, Elise de Doncker-Kapenger, 
!    Christian Ueberhuber, David Kahaner
!
!  Reference:
!
!    Robert Piessens, Elise de Doncker-Kapenger, 
!    Christian Ueberhuber, David Kahaner,
!    QUADPACK, a Subroutine Package for Automatic Integration,
!    Springer Verlag, 1983
!
!  Parameters:
!
!    Input, external real F, the name of the function routine, of the form
!      function f ( x )
!      real f
!      real x
!    which evaluates the integrand function.
!
!    Input, real A, B, the limits of integration.
!
!    Input, integer NPTS2, the number of user-supplied break points within 
!    the integration range, plus 2.  NPTS2 must be at least 2.
!
!    Input/output, real POINTS(NPTS2), contains the user provided interior
!    breakpoints in entries 1 through NPTS2-2.  If these points are not
!    in ascending order on input, they will be sorted.
!
!    Input, real EPSABS, EPSREL, the absolute and relative accuracy requested.
!
!    Output, real RESULT, the estimated value of the integral.
!
!    Output, real ABSERR, an estimate of || I - RESULT ||.
!
!    Output, integer NEVAL, the number of times the integral was evaluated.
!
!    Output, integer IER, return flag.
!                     ier = 0 normal and reliable termination of the
!                             routine. it is assumed that the requested
!                             accuracy has been achieved.
!                     ier > 0 abnormal termination of the routine.
!                             the estimates for integral and error are
!                             less reliable. it is assumed that the
!                             requested accuracy has not been achieved.
!                     ier = 1 maximum number of subdivisions allowed
!                             has been achieved. one can allow more
!                             subdivisions by increasing the data value
!                             of limit in qagp(and taking the according
!                             dimension adjustments into account).
!                             however, if this yields no improvement
!                             it is advised to analyze the integrand
!                             in order to determine the integration
!                             difficulties. if the position of a local
!                             difficulty can be determined (i.e.
!                             singularity, discontinuity within the
!                             interval), it should be supplied to the
!                             routine as an element of the vector
!                             points. if necessary, an appropriate
!                             special-purpose integrator must be used,
!                             which is designed for handling the type
!                             of difficulty involved.
!                         = 2 the occurrence of roundoff error is
!                             detected, which prevents the requested
!                             tolerance from being achieved.
!                             the error may be under-estimated.
!                         = 3 extremely bad integrand behavior occurs
!                             at some points of the integration
!                             interval.
!                         = 4 the algorithm does not converge. roundoff
!                             error is detected in the extrapolation
!                             table. it is presumed that the requested
!                             tolerance cannot be achieved, and that
!                             the returned result is the best which
!                             can be obtained.
!                         = 5 the integral is probably divergent, or
!                             slowly convergent. it must be noted that
!                             divergence can occur with any other value
!                             of ier > 0.
!                         = 6 the input is invalid because
!                             npts2 < 2 or
!                             break points are specified outside
!                             the integration range or
!                             epsabs < 0 and epsrel < 0,
!                             or limit < npts2.
!                             result, abserr, neval are set to zero.
!
!  Local parameters:
!
!            the dimension of rlist2 is determined by the value of
!            limexp in QEXTR (rlist2 should be of dimension
!            (limexp+2) at least).
!
!           alist     - list of left end points of all subintervals
!                       considered up to now
!           blist     - list of right end points of all subintervals
!                       considered up to now
!           rlist(i)  - approximation to the integral over
!                       (alist(i),blist(i))
!           rlist2    - array of dimension at least limexp+2
!                       containing the part of the epsilon table which
!                       is still needed for further computations
!           elist(i)  - error estimate applying to rlist(i)
!           maxerr    - pointer to the interval with largest error
!                       estimate
!           errmax    - elist(maxerr)
!           erlast    - error on the interval currently subdivided
!                       (before that subdivision has taken place)
!           area      - sum of the integrals over the subintervals
!           errsum    - sum of the errors over the subintervals
!           errbnd    - requested accuracy max(epsabs,epsrel*
!                       abs(result))
!           *****1    - variable for the left subinterval
!           *****2    - variable for the right subinterval
!           last      - index for subdivision
!           nres      - number of calls to the extrapolation routine
!           numrl2    - number of elements in rlist2. if an appropriate
!                       approximation to the compounded integral has
!                       obtained, it is put in rlist2(numrl2) after
!                       numrl2 has been increased by one.
!           erlarg    - sum of the errors over the intervals larger
!                       than the smallest interval considered up to now
!           extrap    - logical variable denoting that the routine
!                       is attempting to perform extrapolation. i.e.
!                       before subdividing the smallest interval we
!                       try to decrease the value of erlarg.
!           noext     - logical variable denoting that extrapolation is
!                       no longer allowed (true-value)
!
  implicit none

  integer, parameter :: limit = 500

  real a
  real abseps
  real abserr
  real alist(limit)
  real area
  real area1
  real area12
  real area2
  real a1
  real a2
  real b
  real blist(limit)
  real b1
  real b2
  real correc
  real defabs
  real defab1
  real defab2
  real dres
  real elist(limit)
  real epsabs
  real epsrel
  real erlarg
  real erlast
  real errbnd
  real errmax
  real error1
  real erro12
  real error2
  real errsum
  real ertest
  logical extrap
  procedure(scalar_func) :: f
  integer i
  integer id
  integer ier
  integer ierro
  integer ind1
  integer ind2
  integer iord(limit)
  integer iroff1
  integer iroff2
  integer iroff3
  integer j
  integer jlow
  integer jupbnd
  integer k
  integer ksgn
  integer ktmin
  integer last
  integer levcur
  integer level(limit)
  integer levmax
  integer maxerr
  integer ndin(40)
  integer neval
  integer nint
  logical noext
  integer npts
  integer npts2
  integer nres
  integer nrmax
  integer numrl2
  real points(40)
  real pts(40)
  real resa
  real resabs
  real reseps
  real result
  real res3la(3)
  real rlist(limit)
  real rlist2(52)
  real sign
  real temp
!
!  Test on validity of parameters.
!
  ier = 0
  neval = 0
  last = 0
  result = 0.0e+00
  abserr = 0.0e+00
  alist(1) = a
  blist(1) = b
  rlist(1) = 0.0e+00
  elist(1) = 0.0e+00
  iord(1) = 0
  level(1) = 0
  npts = npts2 - 2

  if ( npts2 < 2 ) then
    ier = 6
    return
  else if ( limit <= npts .or. ( epsabs < 0.0e+00 .and. &
    epsrel < 0.0e+00) ) then
    ier = 6
    return
  end if
!
!  If any break points are provided, sort them into an
!  ascending sequence.
!
  if ( b < a ) then
    sign = -1.0e+00
  else
    sign = +1.0E+00
  end if

  pts(1) = min ( a, b )

  do i = 1, npts
    pts(i+1) = points(i)
  end do

  pts(npts+2) = max ( a, b )
  nint = npts+1
  a1 = pts(1)

  if ( npts /= 0 ) then

    do i = 1, nint
      do j = i+1, nint+1
        if ( pts(j) < pts(i) ) then
          temp   = pts(i)
          pts(i) = pts(j)
          pts(j) = temp
        end if
      end do
    end do

    if (not_exactly_equal(pts(1), min(a, b)) .or. not_exactly_equal(pts(nint+1), max(a, b))) then
      ier = 6
      return
    end if

  end if
!
!  Compute first integral and error approximations.
!
  resabs = 0.0e+00

  do i = 1, nint

    b1 = pts(i+1)
    call qk21 ( f, a1, b1, area1, error1, defabs, resa )
    abserr = abserr + error1
    result = result + area1
    ndin(i) = 0

    if (exactly_equal(error1, resa) .and. is_not_zero(error1)) then
      ndin(i) = 1
    end if

    resabs = resabs + defabs
    level(i) = 0
    elist(i) = error1
    alist(i) = a1
    blist(i) = b1
    rlist(i) = area1
    iord(i) = i
    a1 = b1

  end do

  errsum = 0.0e+00

  do i = 1, nint
    if ( ndin(i) == 1 ) then
      elist(i) = abserr
    end if
    errsum = errsum + elist(i)
  end do
!
!  Test on accuracy.
!
  last = nint
  neval = 21 * nint
  dres = abs ( result )
  errbnd = max ( epsabs, epsrel * dres )

  if ( abserr <= 1.0e+02 * epsilon ( resabs ) * resabs .and. &
    abserr > errbnd ) then
    ier = 2
  end if

  if ( nint /= 1 ) then

    do i = 1, npts

      jlow = i+1
      ind1 = iord(i)

      do j = jlow, nint
        ind2 = iord(j)
        if ( elist(ind1) <= elist(ind2) ) then
          ind1 = ind2
          k = j
        end if
      end do

      if ( ind1 /= iord(i) ) then
        iord(k) = iord(i)
        iord(i) = ind1
      end if

    end do

    if ( limit < npts2 ) then
      ier = 1
    end if

  end if

  if ( ier /= 0 .or. abserr <= errbnd ) then
    return
  end if
!
!  Initialization
!
  rlist2(1) = result
  maxerr = iord(1)
  errmax = elist(maxerr)
  area = result
  nrmax = 1
  nres = 0
  numrl2 = 1
  ktmin = 0
  extrap = .false.
  noext = .false.
  erlarg = errsum
  ertest = errbnd
  levmax = 1
  iroff1 = 0
  iroff2 = 0
  iroff3 = 0
  ierro = 0
  abserr = huge ( abserr )

  if ( dres >= ( 1.0e+00 - 0.5E+00 * epsilon ( resabs ) ) * resabs ) then
    ksgn = 1
  else
    ksgn = -1
  end if

  do last = npts2, limit
!
!  Bisect the subinterval with the NRMAX-th largest error estimate.
!
    levcur = level(maxerr) + 1
    a1 = alist(maxerr)
    b1 = 0.5E+00 * ( alist(maxerr) + blist(maxerr) )
    a2 = b1
    b2 = blist(maxerr)
    erlast = errmax
    call qk21 ( f, a1, b1, area1, error1, resa, defab1 )
    call qk21 ( f, a2, b2, area2, error2, resa, defab2 )
!
!  Improve previous approximations to integral and error
!  and test for accuracy.
!
    neval = neval + 42
    area12 = area1 + area2
    erro12 = error1 + error2
    errsum = errsum + erro12 -errmax
    area = area + area12 - rlist(maxerr)

    if (not_exactly_equal(defab1, error1) .and. not_exactly_equal(defab2, error2)) then

      if ( abs ( rlist ( maxerr ) - area12 ) <= 1.0e-05 * abs(area12) .and. &
        erro12 >= 9.9e-01 * errmax ) then

        if ( extrap ) then
          iroff2 = iroff2+1
        else
          iroff1 = iroff1+1
        end if

      end if

      if ( last > 10 .and. erro12 > errmax ) then
        iroff3 = iroff3 + 1
      end if

    end if

    level(maxerr) = levcur
    level(last) = levcur
    rlist(maxerr) = area1
    rlist(last) = area2
    errbnd = max ( epsabs, epsrel * abs ( area ) )
!
!  Test for roundoff error and eventually set error flag.
!
    if ( 10 <= iroff1 + iroff2 .or. 20 <= iroff3 ) then
      ier = 2
    end if

    if ( 5 <= iroff2 ) then
      ierro = 3
    end if
!
!  Set error flag in the case that the number of subintervals
!  equals limit.
!
    if ( last == limit ) then
      ier = 1
    end if
!
!  Set error flag in the case of bad integrand behavior
!  at a point of the integration range
!
    if ( max ( abs(a1), abs(b2)) <= ( 1.0e+00 + 1.0e+03 * epsilon ( a1 ) )* &
    ( abs ( a2 ) + 1.0e+03 * tiny ( a2 ) ) ) then
      ier = 4
    end if
!
!  Append the newly-created intervals to the list.
!
    if ( error2 <= error1 ) then
      alist(last) = a2
      blist(maxerr) = b1
      blist(last) = b2
      elist(maxerr) = error1
      elist(last) = error2
    else
      alist(maxerr) = a2
      alist(last) = a1
      blist(last) = b1
      rlist(maxerr) = area2
      rlist(last) = area1
      elist(maxerr) = error2
      elist(last) = error1
    end if
!
!  Call QSORT to maintain the descending ordering
!  in the list of error estimates and select the subinterval
!  with nrmax-th largest error estimate (to be bisected next).
!
    call qsort ( limit, last, maxerr, errmax, elist, iord, nrmax )

    if ( errsum <= errbnd ) then
      go to 190
    end if

    if ( ier /= 0 ) then
      exit
    end if

    if ( noext ) then
      cycle
    end if

    erlarg = erlarg - erlast

    if ( levcur+1 <= levmax ) then
      erlarg = erlarg + erro12
    end if
!
!  Test whether the interval to be bisected next is the
!  smallest interval.
!
    if ( .not. extrap ) then

      if ( level(maxerr)+1 <= levmax ) then
        cycle
      end if

      extrap = .true.
      nrmax = 2

    end if
!
!  The smallest interval has the largest error.
!  Before bisecting decrease the sum of the errors over the
!  larger intervals (erlarg) and perform extrapolation.
!
    if ( ierro /= 3 .and. erlarg > ertest ) then

      id = nrmax
      jupbnd = last
      if ( last > (2+limit/2) ) then
        jupbnd = limit+3-last
      end if

      do k = id, jupbnd
        maxerr = iord(nrmax)
        errmax = elist(maxerr)
        if ( level(maxerr)+1 <= levmax ) then
          go to 160
        end if
        nrmax = nrmax + 1
      end do

    end if
!
!  Perform extrapolation.
!
    numrl2 = numrl2 + 1
    rlist2(numrl2) = area

    if ( numrl2 <= 2 ) then
      go to 155
    end if

    call qextr ( numrl2, rlist2, reseps, abseps, res3la, nres )
    ktmin = ktmin+1

    if ( 5 < ktmin .and. abserr < 1.0e-03 * errsum ) then
      ier = 5
    end if

    if ( abseps < abserr ) then

      ktmin = 0
      abserr = abseps
      result = reseps
      correc = erlarg
      ertest = max ( epsabs, epsrel * abs(reseps) )

      if ( abserr < ertest ) then
        exit
      end if

    end if
!
!  Prepare bisection of the smallest interval.
!
    if ( numrl2 == 1 ) then
      noext = .true.
    end if

    if ( 5 <= ier ) then
      exit
    end if

155 continue

    maxerr = iord(1)
    errmax = elist(maxerr)
    nrmax = 1
    extrap = .false.
    levmax = levmax + 1
    erlarg = errsum

160 continue

  end do
!
!  Set the final result.
!
  if (exactly_equal(abserr, huge(abserr))) then
    go to 190
  end if

  if ( ( ier + ierro ) == 0 ) then
    go to 180
  end if

  if ( ierro == 3 ) then
    abserr = abserr + correc
  end if

  if ( ier == 0 ) then
    ier = 3
  end if

  if (is_not_zero(result) .and. is_not_zero(area)) then
    go to 175
  end if

  if ( errsum < abserr ) then
    go to 190
  end if

  if (is_zero(area)) then
    go to 210
  end if

  go to 180

175 continue

  if ( abserr / abs(result) > errsum / abs(area) ) then
    go to 190
  end if
!
!  Test on divergence.
!
  180 continue

  if ( ksgn == (-1) .and. max ( abs(result),abs(area)) <=  &
    resabs*1.0e-02 ) go to 210

  if ( 1.0e-02 > (result/area) .or. (result/area) > 1.0e+02 .or. &
    errsum > abs(area) ) then
    ier = 6
  end if

  go to 210
!
!  Compute global integral sum.
!
190 continue

  result = sum ( rlist(1:last) )

  abserr = errsum

210 continue

  if ( 2 < ier ) then
    ier = ier - 1
  end if

  result = result * sign

  return
end subroutine
subroutine qags ( f, a, b, epsabs, epsrel, result, abserr, neval, ier )

!*****************************************************************************80
!
!! QAGS estimates the integral of a function.
!
!  Discussion:
!
!    The routine calculates an approximation RESULT to a definite integral   
!      I = integral of F over (A,B),
!    hopefully satisfying
!      || I - RESULT || <= max ( EPSABS, EPSREL * ||I|| ).
!
!  Author:
!
!    Robert Piessens, Elise de Doncker-Kapenger, 
!    Christian Ueberhuber, David Kahaner
!
!  Reference:
!
!    Robert Piessens, Elise de Doncker-Kapenger, 
!    Christian Ueberhuber, David Kahaner,
!    QUADPACK, a Subroutine Package for Automatic Integration,
!    Springer Verlag, 1983
!
!  Parameters:
!
!    Input, external real F, the name of the function routine, of the form
!      function f ( x )
!      real f
!      real x
!    which evaluates the integrand function.
!
!    Input, real A, B, the limits of integration.
!
!    Input, real EPSABS, EPSREL, the absolute and relative accuracy requested.
!
!    Output, real RESULT, the estimated value of the integral.
!
!    Output, real ABSERR, an estimate of || I - RESULT ||.
!
!    Output, integer NEVAL, the number of times the integral was evaluated.
!
!    Output, integer IER, error flag.
!                     ier = 0 normal and reliable termination of the
!                             routine. it is assumed that the requested
!                             accuracy has been achieved.
!                     ier > 0 abnormal termination of the routine
!                             the estimates for integral and error are
!                             less reliable. it is assumed that the
!                             requested accuracy has not been achieved.
!                         = 1 maximum number of subdivisions allowed
!                             has been achieved. one can allow more sub-
!                             divisions by increasing the data value of
!                             limit in qags (and taking the according
!                             dimension adjustments into account).
!                             however, if this yields no improvement
!                             it is advised to analyze the integrand
!                             in order to determine the integration
!                             difficulties. if the position of a
!                             local difficulty can be determined (e.g.
!                             singularity, discontinuity within the
!                             interval) one will probably gain from
!                             splitting up the interval at this point
!                             and calling the integrator on the sub-
!                             ranges. if possible, an appropriate
!                             special-purpose integrator should be used,
!                             which is designed for handling the type
!                             of difficulty involved.
!                         = 2 the occurrence of roundoff error is detec-
!                             ted, which prevents the requested
!                             tolerance from being achieved.
!                             the error may be under-estimated.
!                         = 3 extremely bad integrand behavior occurs
!                             at some  points of the integration
!                             interval.
!                         = 4 the algorithm does not converge. roundoff
!                             error is detected in the extrapolation
!                             table. it is presumed that the requested
!                             tolerance cannot be achieved, and that the
!                             returned result is the best which can be
!                             obtained.
!                         = 5 the integral is probably divergent, or
!                             slowly convergent. it must be noted that
!                             divergence can occur with any other value
!                             of ier.
!                         = 6 the input is invalid, because
!                             epsabs < 0 and epsrel < 0,
!                             result, abserr and neval are set to zero.
!
!  Local Parameters:
!
!           alist     - list of left end points of all subintervals
!                       considered up to now
!           blist     - list of right end points of all subintervals
!                       considered up to now
!           rlist(i)  - approximation to the integral over
!                       (alist(i),blist(i))
!           rlist2    - array of dimension at least limexp+2 containing
!                       the part of the epsilon table which is still
!                       needed for further computations
!           elist(i)  - error estimate applying to rlist(i)
!           maxerr    - pointer to the interval with largest error
!                       estimate
!           errmax    - elist(maxerr)
!           erlast    - error on the interval currently subdivided
!                       (before that subdivision has taken place)
!           area      - sum of the integrals over the subintervals
!           errsum    - sum of the errors over the subintervals
!           errbnd    - requested accuracy max(epsabs,epsrel*
!                       abs(result))
!           *****1    - variable for the left interval
!           *****2    - variable for the right interval
!           last      - index for subdivision
!           nres      - number of calls to the extrapolation routine
!           numrl2    - number of elements currently in rlist2. if an
!                       appropriate approximation to the compounded
!                       integral has been obtained it is put in
!                       rlist2(numrl2) after numrl2 has been increased
!                       by one.
!           small     - length of the smallest interval considered
!                       up to now, multiplied by 1.5
!           erlarg    - sum of the errors over the intervals larger
!                       than the smallest interval considered up to now
!           extrap    - logical variable denoting that the routine is
!                       attempting to perform extrapolation i.e. before
!                       subdividing the smallest interval we try to
!                       decrease the value of erlarg.
!           noext     - logical variable denoting that extrapolation
!                       is no longer allowed (true value)
!
  implicit none

  integer, parameter :: limit = 500

  real a
  real abseps
  real abserr
  real alist(limit)
  real area
  real area1
  real area12
  real area2
  real a1
  real a2
  real b
  real blist(limit)
  real b1
  real b2
  real correc
  real defabs
  real defab1
  real defab2
  real dres
  real elist(limit)
  real epsabs
  real epsrel
  real erlarg
  real erlast
  real errbnd
  real errmax
  real error1
  real error2
  real erro12
  real errsum
  real ertest
  logical extrap
  procedure(scalar_func) :: f
  integer id
  integer ier
  integer ierro
  integer iord(limit)
  integer iroff1
  integer iroff2
  integer iroff3
  integer jupbnd
  integer k
  integer ksgn
  integer ktmin
  integer last
  logical noext
  integer maxerr
  integer neval
  integer nres
  integer nrmax
  integer numrl2
  real resabs
  real reseps
  real result
  real res3la(3)
  real rlist(limit)
  real rlist2(52)
  real small
!
!  The dimension of rlist2 is determined by the value of
!  limexp in QEXTR (rlist2 should be of dimension
!  (limexp+2) at least).
!
!  Test on validity of parameters.
!
  ier = 0
  neval = 0
  last = 0
  result = 0.0e+00
  abserr = 0.0e+00
  alist(1) = a
  blist(1) = b
  rlist(1) = 0.0e+00
  elist(1) = 0.0e+00

  if ( epsabs < 0.0e+00 .and. epsrel < 0.0e+00 ) then
    ier = 6
    return
  end if
!
!  First approximation to the integral.
!
  ierro = 0
  call qk21 ( f, a, b, result, abserr, defabs, resabs )
!
!  Test on accuracy.
!
  dres = abs ( result )
  errbnd = max ( epsabs, epsrel * dres )
  last = 1
  rlist(1) = result
  elist(1) = abserr
  iord(1) = 1

  if ( abserr <= 1.0e+02 * epsilon ( defabs ) * defabs .and. &
    abserr > errbnd ) then
    ier = 2
  end if

  if ( limit == 1 ) then
    ier = 1
  end if

  if ( ier /= 0 .or. (abserr <= errbnd .and. not_exactly_equal(abserr, resabs)) &
       .or. is_zero(abserr)) go to 140
!
!  Initialization.
!
  rlist2(1) = result
  errmax = abserr
  maxerr = 1
  area = result
  errsum = abserr
  abserr = huge ( abserr )
  nrmax = 1
  nres = 0
  numrl2 = 2
  ktmin = 0
  extrap = .false.
  noext = .false.
  iroff1 = 0
  iroff2 = 0
  iroff3 = 0

  if ( dres >= (1.0e+00-5.0e+01* epsilon ( defabs ) ) * defabs ) then
    ksgn = 1
  else
    ksgn = -1
  end if

  do last = 2, limit
!
!  Bisect the subinterval with the nrmax-th largest error estimate.
!
    a1 = alist(maxerr)
    b1 = 5.0e-01 * ( alist(maxerr) + blist(maxerr) )
    a2 = b1
    b2 = blist(maxerr)
    erlast = errmax
    call qk21 ( f, a1, b1, area1, error1, resabs, defab1 )
    call qk21 ( f, a2, b2, area2, error2, resabs, defab2 )
!
!  Improve previous approximations to integral and error
!  and test for accuracy.
!
    area12 = area1+area2
    erro12 = error1+error2
    errsum = errsum+erro12-errmax
    area = area+area12-rlist(maxerr)

    if (exactly_equal(defab1, error1) .or. exactly_equal(defab2, error2)) go to 15

    if ( abs ( rlist(maxerr) - area12) > 1.0e-05 * abs(area12) &
      .or. erro12 < 9.9e-01 * errmax ) go to 10

    if ( extrap ) then
      iroff2 = iroff2+1
    else
      iroff1 = iroff1+1
    end if

10  continue

    if ( last > 10 .and. erro12 > errmax ) then
      iroff3 = iroff3+1
    end if

15  continue

    rlist(maxerr) = area1
    rlist(last) = area2
    errbnd = max ( epsabs, epsrel*abs(area) )
!
!  Test for roundoff error and eventually set error flag.
!
    if ( iroff1+iroff2 >= 10 .or. iroff3 >= 20 ) then
      ier = 2
    end if

    if ( iroff2 >= 5 ) then
      ierro = 3
    end if
!
!  Set error flag in the case that the number of subintervals
!  equals limit.
!
    if ( last == limit ) then
      ier = 1
    end if
!
!  Set error flag in the case of bad integrand behavior
!  at a point of the integration range.
!
    if ( max ( abs(a1),abs(b2)) <= (1.0e+00+1.0e+03* epsilon ( a1 ) )* &
      (abs(a2)+1.0e+03* tiny ( a2 ) ) ) then
      ier = 4
    end if
!
!  Append the newly-created intervals to the list.
!
    if ( error2 <= error1 ) then
      alist(last) = a2
      blist(maxerr) = b1
      blist(last) = b2
      elist(maxerr) = error1
      elist(last) = error2
    else
      alist(maxerr) = a2
      alist(last) = a1
      blist(last) = b1
      rlist(maxerr) = area2
      rlist(last) = area1
      elist(maxerr) = error2
      elist(last) = error1
    end if
!
!  Call QSORT to maintain the descending ordering
!  in the list of error estimates and select the subinterval
!  with nrmax-th largest error estimate (to be bisected next).
!
    call qsort ( limit, last, maxerr, errmax, elist, iord, nrmax )

    if ( errsum <= errbnd ) go to 115

    if ( ier /= 0 ) then
      exit
    end if

    if ( last == 2 ) go to 80
    if ( noext ) go to 90

    erlarg = erlarg-erlast

    if ( abs(b1-a1) > small ) then
      erlarg = erlarg+erro12
    end if
!
!  Test whether the interval to be bisected next is the
!  smallest interval.
!
    if ( .not. extrap ) then
      if ( abs(blist(maxerr)-alist(maxerr)) > small ) go to 90
      extrap = .true.
      nrmax = 2
    end if

!40  continue
!
!  The smallest interval has the largest error.
!  Before bisecting decrease the sum of the errors over the
!  larger intervals (erlarg) and perform extrapolation.
!
    if ( ierro /= 3 .and. erlarg > ertest ) then

      id = nrmax
      jupbnd = last

      if ( last > (2+limit/2) ) then
        jupbnd = limit+3-last
      end if

      do k = id, jupbnd
        maxerr = iord(nrmax)
        errmax = elist(maxerr)
        if ( abs(blist(maxerr)-alist(maxerr)) > small ) then
          go to 90
        end if
        nrmax = nrmax+1
      end do

    end if
!
!  Perform extrapolation.
!
!60  continue

    numrl2 = numrl2+1
    rlist2(numrl2) = area
    call qextr ( numrl2, rlist2, reseps, abseps, res3la, nres )
    ktmin = ktmin+1

    if ( ktmin > 5 .and. abserr < 1.0e-03 * errsum ) then
      ier = 5
    end if

    if ( abseps < abserr ) then

      ktmin = 0
      abserr = abseps
      result = reseps
      correc = erlarg
      ertest = max ( epsabs,epsrel*abs(reseps))

      if ( abserr <= ertest ) then
        exit
      end if

    end if
!
!  Prepare bisection of the smallest interval.
!
    if ( numrl2 == 1 ) then
      noext = .true.
    end if

    if ( ier == 5 ) then
      exit
    end if

    maxerr = iord(1)
    errmax = elist(maxerr)
    nrmax = 1
    extrap = .false.
    small = small * 5.0e-01
    erlarg = errsum
    go to 90

80  continue

    small = abs ( b - a ) * 3.75e-01
    erlarg = errsum
    ertest = errbnd
    rlist2(2) = area

90  continue

  end do
!
!  Set final result and error estimate.
!
  if (exactly_equal(abserr, huge(abserr))) then
    go to 115
  end if

  if ( ier + ierro == 0 ) then
    go to 110
  end if

  if ( ierro == 3 ) then
    abserr = abserr + correc
  end if

  if ( ier == 0 ) then
    ier = 3
  end if

  if (is_not_zero(result) .and. is_not_zero(area)) then
    go to 105
  end if

  if ( abserr > errsum ) go to 115
  if (is_zero(area)) go to 130
  go to 110

105 continue

  if ( abserr/abs(result) > errsum/abs(area) ) go to 115
!
!  Test on divergence.
!
110 continue

  if ( ksgn == (-1).and.max ( abs(result),abs(area)) <=  &
   defabs*1.0e-02 ) go to 130

  if ( 1.0e-02 > (result/area) .or. (result/area) > 1.0e+02 &
   .or. errsum > abs(area) ) then
    ier = 6
  end if

  go to 130
!
!  Compute global integral sum.
!
115 continue

  result = sum ( rlist(1:last) )

  abserr = errsum

130 continue
 
  if ( 2 < ier ) then
    ier = ier - 1
  end if

140 continue

  neval = 42*last-21

  return
end subroutine
subroutine qawc ( f, a, b, c, epsabs, epsrel, result, abserr, neval, ier )

!*****************************************************************************80
!
!! QAWC computes a Cauchy principal value.
!
!  Discussion:
!
!    The routine calculates an approximation RESULT to a Cauchy principal
!    value 
!      I = integral of F*W over (A,B),
!    with
!      W(X) = 1 / (X-C),
!    with C distinct from A and B, hopefully satisfying
!      || I - RESULT || <= max ( EPSABS, EPSREL * ||I|| ).
!
!  Author:
!
!    Robert Piessens, Elise de Doncker-Kapenger, 
!    Christian Ueberhuber, David Kahaner
!
!  Reference:
!
!    Robert Piessens, Elise de Doncker-Kapenger, 
!    Christian Ueberhuber, David Kahaner,
!    QUADPACK, a Subroutine Package for Automatic Integration,
!    Springer Verlag, 1983
!
!  Parameters:
!
!    Input, external real F, the name of the function routine, of the form
!      function f ( x )
!      real f
!      real x
!    which evaluates the integrand function.
!
!    Input, real A, B, the limits of integration.
!
!    Input, real C, a parameter in the weight function, which must
!    not be equal to A or B.
!
!    Input, real EPSABS, EPSREL, the absolute and relative accuracy requested.
!
!    Output, real RESULT, the estimated value of the integral.
!
!    Output, real ABSERR, an estimate of || I - RESULT ||.
!
!    Output, integer NEVAL, the number of times the integral was evaluated.
!
!            ier    - integer
!                     ier = 0 normal and reliable termination of the
!                             routine. it is assumed that the requested
!                             accuracy has been achieved.
!                     ier > 0 abnormal termination of the routine
!                             the estimates for integral and error are
!                             less reliable. it is assumed that the
!                             requested accuracy has not been achieved.
!                     ier = 1 maximum number of subdivisions allowed
!                             has been achieved. one can allow more sub-
!                             divisions by increasing the data value of
!                             limit in qawc (and taking the according
!                             dimension adjustments into account).
!                             however, if this yields no improvement it
!                             is advised to analyze the integrand in
!                             order to determine the integration
!                             difficulties. if the position of a local
!                             difficulty can be determined (e.g.
!                             singularity, discontinuity within the
!                             interval one will probably gain from
!                             splitting up the interval at this point
!                             and calling appropriate integrators on the
!                             subranges.
!                         = 2 the occurrence of roundoff error is detec-
!                             ted, which prevents the requested
!                             tolerance from being achieved.
!                         = 3 extremely bad integrand behavior occurs
!                             at some points of the integration
!                             interval.
!                         = 6 the input is invalid, because
!                             c = a or c = b or
!                             epsabs < 0 and epsrel < 0,
!                             result, abserr, neval are set to zero.
!
!  Local parameters:
!
!    LIMIT is the maximum number of subintervals allowed in the
!    subdivision process of qawce. take care that limit >= 1.
!
  implicit none

  integer, parameter :: limit = 500

  real a
  real abserr
  real alist(limit)
  real b
  real blist(limit)
  real elist(limit)
  real c
  real epsabs
  real epsrel
  procedure(scalar_func) :: f
  integer ier
  integer iord(limit)
  integer last
  integer neval
  real result
  real rlist(limit)

  call qawce ( f, a, b, c, epsabs, epsrel, limit, result, abserr, neval, ier, &
    alist, blist, rlist, elist, iord, last )

  return
end subroutine
subroutine qawce ( f, a, b, c, epsabs, epsrel, limit, result, abserr, neval, &
  ier, alist, blist, rlist, elist, iord, last )

!*****************************************************************************80
!
!! QAWCE computes a Cauchy principal value.
!
!  Discussion:
!
!    The routine calculates an approximation RESULT to a Cauchy principal
!    value   
!      I = integral of F*W over (A,B),
!    with
!      W(X) = 1 / ( X - C ),
!    with C distinct from A and B, hopefully satisfying
!      | I - RESULT | <= max ( EPSABS, EPSREL * |I| ).
!
!  Author:
!
!    Robert Piessens, Elise de Doncker-Kapenger, 
!    Christian Ueberhuber, David Kahaner
!
!  Reference:
!
!    Robert Piessens, Elise de Doncker-Kapenger, 
!    Christian Ueberhuber, David Kahaner,
!    QUADPACK, a Subroutine Package for Automatic Integration,
!    Springer Verlag, 1983
!
!  Parameters:
!
!    Input, external real F, the name of the function routine, of the form
!      function f ( x )
!      real f
!      real x
!    which evaluates the integrand function.
!
!    Input, real A, B, the limits of integration.
!
!    Input, real C, a parameter in the weight function, which cannot be
!    equal to A or B.
!
!    Input, real EPSABS, EPSREL, the absolute and relative accuracy requested.
!
!    Input, integer LIMIT, the upper bound on the number of subintervals that
!    will be used in the partition of [A,B].  LIMIT is typically 500.
!
!    Output, real RESULT, the estimated value of the integral.
!
!    Output, real ABSERR, an estimate of || I - RESULT ||.
!
!    Output, integer NEVAL, the number of times the integral was evaluated.
!
!            ier    - integer
!                     ier = 0 normal and reliable termination of the
!                             routine. it is assumed that the requested
!                             accuracy has been achieved.
!                     ier > 0 abnormal termination of the routine
!                             the estimates for integral and error are
!                             less reliable. it is assumed that the
!                             requested accuracy has not been achieved.
!                     ier = 1 maximum number of subdivisions allowed
!                             has been achieved. one can allow more sub-
!                             divisions by increasing the value of
!                             limit. however, if this yields no
!                             improvement it is advised to analyze the
!                             integrand, in order to determine the
!                             integration difficulties.  if the position
!                             of a local difficulty can be determined
!                             (e.g. singularity, discontinuity within
!                             the interval) one will probably gain
!                             from splitting up the interval at this
!                             point and calling appropriate integrators
!                             on the subranges.
!                         = 2 the occurrence of roundoff error is detec-
!                             ted, which prevents the requested
!                             tolerance from being achieved.
!                         = 3 extremely bad integrand behavior occurs
!                             at some interior points of the integration
!                             interval.
!                         = 6 the input is invalid, because
!                             c = a or c = b or
!                             epsabs < 0 and epsrel < 0,
!                             or limit < 1.
!                             result, abserr, neval, rlist(1), elist(1),
!                             iord(1) and last are set to zero.
!                             alist(1) and blist(1) are set to a and b
!                             respectively.
!
!    Workspace, real ALIST(LIMIT), BLIST(LIMIT), contains in entries 1 
!    through LAST the left and right ends of the partition subintervals.
!
!    Workspace, real RLIST(LIMIT), contains in entries 1 through LAST
!    the integral approximations on the subintervals.
!
!    Workspace, real ELIST(LIMIT), contains in entries 1 through LAST
!    the absolute error estimates on the subintervals.
!
!            iord    - integer
!                      vector of dimension at least limit, the first k
!                      elements of which are pointers to the error
!                      estimates over the subintervals, so that
!                      elist(iord(1)), ...,  elist(iord(k)) with
!                      k = last if last <= (limit/2+2), and
!                      k = limit+1-last otherwise, form a decreasing
!                      sequence.
!
!            last    - integer
!                      number of subintervals actually produced in
!                      the subdivision process
!
!  Local parameters:
!
!           alist     - list of left end points of all subintervals
!                       considered up to now
!           blist     - list of right end points of all subintervals
!                       considered up to now
!           rlist(i)  - approximation to the integral over
!                       (alist(i),blist(i))
!           elist(i)  - error estimate applying to rlist(i)
!           maxerr    - pointer to the interval with largest error
!                       estimate
!           errmax    - elist(maxerr)
!           area      - sum of the integrals over the subintervals
!           errsum    - sum of the errors over the subintervals
!           errbnd    - requested accuracy max(epsabs,epsrel*
!                       abs(result))
!           *****1    - variable for the left subinterval
!           *****2    - variable for the right subinterval
!           last      - index for subdivision
!
  implicit none

  integer limit

  real a
  real aa
  real abserr
  real alist(limit)
  real area
  real area1
  real area12
  real area2
  real a1
  real a2
  real b
  real bb
  real blist(limit)
  real b1
  real b2
  real c
  real elist(limit)
  real epsabs
  real epsrel
  real errbnd
  real errmax
  real error1
  real error2
  real erro12
  real errsum
  procedure(scalar_func) :: f
  integer ier
  integer iord(limit)
  integer iroff1
  integer iroff2
  integer krule
  integer last
  integer maxerr
  integer nev
  integer neval
  integer nrmax
  real result
  real rlist(limit)
!
!  Test on validity of parameters.
!
  ier = 0
  neval = 0
  last = 0
  alist(1) = a
  blist(1) = b
  rlist(1) = 0.0e+00
  elist(1) = 0.0e+00
  iord(1) = 0
  result = 0.0e+00
  abserr = 0.0e+00

  if (exactly_equal(c, a)) then
    ier = 6
    return
  else if (exactly_equal(c, b)) then
    ier = 6
    return
  else if ( epsabs < 0.0e+00 .and. epsrel < 0.0e+00 ) then
    ier = 6
    return
  end if
!
!  First approximation to the integral.
!
  if ( a <= b ) then
    aa = a
    bb = b
  else
    aa = b
    bb = a
  end if

  krule = 1
  call qc25c ( f, aa, bb, c, result, abserr, krule, neval )
  last = 1
  rlist(1) = result
  elist(1) = abserr
  iord(1) = 1
  alist(1) = a
  blist(1) = b
!
!  Test on accuracy.
!
  errbnd = max ( epsabs, epsrel * abs(result) )

  if ( limit == 1 ) then
    ier = 1
    go to 70
  end if

  if ( abserr < min ( 1.0e-02 * abs(result), errbnd)  ) then
    go to 70
  end if
!
!  Initialization
!
  alist(1) = aa
  blist(1) = bb
  rlist(1) = result
  errmax = abserr
  maxerr = 1
  area = result
  errsum = abserr
  nrmax = 1
  iroff1 = 0
  iroff2 = 0

  do last = 2, limit
!
!  Bisect the subinterval with nrmax-th largest error estimate.
!
    a1 = alist(maxerr)
    b1 = 5.0e-01*(alist(maxerr)+blist(maxerr))
    b2 = blist(maxerr)

    if ( c <= b1 .and. a1 < c ) then
      b1 = 5.0e-01*(c+b2)
    end if

    if ( b1 < c .and. c < b2 ) then
      b1 = 5.0e-01 * ( a1 + c )
    end if

    a2 = b1
    krule = 2

    call qc25c ( f, a1, b1, c, area1, error1, krule, nev )
    neval = neval+nev

    call qc25c ( f, a2, b2, c, area2, error2, krule, nev )
    neval = neval+nev
!
!  Improve previous approximations to integral and error
!  and test for accuracy.
!
    area12 = area1 + area2
    erro12 = error1 + error2
    errsum = errsum + erro12 - errmax
    area = area + area12 - rlist(maxerr)

    if ( abs ( rlist(maxerr)-area12) < 1.0e-05 * abs(area12) &
      .and. erro12 >= 9.9e-01 * errmax .and. krule == 0 ) &
      iroff1 = iroff1+1

    if ( last > 10.and.erro12 > errmax .and. krule == 0 ) then
      iroff2 = iroff2+1
    end if

    rlist(maxerr) = area1
    rlist(last) = area2
    errbnd = max ( epsabs, epsrel * abs(area) )

    if ( errsum > errbnd ) then
!
!  Test for roundoff error and eventually set error flag.
!
      if ( iroff1 >= 6 .and. iroff2 > 20 ) then
        ier = 2
      end if
!
!  Set error flag in the case that number of interval
!  bisections exceeds limit.
!
      if ( last == limit ) then
        ier = 1
      end if
!
!  Set error flag in the case of bad integrand behavior at
!  a point of the integration range.
!
      if ( max ( abs(a1), abs(b2) ) <= ( 1.0e+00 + 1.0e+03 * epsilon ( a1 ) ) &
        *( abs(a2)+1.0e+03* tiny ( a2 ) )) then
        ier = 3
      end if

    end if
!
!  Append the newly-created intervals to the list.
!
    if ( error2 <= error1 ) then
      alist(last) = a2
      blist(maxerr) = b1
      blist(last) = b2
      elist(maxerr) = error1
      elist(last) = error2
    else
      alist(maxerr) = a2
      alist(last) = a1
      blist(last) = b1
      rlist(maxerr) = area2
      rlist(last) = area1
      elist(maxerr) = error2
      elist(last) = error1
    end if
!
!  Call QSORT to maintain the descending ordering
!  in the list of error estimates and select the subinterval
!  with NRMAX-th largest error estimate (to be bisected next).
!
    call qsort ( limit, last, maxerr, errmax, elist, iord, nrmax )

    if ( ier /= 0 .or. errsum <= errbnd ) then
      exit
    end if

  end do
!
!  Compute final result.
!
  result = sum ( rlist(1:last) )

  abserr = errsum

70 continue 

  if (exactly_equal(aa, b)) then
    result = - result
  end if

  return
end subroutine
subroutine qawf ( f, a, omega, integr, epsabs, result, abserr, neval, ier )

!*****************************************************************************80
!
!! QAWF computes Fourier integrals over the interval [ A, +Infinity ).
!
!  Discussion:
!
!    The routine calculates an approximation RESULT to a definite integral  
! 
!      I = integral of F*COS(OMEGA*X) 
!    or 
!      I = integral of F*SIN(OMEGA*X) 
!
!    over the interval [A,+Infinity), hopefully satisfying
!
!      || I - RESULT || <= EPSABS.
!
!    If OMEGA = 0 and INTEGR = 1, the integral is calculated by means 
!    of QAGI, and IER has the meaning as described in the comments of QAGI.
!
!  Author:
!
!    Robert Piessens, Elise de Doncker-Kapenger, 
!    Christian Ueberhuber, David Kahaner
!
!  Reference:
!
!    Robert Piessens, Elise de Doncker-Kapenger, 
!    Christian Ueberhuber, David Kahaner,
!    QUADPACK, a Subroutine Package for Automatic Integration,
!    Springer Verlag, 1983
!
!  Parameters:
!
!    Input, external real F, the name of the function routine, of the form
!      function f ( x )
!      real f
!      real x
!    which evaluates the integrand function.
!
!    Input, real A, the lower limit of integration.
!
!    Input, real OMEGA, the parameter in the weight function.
!
!    Input, integer INTEGR, indicates which weight functions is used
!    = 1, w(x) = cos(omega*x)
!    = 2, w(x) = sin(omega*x)
!
!    Input, real EPSABS, the absolute accuracy requested.
!
!    Output, real RESULT, the estimated value of the integral.
!
!    Output, real ABSERR, an estimate of || I - RESULT ||.
!
!    Output, integer NEVAL, the number of times the integral was evaluated.
!
!            ier    - integer
!                     ier = 0 normal and reliable termination of the
!                             routine. it is assumed that the
!                             requested accuracy has been achieved.
!                     ier > 0 abnormal termination of the routine.
!                             the estimates for integral and error are
!                             less reliable. it is assumed that the
!                             requested accuracy has not been achieved.
!                    if omega /= 0
!                     ier = 6 the input is invalid because
!                             (integr /= 1 and integr /= 2) or
!                              epsabs <= 0
!                              result, abserr, neval, lst are set to
!                              zero.
!                         = 7 abnormal termination of the computation
!                             of one or more subintegrals
!                         = 8 maximum number of cycles allowed
!                             has been achieved, i.e. of subintervals
!                             (a+(k-1)c,a+kc) where
!                             c = (2*int(abs(omega))+1)*pi/abs(omega),
!                             for k = 1, 2, ...
!                         = 9 the extrapolation table constructed for
!                             convergence acceleration of the series
!                             formed by the integral contributions
!                             over the cycles, does not converge to
!                             within the requested accuracy.
!
!  Local parameters:
!
!    Integer LIMLST, gives an upper bound on the number of cycles, LIMLST >= 3.
!    if limlst < 3, the routine will end with ier = 6.
!
!    Integer MAXP1, an upper bound on the number of Chebyshev moments which 
!    can be stored, i.e. for the intervals of lengths abs(b-a)*2**(-l), 
!    l = 0,1, ..., maxp1-2, maxp1 >= 1.  if maxp1 < 1, the routine will end
!    with ier = 6.
!
  implicit none

  integer, parameter :: limit = 500
  integer, parameter :: limlst = 50
  integer, parameter :: maxp1 = 21

  real a
  real abserr
  real alist(limit)
  real blist(limit)
  real chebmo(maxp1,25)
  real elist(limit)
  real epsabs
  real erlst(limlst)
  procedure(scalar_func) :: f
  integer ier
  integer integr
  integer iord(limit)
  integer ierlst(limlst)
  integer lst
  integer neval
  integer nnlog(limit)
  real omega
  real result
  real rlist(limit)
  real rslst(limlst)

  ier = 6
  neval = 0
  result = 0.0e+00
  abserr = 0.0e+00

  if ( limlst < 3 .or. maxp1 < 1 ) then
    return
  end if

  call qawfe ( f, a, omega, integr, epsabs, limlst, limit, maxp1, &
    result, abserr, neval, ier, rslst, erlst, ierlst, lst, alist, blist, &
    rlist, elist, iord, nnlog, chebmo )

  return
end subroutine
subroutine qawfe ( f, a, omega, integr, epsabs, limlst, limit, maxp1, &
  result, abserr, neval, ier, rslst, erlst, ierlst, lst, alist, blist, &
  rlist, elist, iord, nnlog, chebmo )

!*****************************************************************************80
!
!! QAWFE computes Fourier integrals.
!
!  Discussion:
!
!    The routine calculates an approximation RESULT to a definite integral   
!      I = integral of F*COS(OMEGA*X) or F*SIN(OMEGA*X) over (A,+Infinity),
!    hopefully satisfying
!      || I - RESULT || <= EPSABS.
!
!  Author:
!
!    Robert Piessens, Elise de Doncker-Kapenger, 
!    Christian Ueberhuber, David Kahaner
!
!  Reference:
!
!    Robert Piessens, Elise de Doncker-Kapenger, 
!    Christian Ueberhuber, David Kahaner,
!    QUADPACK, a Subroutine Package for Automatic Integration,
!    Springer Verlag, 1983
!
!  Parameters:
!
!    Input, external real F, the name of the function routine, of the form
!      function f ( x )
!      real f
!      real x
!    which evaluates the integrand function.
!
!    Input, real A, the lower limit of integration.
!
!    Input, real OMEGA, the parameter in the weight function.
!
!    Input, integer INTEGR, indicates which weight function is used
!    = 1      w(x) = cos(omega*x)
!    = 2      w(x) = sin(omega*x)
!
!    Input, real EPSABS, the absolute accuracy requested.
!
!    Input, integer LIMLST, an upper bound on the number of cycles.
!    LIMLST must be at least 1.  In fact, if LIMLST < 3, the routine 
!    will end with IER= 6.
!
!    Input, integer LIMIT, an upper bound on the number of subintervals 
!    allowed in the partition of each cycle, limit >= 1.
!
!            maxp1  - integer
!                     gives an upper bound on the number of
!                     Chebyshev moments which can be stored, i.e.
!                     for the intervals of lengths abs(b-a)*2**(-l),
!                     l=0,1, ..., maxp1-2, maxp1 >= 1
!
!    Output, real RESULT, the estimated value of the integral.
!
!    Output, real ABSERR, an estimate of || I - RESULT ||.
!
!    Output, integer NEVAL, the number of times the integral was evaluated.
!
!            ier    - ier = 0 normal and reliable termination of
!                             the routine. it is assumed that the
!                             requested accuracy has been achieved.
!                     ier > 0 abnormal termination of the routine
!                             the estimates for integral and error
!                             are less reliable. it is assumed that
!                             the requested accuracy has not been
!                             achieved.
!                    if omega /= 0
!                     ier = 6 the input is invalid because
!                             (integr /= 1 and integr /= 2) or
!                              epsabs <= 0 or limlst < 3.
!                              result, abserr, neval, lst are set
!                              to zero.
!                         = 7 bad integrand behavior occurs within
!                             one or more of the cycles. location
!                             and type of the difficulty involved
!                             can be determined from the vector ierlst.
!                             here lst is the number of cycles actually
!                             needed (see below).
!                             ierlst(k) = 1 the maximum number of
!                                           subdivisions (= limit)
!                                           has been achieved on the
!                                           k th cycle.
!                                       = 2 occurence of roundoff
!                                           error is detected and
!                                           prevents the tolerance
!                                           imposed on the k th cycle
!                                           from being acheived.
!                                       = 3 extremely bad integrand
!                                           behavior occurs at some
!                                           points of the k th cycle.
!                                       = 4 the integration procedure
!                                           over the k th cycle does
!                                           not converge (to within the
!                                           required accuracy) due to
!                                           roundoff in the
!                                           extrapolation procedure
!                                           invoked on this cycle. it
!                                           is assumed that the result
!                                           on this interval is the
!                                           best which can be obtained.
!                                       = 5 the integral over the k th
!                                           cycle is probably divergent
!                                           or slowly convergent. it
!                                           must be noted that
!                                           divergence can occur with
!                                           any other value of
!                                           ierlst(k).
!                         = 8 maximum number of  cycles  allowed
!                             has been achieved, i.e. of subintervals
!                             (a+(k-1)c,a+kc) where
!                             c = (2*int(abs(omega))+1)*pi/abs(omega),
!                             for k = 1, 2, ..., lst.
!                             one can allow more cycles by increasing
!                             the value of limlst (and taking the
!                             according dimension adjustments into
!                             account).
!                             examine the array iwork which contains
!                             the error flags over the cycles, in order
!                             to eventual look for local integration
!                             difficulties.
!                             if the position of a local difficulty can
!                             be determined (e.g. singularity,
!                             discontinuity within the interval)
!                             one will probably gain from splitting
!                             up the interval at this point and
!                             calling appopriate integrators on the
!                             subranges.
!                         = 9 the extrapolation table constructed for
!                             convergence acceleration of the series
!                             formed by the integral contributions
!                             over the cycles, does not converge to
!                             within the required accuracy.
!                             as in the case of ier = 8, it is advised
!                             to examine the array iwork which contains
!                             the error flags on the cycles.
!                    if omega = 0 and integr = 1,
!                    the integral is calculated by means of qagi
!                    and ier = ierlst(1) (with meaning as described
!                    for ierlst(k), k = 1).
!
!            rslst  - real
!                     vector of dimension at least limlst
!                     rslst(k) contains the integral contribution
!                     over the interval (a+(k-1)c,a+kc) where
!                     c = (2*int(abs(omega))+1)*pi/abs(omega),
!                     k = 1, 2, ..., lst.
!                     note that, if omega = 0, rslst(1) contains
!                     the value of the integral over (a,infinity).
!
!            erlst  - real
!                     vector of dimension at least limlst
!                     erlst(k) contains the error estimate
!                     corresponding with rslst(k).
!
!            ierlst - integer
!                     vector of dimension at least limlst
!                     ierlst(k) contains the error flag corresponding
!                     with rslst(k). for the meaning of the local error
!                     flags see description of output parameter ier.
!
!            lst    - integer
!                     number of subintervals needed for the integration
!                     if omega = 0 then lst is set to 1.
!
!            alist, blist, rlist, elist - real
!                     vector of dimension at least limit,
!
!            iord, nnlog - integer
!                     vector of dimension at least limit, providing
!                     space for the quantities needed in the
!                     subdivision process of each cycle
!
!            chebmo - real
!                     array of dimension at least (maxp1,25),
!                     providing space for the Chebyshev moments
!                     needed within the cycles
!
!  Local parameters:
!
!           c1, c2    - end points of subinterval (of length
!                       cycle)
!           cycle     - (2*int(abs(omega))+1)*pi/abs(omega)
!           psum      - vector of dimension at least (limexp+2)
!                       (see routine qextr)
!                       psum contains the part of the epsilon table
!                       which is still needed for further computations.
!                       each element of psum is a partial sum of
!                       the series which should sum to the value of
!                       the integral.
!           errsum    - sum of error estimates over the
!                       subintervals, calculated cumulatively
!           epsa      - absolute tolerance requested over current
!                       subinterval
!           chebmo    - array containing the modified Chebyshev
!                       moments (see also routine qc25o)
!
  implicit none

  integer limit
  integer limlst
  integer maxp1

  real a
  real abseps
  real abserr
  real alist(limit)
  real blist(limit)
  real chebmo(maxp1,25)
  real correc
  real cycle
  real c1
  real c2
  real dl
! real dla
  real drl
  real elist(limit)
  real ep
  real eps
  real epsa
  real epsabs
  real erlst(limlst)
  real errsum
  procedure(scalar_func) :: f
  real fact
  integer ier
  integer ierlst(limlst)
  integer integr
  integer iord(limit)
  integer ktmin
  integer l
  integer ll
  integer lst
  integer momcom
  integer nev
  integer neval
  integer nnlog(limit)
  integer nres
  integer numrl2
  real omega
  real, parameter :: p = 0.9E+00
  real, parameter :: pi = 3.1415926535897932E+00
  real p1
  real psum(52)
  real reseps
  real result
  real res3la(3)
  real rlist(limit)
  real rslst(limlst)
!
!  The dimension of  psum  is determined by the value of
!  limexp in QEXTR (psum must be
!  of dimension (limexp+2) at least).
!
!  Test on validity of parameters.
!
  result = 0.0e+00
  abserr = 0.0e+00
  neval = 0
  lst = 0
  ier = 0

  if ( (integr /= 1 .and. integr /= 2 ) .or. &
    epsabs <= 0.0e+00 .or. &
    limlst < 3 ) then
    ier = 6
    return
  end if

  if (is_zero(omega)) then

    if ( integr == 1 ) then
      call qagi ( f, 0.0e+00, 1, epsabs, 0.0e+00, result, abserr, neval, ier )
    else
      result = 0.0E+00
      abserr = 0.0E+00
      neval = 0
      ier = 0
    end if

    rslst(1) = result
    erlst(1) = abserr
    ierlst(1) = ier
    lst = 1

    return
  end if
!
!  Initializations.
!
  l = int ( abs ( omega ) )
  dl = 2 * l + 1
  cycle = dl * pi / abs ( omega )
  ier = 0
  ktmin = 0
  neval = 0
  numrl2 = 0
  nres = 0
  c1 = a
  c2 = cycle+a
  p1 = 1.0e+00-p
  eps = epsabs

  if ( epsabs > tiny ( epsabs ) / p1 ) then
    eps = epsabs * p1
  end if

  ep = eps
  fact = 1.0e+00
  correc = 0.0e+00
  abserr = 0.0e+00
  errsum = 0.0e+00

  do lst = 1, limlst
!
!  Integrate over current subinterval.
!
!   dla = lst
    epsa = eps * fact

    call qfour ( f, c1, c2, omega, integr, epsa, 0.0e+00, limit, lst, maxp1, &
      rslst(lst), erlst(lst), nev, ierlst(lst), alist, blist, rlist, elist, &
      iord, nnlog, momcom, chebmo )

    neval = neval + nev
    fact = fact * p
    errsum = errsum + erlst(lst)
    drl = 5.0e+01 * abs(rslst(lst))
!
!  Test on accuracy with partial sum.
!
    if ((errsum+drl) <= epsabs.and.lst >= 6) then
      go to 80
    end if

    correc = max ( correc,erlst(lst))

    if ( ierlst(lst) /= 0 ) then
      eps = max ( ep,correc*p1)
      ier = 7
    end if

    if ( ier == 7 .and. (errsum+drl) <= correc*1.0e+01.and. lst > 5) go to 80

    numrl2 = numrl2+1

    if ( lst <= 1 ) then
      psum(1) = rslst(1)
      go to 40
    end if

    psum(numrl2) = psum(ll) + rslst(lst)

    if ( lst == 2 ) then
      go to 40
    end if
!
!  Test on maximum number of subintervals
!
    if ( lst == limlst ) then
      ier = 8
    end if
!
!  Perform new extrapolation
!
    call qextr ( numrl2, psum, reseps, abseps, res3la, nres )
!
!  Test whether extrapolated result is influenced by roundoff
!
    ktmin = ktmin + 1

    if ( ktmin >= 15 .and. abserr <= 1.0e-03 * (errsum+drl) ) then
      ier = 9
    end  if

    if ( abseps <= abserr .or. lst == 3 ) then

      abserr = abseps
      result = reseps
      ktmin = 0
!
!  If IER is not 0, check whether direct result (partial
!  sum) or extrapolated result yields the best integral
!  approximation
!
      if ( ( abserr + 1.0e+01 * correc ) <= epsabs ) then
        exit
      end if

      if ( abserr <= epsabs .and. 1.0e+01 * correc >= epsabs ) then
        exit
      end if

    end if

    if ( ier /= 0 .and. ier /= 7 ) then
      exit
    end if

40  continue

    ll = numrl2
    c1 = c2
    c2 = c2+cycle

  end do
!
!  Set final result and error estimate.
!
!60 continue

  abserr = abserr + 1.0e+01 * correc

  if ( ier == 0 ) then
    return
  end if

  if (is_not_zero(result) .and. is_not_zero(psum(numrl2))) go to 70

  if ( abserr > errsum ) then
    go to 80
  end if

  if (is_zero(psum(numrl2))) then
    return
  end if

70 continue

  if ( abserr / abs(result) <= (errsum+drl)/abs(psum(numrl2)) ) then

    if ( ier >= 1 .and. ier /= 7 ) then
      abserr = abserr + drl
    end if

    return

  end if

80 continue

  result = psum(numrl2)
  abserr = errsum + drl

  return
end subroutine
subroutine qawo ( f, a, b, omega, integr, epsabs, epsrel, result, abserr, &
  neval, ier )

!*****************************************************************************80
!
!! QAWO computes the integrals of oscillatory integrands.
!
!  Discussion:
!
!    The routine calculates an approximation RESULT to a given
!    definite integral
!      I = Integral ( A <= X <= B ) F(X) * cos ( OMEGA * X ) dx
!    or 
!      I = Integral ( A <= X <= B ) F(X) * sin ( OMEGA * X ) dx
!    hopefully satisfying following claim for accuracy
!      | I - RESULT | <= max ( epsabs, epsrel * |I| ).
!
!  Author:
!
!    Robert Piessens, Elise de Doncker-Kapenger, 
!    Christian Ueberhuber, David Kahaner
!
!  Reference:
!
!    Robert Piessens, Elise de Doncker-Kapenger, 
!    Christian Ueberhuber, David Kahaner,
!    QUADPACK, a Subroutine Package for Automatic Integration,
!    Springer Verlag, 1983
!
!  Parameters:
!
!    Input, external real F, the name of the function routine, of the form
!      function f ( x )
!      real f
!      real x
!    which evaluates the integrand function.
!
!    Input, real A, B, the limits of integration.
!
!    Input, real OMEGA, the parameter in the weight function.
!
!    Input, integer INTEGR, specifies the weight function:
!    1, W(X) = cos ( OMEGA * X )
!    2, W(X) = sin ( OMEGA * X )
!
!    Input, real EPSABS, EPSREL, the absolute and relative accuracy requested.
!
!    Output, real RESULT, the estimated value of the integral.
!
!    Output, real ABSERR, an estimate of || I - RESULT ||.
!
!    Output, integer NEVAL, the number of times the integral was evaluated.
!
!            ier    - integer
!                     ier = 0 normal and reliable termination of the
!                             routine. it is assumed that the
!                             requested accuracy has been achieved.
!                   - ier > 0 abnormal termination of the routine.
!                             the estimates for integral and error are
!                             less reliable. it is assumed that the
!                             requested accuracy has not been achieved.
!                     ier = 1 maximum number of subdivisions allowed
!                             (= leniw/2) has been achieved. one can
!                             allow more subdivisions by increasing the
!                             value of leniw (and taking the according
!                             dimension adjustments into account).
!                             however, if this yields no improvement it
!                             is advised to analyze the integrand in
!                             order to determine the integration
!                             difficulties. if the position of a local
!                             difficulty can be determined (e.g.
!                             singularity, discontinuity within the
!                             interval) one will probably gain from
!                             splitting up the interval at this point
!                             and calling the integrator on the
!                             subranges. if possible, an appropriate
!                             special-purpose integrator should
!                             be used which is designed for handling
!                             the type of difficulty involved.
!                         = 2 the occurrence of roundoff error is
!                             detected, which prevents the requested
!                             tolerance from being achieved.
!                             the error may be under-estimated.
!                         = 3 extremely bad integrand behavior occurs
!                             at some interior points of the integration
!                             interval.
!                         = 4 the algorithm does not converge. roundoff
!                             error is detected in the extrapolation
!                             table. it is presumed that the requested
!                             tolerance cannot be achieved due to
!                             roundoff in the extrapolation table,
!                             and that the returned result is the best
!                             which can be obtained.
!                         = 5 the integral is probably divergent, or
!                             slowly convergent. it must be noted that
!                             divergence can occur with any other value
!                             of ier.
!                         = 6 the input is invalid, because
!                             epsabs < 0 and epsrel < 0,
!                             result, abserr, neval are set to zero.
!
!  Local parameters:
!
!    limit is the maximum number of subintervals allowed in the
!    subdivision process of QFOUR. take care that limit >= 1.
!
!    maxp1 gives an upper bound on the number of Chebyshev moments
!    which can be stored, i.e. for the intervals of lengths
!    abs(b-a)*2**(-l), l = 0, 1, ... , maxp1-2. take care that
!    maxp1 >= 1.

  implicit none

  integer, parameter :: limit = 500
  integer, parameter :: maxp1 = 21

  real a
  real abserr
  real alist(limit)
  real b
  real blist(limit)
  real chebmo(maxp1,25)
  real elist(limit)
  real epsabs
  real epsrel
  procedure(scalar_func) :: f
  integer ier
  integer integr
  integer iord(limit)
  integer momcom
  integer neval
  integer nnlog(limit)
  real omega
  real result
  real rlist(limit)

  call qfour ( f, a, b, omega, integr, epsabs, epsrel, limit, 1, maxp1, &
    result, abserr, neval, ier, alist, blist, rlist, elist, iord, nnlog, &
    momcom, chebmo )

  return
end subroutine
subroutine qaws ( f, a, b, alfa, beta, integr, epsabs, epsrel, result, &
  abserr, neval, ier )

!*****************************************************************************80
!
!! QAWS estimates integrals with algebraico-logarithmic endpoint singularities.
!
!  Discussion:
!
!    This routine calculates an approximation RESULT to a given
!    definite integral   
!      I = integral of f*w over (a,b) 
!    where w shows a singular behavior at the end points, see parameter
!    integr, hopefully satisfying following claim for accuracy
!      abs(i-result) <= max(epsabs,epsrel*abs(i)).
!
!  Author:
!
!    Robert Piessens, Elise de Doncker-Kapenger, 
!    Christian Ueberhuber, David Kahaner
!
!  Reference:
!
!    Robert Piessens, Elise de Doncker-Kapenger, 
!    Christian Ueberhuber, David Kahaner,
!    QUADPACK, a Subroutine Package for Automatic Integration,
!    Springer Verlag, 1983
!
!  Parameters:
!
!    Input, external real F, the name of the function routine, of the form
!      function f ( x )
!      real f
!      real x
!    which evaluates the integrand function.
!
!    Input, real A, B, the limits of integration.
!
!    Input, real ALFA, BETA, parameters used in the weight function.
!    ALFA and BETA should be greater than -1.
!
!    Input, integer INTEGR, indicates which weight function is to be used
!    = 1  (x-a)**alfa*(b-x)**beta
!    = 2  (x-a)**alfa*(b-x)**beta*log(x-a)
!    = 3  (x-a)**alfa*(b-x)**beta*log(b-x)
!    = 4  (x-a)**alfa*(b-x)**beta*log(x-a)*log(b-x)
!
!    Input, real EPSABS, EPSREL, the absolute and relative accuracy requested.
!
!    Output, real RESULT, the estimated value of the integral.
!
!    Output, real ABSERR, an estimate of || I - RESULT ||.
!
!    Output, integer NEVAL, the number of times the integral was evaluated.
!
!            ier    - integer
!                     ier = 0 normal and reliable termination of the
!                             routine. it is assumed that the requested
!                             accuracy has been achieved.
!                     ier > 0 abnormal termination of the routine
!                             the estimates for the integral and error
!                             are less reliable. it is assumed that the
!                             requested accuracy has not been achieved.
!                     ier = 1 maximum number of subdivisions allowed
!                             has been achieved. one can allow more
!                             subdivisions by increasing the data value
!                             of limit in qaws (and taking the according
!                             dimension adjustments into account).
!                             however, if this yields no improvement it
!                             is advised to analyze the integrand, in
!                             order to determine the integration
!                             difficulties which prevent the requested
!                             tolerance from being achieved. in case of
!                             a jump discontinuity or a local
!                             singularity of algebraico-logarithmic type
!                             at one or more interior points of the
!                             integration range, one should proceed by
!                             splitting up the interval at these points
!                             and calling the integrator on the
!                             subranges.
!                         = 2 the occurrence of roundoff error is
!                             detected, which prevents the requested
!                             tolerance from being achieved.
!                         = 3 extremely bad integrand behavior occurs
!                             at some points of the integration
!                             interval.
!                         = 6 the input is invalid, because
!                             b <= a or alfa <= (-1) or beta <= (-1) or
!                             integr < 1 or integr > 4 or
!                             epsabs < 0 and epsrel < 0,
!                             result, abserr, neval are set to zero.
!
!  Local parameters:
!
!    LIMIT is the maximum number of subintervals allowed in the
!    subdivision process of qawse. take care that limit >= 2.
!
  implicit none

  integer, parameter :: limit = 500

  real a
  real abserr
  real alfa
  real alist(limit)
  real b
  real blist(limit)
  real beta
  real elist(limit)
  real epsabs
  real epsrel
  procedure(scalar_func) :: f
  integer ier
  integer integr
  integer iord(limit)
  integer last
  integer neval
  real result
  real rlist(limit)

  call qawse ( f, a, b, alfa, beta, integr, epsabs, epsrel, limit, result, &
    abserr, neval, ier, alist, blist, rlist, elist, iord, last )

  return
end subroutine
subroutine qawse ( f, a, b, alfa, beta, integr, epsabs, epsrel, limit, &
  result, abserr, neval, ier, alist, blist, rlist, elist, iord, last )

!*****************************************************************************80
!
!! QAWSE estimates integrals with algebraico-logarithmic endpoint singularities.
!
!  Discussion:
!
!    This routine calculates an approximation RESULT to an integral
!      I = integral of F(X) * W(X) over (a,b), 
!    where W(X) shows a singular behavior at the endpoints, hopefully 
!    satisfying:
!      | I - RESULT | <= max ( epsabs, epsrel * |I| ).
!
!  Author:
!
!    Robert Piessens, Elise de Doncker-Kapenger, 
!    Christian Ueberhuber, David Kahaner
!
!  Reference:
!
!    Robert Piessens, Elise de Doncker-Kapenger, 
!    Christian Ueberhuber, David Kahaner,
!    QUADPACK, a Subroutine Package for Automatic Integration,
!    Springer Verlag, 1983
!
!  Parameters:
!
!    Input, external real F, the name of the function routine, of the form
!      function f ( x )
!      real f
!      real x
!    which evaluates the integrand function.
!
!    Input, real A, B, the limits of integration.
!
!    Input, real ALFA, BETA, parameters used in the weight function.
!    ALFA and BETA should be greater than -1.
!
!    Input, integer INTEGR, indicates which weight function is used:
!    = 1  (x-a)**alfa*(b-x)**beta
!    = 2  (x-a)**alfa*(b-x)**beta*log(x-a)
!    = 3  (x-a)**alfa*(b-x)**beta*log(b-x)
!    = 4  (x-a)**alfa*(b-x)**beta*log(x-a)*log(b-x)
!
!    Input, real EPSABS, EPSREL, the absolute and relative accuracy requested.
!
!    Input, integer LIMIT, an upper bound on the number of subintervals
!    in the partition of (A,B), LIMIT >= 2.  If LIMIT < 2, the routine 
!     will end with IER = 6.
!
!    Output, real RESULT, the estimated value of the integral.
!
!    Output, real ABSERR, an estimate of || I - RESULT ||.
!
!    Output, integer NEVAL, the number of times the integral was evaluated.
!
!            ier    - integer
!                     ier = 0 normal and reliable termination of the
!                             routine. it is assumed that the requested
!                             accuracy has been achieved.
!                     ier > 0 abnormal termination of the routine
!                             the estimates for the integral and error
!                             are less reliable. it is assumed that the
!                             requested accuracy has not been achieved.
!                         = 1 maximum number of subdivisions allowed
!                             has been achieved. one can allow more
!                             subdivisions by increasing the value of
!                             limit. however, if this yields no
!                             improvement it is advised to analyze the
!                             integrand, in order to determine the
!                             integration difficulties which prevent
!                             the requested tolerance from being
!                             achieved. in case of a jump discontinuity
!                             or a local singularity of algebraico-
!                             logarithmic type at one or more interior
!                             points of the integration range, one
!                             should proceed by splitting up the
!                             interval at these points and calling the
!                             integrator on the subranges.
!                         = 2 the occurrence of roundoff error is
!                             detected, which prevents the requested
!                             tolerance from being achieved.
!                         = 3 extremely bad integrand behavior occurs
!                             at some points of the integration
!                             interval.
!                         = 6 the input is invalid, because
!                             b <= a or alfa <= (-1) or beta <= (-1) or
!                             integr < 1 or integr > 4, or
!                             epsabs < 0 and epsrel < 0,
!                             or limit < 2.
!                             result, abserr, neval, rlist(1), elist(1),
!                             iord(1) and last are set to zero.
!                             alist(1) and blist(1) are set to a and b
!                             respectively.
!
!    Workspace, real ALIST(LIMIT), BLIST(LIMIT), contains in entries 1 
!    through LAST the left and right ends of the partition subintervals.
!
!    Workspace, real RLIST(LIMIT), contains in entries 1 through LAST
!    the integral approximations on the subintervals.
!
!    Workspace, real ELIST(LIMIT), contains in entries 1 through LAST
!    the absolute error estimates on the subintervals.
!
!            iord   - integer
!                     vector of dimension at least limit, the first k
!                     elements of which are pointers to the error
!                     estimates over the subintervals, so that
!                     elist(iord(1)), ..., elist(iord(k)) with k = last
!                     if last <= (limit/2+2), and k = limit+1-last
!                     otherwise, form a decreasing sequence.
!
!    Output, integer LAST, the number of subintervals actually produced in 
!    the subdivision process.
!
!  Local parameters:
!
!           alist     - list of left end points of all subintervals
!                       considered up to now
!           blist     - list of right end points of all subintervals
!                       considered up to now
!           rlist(i)  - approximation to the integral over
!                       (alist(i),blist(i))
!           elist(i)  - error estimate applying to rlist(i)
!           maxerr    - pointer to the interval with largest error
!                       estimate
!           errmax    - elist(maxerr)
!           area      - sum of the integrals over the subintervals
!           errsum    - sum of the errors over the subintervals
!           errbnd    - requested accuracy max(epsabs,epsrel*
!                       abs(result))
!           *****1    - variable for the left subinterval
!           *****2    - variable for the right subinterval
!           last      - index for subdivision
!
  implicit none

  integer limit

  real a
  real abserr
  real alfa
  real alist(limit)
  real area
  real area1
  real area12
  real area2
  real a1
  real a2
  real b
  real beta
  real blist(limit)
  real b1
  real b2
  real centre
  real elist(limit)
  real epsabs
  real epsrel
  real errbnd
  real errmax
  real error1
  real erro12
  real error2
  real errsum
  procedure(scalar_func) :: f
  integer ier
  integer integr
  integer iord(limit)
  integer iroff1
  integer iroff2
  integer last
  integer maxerr
  integer nev
  integer neval
  integer nrmax
  real resas1
  real resas2
  real result
  real rg(25)
  real rh(25)
  real ri(25)
  real rj(25)
  real rlist(limit)
!
!  Test on validity of parameters.
!
  ier = 0
  neval = 0
  last = 0
  rlist(1) = 0.0e+00
  elist(1) = 0.0e+00
  iord(1) = 0
  result = 0.0e+00
  abserr = 0.0e+00

  if ( b <= a .or. &
    (epsabs < 0.0e+00 .and. epsrel < 0.0e+00) .or. &
    alfa <= (-1.0e+00) .or. &
    beta <= (-1.0e+00) .or. &
    integr < 1 .or. &
    integr > 4 .or. &
    limit < 2 ) then
    ier = 6
    return
  end if
!
!  Compute the modified Chebyshev moments.
!
  call qmomo ( alfa, beta, ri, rj, rg, rh, integr )
!
!  Integrate over the intervals (a,(a+b)/2) and ((a+b)/2,b).
!
  centre = 5.0e-01 * ( b + a )

  call qc25s ( f, a, b, a, centre, alfa, beta, ri, rj, rg, rh, area1, &
    error1, resas1, integr, nev )

  neval = nev

  call qc25s ( f, a, b, centre, b, alfa, beta, ri, rj, rg, rh, area2, &
    error2, resas2, integr, nev )

  last = 2
  neval = neval+nev
  result = area1+area2
  abserr = error1+error2
!
!  Test on accuracy.
!
  errbnd = max ( epsabs, epsrel * abs ( result ) )
!
!  Initialization.
!
  if ( error2 <= error1 ) then
    alist(1) = a
    alist(2) = centre
    blist(1) = centre
    blist(2) = b
    rlist(1) = area1
    rlist(2) = area2
    elist(1) = error1
    elist(2) = error2
  else
    alist(1) = centre
    alist(2) = a
    blist(1) = b
    blist(2) = centre
    rlist(1) = area2
    rlist(2) = area1
    elist(1) = error2
    elist(2) = error1
  end if

  iord(1) = 1
  iord(2) = 2

  if ( limit == 2 ) then
    ier = 1
    return
  end if

  if ( abserr <= errbnd ) then
    return
  end if

  errmax = elist(1)
  maxerr = 1
  nrmax = 1
  area = result
  errsum = abserr
  iroff1 = 0
  iroff2 = 0

  do last = 3, limit
!
!  Bisect the subinterval with largest error estimate.
!
    a1 = alist(maxerr)
    b1 = 5.0e-01 * ( alist(maxerr) + blist(maxerr) )
    a2 = b1
    b2 = blist(maxerr)

    call qc25s ( f, a, b, a1, b1, alfa, beta, ri, rj, rg, rh, area1, &
      error1, resas1, integr, nev )

    neval = neval + nev

    call qc25s ( f, a, b, a2, b2, alfa, beta, ri, rj, rg, rh, area2, &
      error2, resas2, integr, nev )

    neval = neval + nev
!
!  Improve previous approximations integral and error and
!  test for accuracy.
!
    area12 = area1+area2
    erro12 = error1+error2
    errsum = errsum+erro12-errmax
    area = area+area12-rlist(maxerr)
!
!  Test for roundoff error.
!
    if (not_exactly_equal(a, a1) .and. not_exactly_equal(b, b2)) then
      if (not_exactly_equal(resas1, error1) .and. not_exactly_equal(resas2, error2)) then

        if ( abs ( rlist(maxerr) - area12 ) < 1.0e-05 * abs ( area12 ) &
          .and.erro12 >= 9.9e-01*errmax ) then
          iroff1 = iroff1 + 1
        end if

        if ( last > 10 .and. erro12 > errmax ) then
          iroff2 = iroff2 + 1
        end if

      end if

    end if

    rlist(maxerr) = area1
    rlist(last) = area2
!
!  Test on accuracy.
!
    errbnd = max ( epsabs, epsrel * abs ( area ) )

    if ( errsum > errbnd ) then
!
!  Set error flag in the case that the number of interval
!  bisections exceeds limit.
!
      if ( last == limit ) then
        ier = 1
      end if
!
!  Set error flag in the case of roundoff error.
!
      if ( iroff1 >= 6 .or. iroff2 >= 20 ) then
        ier = 2
     end if
!
!  Set error flag in the case of bad integrand behavior
!  at interior points of integration range.
!
      if ( max ( abs(a1),abs(b2)) <= (1.0e+00+1.0e+03* epsilon ( a1 ) )* &
        ( abs(a2) + 1.0e+03* tiny ( a2) )) then
        ier = 3
      end if

    end if
!
!  Append the newly-created intervals to the list.
!
    if ( error2 <= error1 ) then
      alist(last) = a2
      blist(maxerr) = b1
      blist(last) = b2
      elist(maxerr) = error1
      elist(last) = error2
    else
      alist(maxerr) = a2
      alist(last) = a1
      blist(last) = b1
      rlist(maxerr) = area2
      rlist(last) = area1
      elist(maxerr) = error2
      elist(last) = error1
    end if
!
!  Call QSORT to maintain the descending ordering
!  in the list of error estimates and select the subinterval
!  with largest error estimate (to be bisected next).
!
    call qsort ( limit, last, maxerr, errmax, elist, iord, nrmax )

    if ( ier /= 0 .or. errsum <= errbnd ) then
      exit
    end if

  end do
!
!  Compute final result.
!
  result = sum ( rlist(1:last) )

  abserr = errsum

  return
end subroutine
subroutine qc25c ( f, a, b, c, result, abserr, krul, neval )

!*****************************************************************************80
!
!! QC25C returns integration rules for Cauchy Principal Value integrals.
!
!  Discussion:
!
!    This routine estimates 
!      I = integral of F(X) * W(X) over (a,b) 
!    with error estimate, where 
!      w(x) = 1/(x-c)
!
!  Author:
!
!    Robert Piessens, Elise de Doncker-Kapenger, 
!    Christian Ueberhuber, David Kahaner
!
!  Reference:
!
!    Robert Piessens, Elise de Doncker-Kapenger, 
!    Christian Ueberhuber, David Kahaner,
!    QUADPACK, a Subroutine Package for Automatic Integration,
!    Springer Verlag, 1983
!
!  Parameters:
!
!    Input, external real F, the name of the function routine, of the form
!      function f ( x )
!      real f
!      real x
!    which evaluates the integrand function.
!
!    Input, real A, B, the limits of integration.
!
!    Input, real C, the parameter in the weight function.
!
!    Output, real RESULT, the estimated value of the integral.
!    RESULT is computed by using a generalized Clenshaw-Curtis method if
!    C lies within ten percent of the integration interval.  In the 
!    other case the 15-point Kronrod rule obtained by optimal addition
!    of abscissae to the 7-point Gauss rule, is applied.
!
!    Output, real ABSERR, an estimate of || I - RESULT ||.
!
!           krul   - integer
!                    key which is decreased by 1 if the 15-point
!                    Gauss-Kronrod scheme has been used
!
!    Output, integer NEVAL, the number of times the integral was evaluated.
!
!  Local parameters:
!
!           fval   - value of the function f at the points
!                    cos(k*pi/24),  k = 0, ..., 24
!           cheb12 - Chebyshev series expansion coefficients, for the
!                    function f, of degree 12
!           cheb24 - Chebyshev series expansion coefficients, for the
!                    function f, of degree 24
!           res12  - approximation to the integral corresponding to the
!                    use of cheb12
!           res24  - approximation to the integral corresponding to the
!                    use of cheb24
!           qwgtc  - external function subprogram defining the weight
!                    function
!           hlgth  - half-length of the interval
!           centr  - mid point of the interval
!
  implicit none

  real a
  real abserr
  real ak22
  real amom0
  real amom1
  real amom2
  real b
  real c
  real cc
  real centr
  real cheb12(13)
  real cheb24(25)
  procedure(scalar_func) :: f
  real fval(25)
  real hlgth
  integer i
  integer isym
  integer k 
  integer kp
  integer krul
  integer neval
  real p2
  real p3
  real p4
  real resabs
  real resasc
  real result
  real res12
  real res24
  real u
  real, parameter, dimension ( 11 ) :: x = (/ &
    9.914448613738104e-01, 9.659258262890683e-01, &
    9.238795325112868e-01, 8.660254037844386e-01, &
    7.933533402912352e-01, 7.071067811865475e-01, &
    6.087614290087206e-01, 5.000000000000000e-01, &
    3.826834323650898e-01, 2.588190451025208e-01, &
    1.305261922200516e-01 /)
!
!  Check the position of C.
!
  cc = ( 2.0e+00 * c - b - a ) / ( b - a )
!
!  Apply the 15-point Gauss-Kronrod scheme.
!
  if ( abs ( cc ) >= 1.1e+00 ) then
    krul = krul - 1
    call qk15w ( f, qwgtc, c, p2, p3, p4, kp, a, b, result, abserr, &
      resabs, resasc )
    neval = 15
    if (exactly_equal(resasc, abserr)) then
      krul = krul+1
    end if
    return
  end if
!
!  Use the generalized Clenshaw-Curtis method.
!
  hlgth = 5.0e-01 * ( b - a )
  centr = 5.0e-01 * ( b + a )
  neval = 25
  fval(1) = 5.0e-01 * f(hlgth+centr)
  fval(13) = f(centr)
  fval(25) = 5.0e-01 * f(centr-hlgth)

  do i = 2, 12
    u = hlgth * x(i-1)
    isym = 26 - i
    fval(i) = f(u+centr)
    fval(isym) = f(centr-u)
  end do
!
!  Compute the Chebyshev series expansion.
!
  call qcheb ( x, fval, cheb12, cheb24 )
!
!  The modified Chebyshev moments are computed by forward
!  recursion, using AMOM0 and AMOM1 as starting values.
!
  amom0 = log ( abs ( ( 1.0e+00 - cc ) / ( 1.0e+00 + cc ) ) )
  amom1 = 2.0e+00 + cc * amom0
  res12 = cheb12(1) * amom0 + cheb12(2) * amom1
  res24 = cheb24(1) * amom0 + cheb24(2) * amom1

  do k = 3, 13
    amom2 = 2.0e+00 * cc * amom1 - amom0
    ak22 = ( k - 2 ) * ( k - 2 )
    if ( ( k / 2 ) * 2 == k ) then
      amom2 = amom2 - 4.0e+00 / ( ak22 - 1.0e+00 )
    end if
    res12 = res12 + cheb12(k) * amom2
    res24 = res24 + cheb24(k) * amom2
    amom0 = amom1
    amom1 = amom2
  end do

  do k = 14, 25
    amom2 = 2.0e+00 * cc * amom1 - amom0
    ak22 = ( k - 2 ) * ( k - 2 )
    if ( ( k / 2 ) * 2 == k ) then
      amom2 = amom2 - 4.0e+00 / ( ak22 - 1.0e+00 )
    end if
    res24 = res24 + cheb24(k) * amom2
    amom0 = amom1
    amom1 = amom2
  end do

  result = res24
  abserr = abs ( res24 - res12 )

  return
end subroutine
subroutine qc25o ( f, a, b, omega, integr, nrmom, maxp1, ksave, result, &
  abserr, neval, resabs, resasc, momcom, chebmo )

!*****************************************************************************80
!
!! QC25O returns integration rules for integrands with a COS or SIN factor.
!
!  Discussion:
!
!    This routine estimates the integral
!      I = integral of f(x) * w(x) over (a,b)
!    where
!      w(x) = cos(omega*x)
!    or 
!      w(x) = sin(omega*x),
!    and estimates
!      J = integral ( A <= X <= B ) |F(X)| dx.
!
!    For small values of OMEGA or small intervals (a,b) the 15-point
!    Gauss-Kronrod rule is used.  In all other cases a generalized
!    Clenshaw-Curtis method is used, that is, a truncated Chebyshev 
!    expansion of the function F is computed on (a,b), so that the 
!    integrand can be written as a sum of terms of the form W(X)*T(K,X), 
!    where T(K,X) is the Chebyshev polynomial of degree K.  The Chebyshev
!    moments are computed with use of a linear recurrence relation.
!
!  Author:
!
!    Robert Piessens, Elise de Doncker-Kapenger, 
!    Christian Ueberhuber, David Kahaner
!
!  Reference:
!
!    Robert Piessens, Elise de Doncker-Kapenger, 
!    Christian Ueberhuber, David Kahaner,
!    QUADPACK, a Subroutine Package for Automatic Integration,
!    Springer Verlag, 1983
!
!  Parameters:
!
!    Input, external real F, the name of the function routine, of the form
!      function f ( x )
!      real f
!      real x
!    which evaluates the integrand function.
!
!    Input, real A, B, the limits of integration.
!
!    Input, real OMEGA, the parameter in the weight function.
!
!    Input, integer INTEGR, indicates which weight function is to be used
!    = 1, w(x) = cos(omega*x)
!    = 2, w(x) = sin(omega*x)
!
!    ?, integer NRMOM, the length of interval (a,b) is equal to the length
!    of the original integration interval divided by
!    2**nrmom (we suppose that the routine is used in an
!    adaptive integration process, otherwise set
!    nrmom = 0).  nrmom must be zero at the first call.
!
!           maxp1  - integer
!                    gives an upper bound on the number of Chebyshev
!                    moments which can be stored, i.e. for the intervals
!                    of lengths abs(bb-aa)*2**(-l), l = 0,1,2, ...,
!                    maxp1-2.
!
!           ksave  - integer
!                    key which is one when the moments for the
!                    current interval have been computed
!
!    Output, real RESULT, the estimated value of the integral.
!
!           abserr - real
!                    estimate of the modulus of the absolute
!                    error, which should equal or exceed abs(i-result)
!
!    Output, integer NEVAL, the number of times the integral was evaluated.
!
!    Output, real RESABS, approximation to the integral J.
!
!    Output, real RESASC, approximation to the integral of abs(F-I/(B-A)).
!
!         on entry and return
!           momcom - integer
!                    for each interval length we need to compute
!                    the Chebyshev moments. momcom counts the number
!                    of intervals for which these moments have already
!                    been computed. if nrmom < momcom or ksave = 1,
!                    the Chebyshev moments for the interval (a,b)
!                    have already been computed and stored, otherwise
!                    we compute them and we increase momcom.
!
!           chebmo - real
!                    array of dimension at least (maxp1,25) containing
!                    the modified Chebyshev moments for the first momcom
!                    interval lengths
!
!  Local parameters:
!
!    maxp1 gives an upper bound
!           on the number of Chebyshev moments which can be
!           computed, i.e. for the interval (bb-aa), ...,
!           (bb-aa)/2**(maxp1-2).
!           should this number be altered, the first dimension of
!           chebmo needs to be adapted.
!
!    x contains the values cos(k*pi/24)
!           k = 1, ...,11, to be used for the Chebyshev expansion of f
!
!           centr  - mid point of the integration interval
!           hlgth  - half length of the integration interval
!           fval   - value of the function f at the points
!                    (b-a)*0.5*cos(k*pi/12) + (b+a)*0.5
!                    k = 0, ...,24
!           cheb12 - coefficients of the Chebyshev series expansion
!                    of degree 12, for the function f, in the
!                    interval (a,b)
!           cheb24 - coefficients of the Chebyshev series expansion
!                    of degree 24, for the function f, in the
!                    interval (a,b)
!           resc12 - approximation to the integral of
!                    cos(0.5*(b-a)*omega*x)*f(0.5*(b-a)*x+0.5*(b+a))
!                    over (-1,+1), using the Chebyshev series
!                    expansion of degree 12
!           resc24 - approximation to the same integral, using the
!                    Chebyshev series expansion of degree 24
!           ress12 - the analogue of resc12 for the sine
!           ress24 - the analogue of resc24 for the sine
!
  implicit none

  integer maxp1

  real a
  real abserr
  real ac
  real an
  real an2
  real as
  real asap
  real ass
  real b
  real centr
  real chebmo(maxp1,25)
  real cheb12(13)
  real cheb24(25)
  real conc
  real cons
  real cospar
  real d(28)
  real d1(28)
  real d2(28)
  real d3(28)
  real estc
  real ests
  procedure(scalar_func) :: f
  real fval(25)
  real hlgth
  integer i
  integer integr
  integer isym
  integer j
  integer k
  integer ksave
  integer m
  integer momcom
  integer neval
  integer, parameter :: nmac = 28
  integer noeq1
  integer noequ
  integer nrmom
  real omega
  real parint
  real par2
  real par22
  real p2
  real p3
  real p4
  real resabs
  real resasc
  real resc12
  real resc24
  real ress12
  real ress24
  real result
  real sinpar
  real v(28)
  real, dimension ( 11 ) :: x = (/ &
    9.914448613738104e-01,     9.659258262890683e-01, &
    9.238795325112868e-01,     8.660254037844386e-01, &
    7.933533402912352e-01,     7.071067811865475e-01, &
    6.087614290087206e-01,     5.000000000000000e-01, &
    3.826834323650898e-01,     2.588190451025208e-01, &
    1.305261922200516e-01 /)

  centr = 5.0e-01 * ( b + a )
  hlgth = 5.0e-01 * ( b - a )
  parint = omega * hlgth
!
!  Compute the integral using the 15-point Gauss-Kronrod
!  formula if the value of the parameter in the integrand
!  is small or if the length of the integration interval
!  is less than (bb-aa)/2**(maxp1-2), where (aa,bb) is the
!  original integration interval.
!
  if ( abs ( parint ) <= 2.0e+00 ) then

    call qk15w ( f, qwgto, omega, p2, p3, p4, integr, a, b, result, &
      abserr, resabs, resasc )

    neval = 15
    return

  end if
!
!  Compute the integral using the generalized clenshaw-curtis method.
!
  conc = hlgth * cos(centr*omega)
  cons = hlgth * sin(centr*omega)
  resasc = huge ( resasc )
  neval = 25
!
!  Check whether the Chebyshev moments for this interval
!  have already been computed.
!
  if ( nrmom < momcom .or. ksave == 1 ) then
    go to 140
  end if
!
!  Compute a new set of Chebyshev moments.
!
  m = momcom + 1
  par2 = parint * parint
  par22 = par2 + 2.0e+00
  sinpar = sin(parint)
  cospar = cos(parint)
!
!  Compute the Chebyshev moments with respect to cosine.
!
  v(1) = 2.0e+00 * sinpar / parint
  v(2) = (8.0e+00*cospar+(par2+par2-8.0e+00)*sinpar/ parint)/par2
  v(3) = (3.2e+01*(par2-1.2e+01)*cospar+(2.0e+00* &
     ((par2-8.0e+01)*par2+1.92e+02)*sinpar)/ &
     parint)/(par2*par2)
  ac = 8.0e+00*cospar
  as = 2.4e+01*parint*sinpar

  if ( abs ( parint ) > 2.4e+01 ) then
    go to 70
  end if
!
!  Compute the Chebyshev moments as the solutions of a boundary value 
!  problem with one initial value (v(3)) and one end value computed
!  using an asymptotic formula.
!
  noequ = nmac-3
  noeq1 = noequ-1
  an = 6.0e+00

  do k = 1, noeq1
    an2 = an*an
    d(k) = -2.0e+00*(an2-4.0e+00) * (par22-an2-an2)
    d2(k) = (an-1.0e+00)*(an-2.0e+00) * par2
    d1(k) = (an+3.0e+00)*(an+4.0e+00) * par2
    v(k+3) = as-(an2-4.0e+00) * ac
    an = an+2.0e+00
  end do

  an2 = an*an
  d(noequ) = -2.0e+00*(an2-4.0e+00) * (par22-an2-an2)
  v(noequ+3) = as - ( an2 - 4.0e+00 ) * ac
  v(4) = v(4) - 5.6e+01 * par2 * v(3)
  ass = parint * sinpar
  asap = (((((2.10e+02*par2-1.0e+00)*cospar-(1.05e+02*par2 &
     -6.3e+01)*ass)/an2-(1.0e+00-1.5e+01*par2)*cospar &
     +1.5e+01*ass)/an2-cospar+3.0e+00*ass)/an2-cospar)/an2
  v(noequ+3) = v(noequ+3)-2.0e+00*asap*par2*(an-1.0e+00)* &
      (an-2.0e+00)
!
!  Solve the tridiagonal system by means of Gaussian
!  elimination with partial pivoting.
!
  d3(1:noequ) = 0.0e+00

  d2(noequ) = 0.0e+00

  do i = 1, noeq1

    if ( abs(d1(i)) > abs(d(i)) ) then
      an = d1(i)
      d1(i) = d(i)
      d(i) = an
      an = d2(i)
      d2(i) = d(i+1)
      d(i+1) = an
      d3(i) = d2(i+1)
      d2(i+1) = 0.0e+00
      an = v(i+4)
      v(i+4) = v(i+3)
      v(i+3) = an
    end if

    d(i+1) = d(i+1)-d2(i)*d1(i)/d(i)
    d2(i+1) = d2(i+1)-d3(i)*d1(i)/d(i)
    v(i+4) = v(i+4)-v(i+3)*d1(i)/d(i)

  end do

  v(noequ+3) = v(noequ+3) / d(noequ)
  v(noequ+2) = (v(noequ+2)-d2(noeq1)*v(noequ+3))/d(noeq1)

  do i = 2, noeq1
    k = noequ-i
    v(k+3) = (v(k+3)-d3(k)*v(k+5)-d2(k)*v(k+4))/d(k)
  end do

  go to 90
!
!  Compute the Chebyshev moments by means of forward recursion
!
70 continue

  an = 4.0e+00

  do i = 4, 13
    an2 = an*an
    v(i) = ((an2-4.0e+00)*(2.0e+00*(par22-an2-an2)*v(i-1)-ac) &
     +as-par2*(an+1.0e+00)*(an+2.0e+00)*v(i-2))/ &
     (par2*(an-1.0e+00)*(an-2.0e+00))
    an = an+2.0e+00
  end do

90 continue

  do j = 1, 13
    chebmo(m,2*j-1) = v(j)
  end do
!
!  Compute the Chebyshev moments with respect to sine.
!
  v(1) = 2.0e+00*(sinpar-parint*cospar)/par2
  v(2) = (1.8e+01-4.8e+01/par2)*sinpar/par2 &
     +(-2.0e+00+4.8e+01/par2)*cospar/parint
  ac = -2.4e+01*parint*cospar
  as = -8.0e+00*sinpar
  chebmo(m,2) = v(1)
  chebmo(m,4) = v(2)

  if ( abs(parint) <= 2.4e+01 ) then

    do k = 3, 12
      an = k
      chebmo(m,2*k) = -sinpar/(an*(2.0e+00*an-2.0e+00)) &
                 -2.5e-01*parint*(v(k+1)/an-v(k)/(an-1.0e+00))
    end do
!
!  Compute the Chebyshev moments by means of forward recursion.
!
  else

    an = 3.0e+00

    do i = 3, 12
      an2 = an*an
      v(i) = ((an2-4.0e+00)*(2.0e+00*(par22-an2-an2)*v(i-1)+as) &
       +ac-par2*(an+1.0e+00)*(an+2.0e+00)*v(i-2)) &
       /(par2*(an-1.0e+00)*(an-2.0e+00))
      an = an+2.0e+00
      chebmo(m,2*i) = v(i)
    end do

  end if

140 continue

  if ( nrmom < momcom ) then
    m = nrmom + 1
  end if

  if ( momcom < maxp1 - 1 .and. nrmom >= momcom ) then
    momcom = momcom + 1
  end if
!
!  Compute the coefficients of the Chebyshev expansions
!  of degrees 12 and 24 of the function F.
!
  fval(1) = 5.0e-01 * f(centr+hlgth)
  fval(13) = f(centr)
  fval(25) = 5.0e-01 * f(centr-hlgth)

  do i = 2, 12
    isym = 26-i
    fval(i) = f(hlgth*x(i-1)+centr)
    fval(isym) = f(centr-hlgth*x(i-1))
  end do

  call qcheb ( x, fval, cheb12, cheb24 )
!
!  Compute the integral and error estimates.
!
  resc12 = cheb12(13) * chebmo(m,13)
  ress12 = 0.0e+00
  estc = abs ( cheb24(25)*chebmo(m,25))+abs((cheb12(13)- &
    cheb24(13))*chebmo(m,13) )
  ests = 0.0e+00
  k = 11

  do j = 1, 6
    resc12 = resc12+cheb12(k)*chebmo(m,k)
    ress12 = ress12+cheb12(k+1)*chebmo(m,k+1)
    estc = estc+abs((cheb12(k)-cheb24(k))*chebmo(m,k))
    ests = ests+abs((cheb12(k+1)-cheb24(k+1))*chebmo(m,k+1))
    k = k-2
  end do

  resc24 = cheb24(25)*chebmo(m,25)
  ress24 = 0.0e+00
  resabs = abs(cheb24(25))
  k = 23

  do j = 1, 12

    resc24 = resc24+cheb24(k)*chebmo(m,k)
    ress24 = ress24+cheb24(k+1)*chebmo(m,k+1)
    resabs = resabs+abs(cheb24(k))+abs(cheb24(k+1))

    if ( j <= 5 ) then
      estc = estc+abs(cheb24(k)*chebmo(m,k))
      ests = ests+abs(cheb24(k+1)*chebmo(m,k+1))
    end if

    k = k-2

  end do

  resabs = resabs * abs ( hlgth )

  if ( integr == 1 ) then
    result = conc * resc24-cons*ress24
    abserr = abs ( conc * estc ) + abs ( cons * ests )
  else
    result = conc*ress24+cons*resc24
    abserr = abs(conc*ests)+abs(cons*estc)
  end if

  return
end subroutine
subroutine qc25s ( f, a, b, bl, br, alfa, beta, ri, rj, rg, rh, result, &
  abserr, resasc, integr, neval )

!*****************************************************************************80
!
!! QC25S returns rules for algebraico-logarithmic end point singularities.
!
!  Discussion:
!
!    This routine computes 
!      i = integral of F(X) * W(X) over (bl,br), 
!    with error estimate, where the weight function W(X) has a singular
!    behavior of algebraico-logarithmic type at the points
!    a and/or b. 
!
!    The interval (bl,br) is a subinterval of (a,b).
!
!  Author:
!
!    Robert Piessens, Elise de Doncker-Kapenger, 
!    Christian Ueberhuber, David Kahaner
!
!  Reference:
!
!    Robert Piessens, Elise de Doncker-Kapenger, 
!    Christian Ueberhuber, David Kahaner,
!    QUADPACK, a Subroutine Package for Automatic Integration,
!    Springer Verlag, 1983
!
!  Parameters:
!
!    Input, external real F, the name of the function routine, of the form
!      function f ( x )
!      real f
!      real x
!    which evaluates the integrand function.
!
!    Input, real A, B, the limits of integration.
!
!    Input, real BL, BR, the lower and upper limits of integration.
!    A <= BL < BR <= B.
!
!    Input, real ALFA, BETA, parameters in the weight function.
!
!    Input, real RI(25), RJ(25), RG(25), RH(25), modified Chebyshev moments 
!    for the application of the generalized Clenshaw-Curtis method,
!    computed in QMOMO.
!
!    Output, real RESULT, the estimated value of the integral, computed by 
!    using a generalized clenshaw-curtis method if b1 = a or br = b.
!    In all other cases the 15-point Kronrod rule is applied, obtained by
!    optimal addition of abscissae to the 7-point Gauss rule.
!
!    Output, real ABSERR, an estimate of || I - RESULT ||.
!
!    Output, real RESASC, approximation to the integral of abs(F*W-I/(B-A)).
!
!    Input, integer INTEGR,  determines the weight function
!    1, w(x) = (x-a)**alfa*(b-x)**beta
!    2, w(x) = (x-a)**alfa*(b-x)**beta*log(x-a)
!    3, w(x) = (x-a)**alfa*(b-x)**beta*log(b-x)
!    4, w(x) = (x-a)**alfa*(b-x)**beta*log(x-a)*log(b-x)
!
!    Output, integer NEVAL, the number of times the integral was evaluated.
!
!  Local Parameters:
!
!           fval   - value of the function f at the points
!                    (br-bl)*0.5*cos(k*pi/24)+(br+bl)*0.5
!                    k = 0, ..., 24
!           cheb12 - coefficients of the Chebyshev series expansion
!                    of degree 12, for the function f, in the interval
!                    (bl,br)
!           cheb24 - coefficients of the Chebyshev series expansion
!                    of degree 24, for the function f, in the interval
!                    (bl,br)
!           res12  - approximation to the integral obtained from cheb12
!           res24  - approximation to the integral obtained from cheb24
!           qwgts  - external function subprogram defining the four
!                    possible weight functions
!           hlgth  - half-length of the interval (bl,br)
!           centr  - mid point of the interval (bl,br)
!
!           the vector x contains the values cos(k*pi/24)
!           k = 1, ..., 11, to be used for the computation of the
!           Chebyshev series expansion of f.
!
  implicit none

  real a
  real abserr
  real alfa
  real b
  real beta
  real bl
  real br
  real centr
  real cheb12(13)
  real cheb24(25)
  real dc
  procedure(scalar_func) :: f
  real factor
  real fix
  real fval(25)
  real hlgth
  integer i
  integer integr
  integer isym
  integer neval
  real resabs
  real resasc
  real result
  real res12
  real res24
  real rg(25)
  real rh(25)
  real ri(25)
  real rj(25)
  real u
  real, dimension ( 11 ) :: x = (/ &
    9.914448613738104e-01,     9.659258262890683e-01, &
    9.238795325112868e-01,     8.660254037844386e-01, &
    7.933533402912352e-01,     7.071067811865475e-01, &
    6.087614290087206e-01,     5.000000000000000e-01, &
    3.826834323650898e-01,     2.588190451025208e-01, &
    1.305261922200516e-01 /)

  neval = 25

  if (exactly_equal(bl, a) .and. (is_not_zero(alfa) .or. integr == 2 .or. integr == 4)) go to 10
  if (exactly_equal(br, b) .and. (is_not_zero(beta) .or. integr == 3 .or. integr == 4)) go to 140
!
!  If a > bl and b < br, apply the 15-point Gauss-Kronrod scheme.
!
  call qk15w ( f, qwgts, a, b, alfa, beta, integr, bl, br, result, abserr, &
    resabs, resasc )

  neval = 15
  return
!
!  This part of the program is executed only if a = bl.
!
!  Compute the Chebyshev series expansion of the function
!  f1 = (0.5*(b+b-br-a)-0.5*(br-a)*x)**beta*f(0.5*(br-a)*x+0.5*(br+a))
!
10 continue

  hlgth = 5.0e-01*(br-bl)
  centr = 5.0e-01*(br+bl)
  fix = b-centr
  fval(1) = 5.0e-01*f(hlgth+centr)*(fix-hlgth)**beta
  fval(13) = f(centr)*(fix**beta)
  fval(25) = 5.0e-01*f(centr-hlgth)*(fix+hlgth)**beta

  do i = 2, 12
    u = hlgth*x(i-1)
    isym = 26-i
    fval(i) = f(u+centr)*(fix-u)**beta
    fval(isym) = f(centr-u)*(fix+u)**beta
  end do

  factor = hlgth**(alfa+1.0e+00)
  result = 0.0e+00
  abserr = 0.0e+00
  res12 = 0.0e+00
  res24 = 0.0e+00

  if ( integr > 2 ) go to 70

  call qcheb ( x, fval, cheb12, cheb24 )
!
!  integr = 1  (or 2)
!
  do i = 1, 13
    res12 = res12+cheb12(i)*ri(i)
    res24 = res24+cheb24(i)*ri(i)
  end do

  do i = 14, 25
    res24 = res24 + cheb24(i) * ri(i)
  end do

  if ( integr == 1 ) go to 130
!
!  integr = 2
!
  dc = log ( br - bl )
  result = res24 * dc
  abserr = abs((res24-res12)*dc)
  res12 = 0.0e+00
  res24 = 0.0e+00

  do i = 1, 13
    res12 = res12+cheb12(i)*rg(i)
    res24 = res24+cheb24(i)*rg(i)
  end do

  do i = 14, 25
    res24 = res24+cheb24(i)*rg(i)
  end do

  go to 130
!
!  Compute the Chebyshev series expansion of the function
!  F4 = f1*log(0.5*(b+b-br-a)-0.5*(br-a)*x)
!
70 continue

  fval(1) = fval(1) * log ( fix - hlgth )
  fval(13) = fval(13) * log ( fix )
  fval(25) = fval(25) * log ( fix + hlgth )

  do i = 2, 12
    u = hlgth*x(i-1)
    isym = 26-i
    fval(i) = fval(i) * log ( fix - u )
    fval(isym) = fval(isym) * log ( fix + u )
  end do

  call qcheb ( x, fval, cheb12, cheb24 )
!
!  integr = 3  (or 4)
!
  do i = 1, 13
    res12 = res12+cheb12(i)*ri(i)
    res24 = res24+cheb24(i)*ri(i)
  end do

  do i = 14, 25
    res24 = res24+cheb24(i)*ri(i)
  end do

  if ( integr == 3 ) then
    go to 130
  end if
!
!  integr = 4
!
  dc = log ( br - bl )
  result = res24*dc
  abserr = abs((res24-res12)*dc)
  res12 = 0.0e+00
  res24 = 0.0e+00

  do i = 1, 13
    res12 = res12+cheb12(i)*rg(i)
    res24 = res24+cheb24(i)*rg(i)
  end do

  do i = 14, 25
    res24 = res24+cheb24(i)*rg(i)
  end do

130 continue

  result = (result+res24)*factor
  abserr = (abserr+abs(res24-res12))*factor
  go to 270
!
!  This part of the program is executed only if b = br.
!
!  Compute the Chebyshev series expansion of the function
!  f2 = (0.5*(b+bl-a-a)+0.5*(b-bl)*x)**alfa*f(0.5*(b-bl)*x+0.5*(b+bl))
!
140 continue

  hlgth = 5.0e-01*(br-bl)
  centr = 5.0e-01*(br+bl)
  fix = centr-a
  fval(1) = 5.0e-01*f(hlgth+centr)*(fix+hlgth)**alfa
  fval(13) = f(centr)*(fix**alfa)
  fval(25) = 5.0e-01*f(centr-hlgth)*(fix-hlgth)**alfa

  do i = 2, 12
    u = hlgth*x(i-1)
    isym = 26-i
    fval(i) = f(u+centr)*(fix+u)**alfa
    fval(isym) = f(centr-u)*(fix-u)**alfa
  end do

  factor = hlgth**(beta+1.0e+00)
  result = 0.0e+00
  abserr = 0.0e+00
  res12 = 0.0e+00
  res24 = 0.0e+00

  if ( integr == 2 .or. integr == 4 ) then
    go to 200
  end if
!
!  integr = 1  (or 3)
!
  call qcheb ( x, fval, cheb12, cheb24 )

  do i = 1, 13
    res12 = res12+cheb12(i)*rj(i)
    res24 = res24+cheb24(i)*rj(i)
  end do

  do i = 14, 25
    res24 = res24+cheb24(i)*rj(i)
  end do

  if ( integr == 1 ) go to 260
!
!  integr = 3
!
  dc = log ( br - bl )
  result = res24*dc
  abserr = abs((res24-res12)*dc)
  res12 = 0.0e+00
  res24 = 0.0e+00

  do i = 1, 13
    res12 = res12+cheb12(i)*rh(i)
    res24 = res24+cheb24(i)*rh(i)
  end do

  do i = 14, 25
    res24 = res24+cheb24(i)*rh(i)
  end do

  go to 260
!
!  Compute the Chebyshev series expansion of the function
!  f3 = f2*log(0.5*(b-bl)*x+0.5*(b+bl-a-a))
!
200 continue

  fval(1) = fval(1) * log ( hlgth + fix )
  fval(13) = fval(13) * log ( fix )
  fval(25) = fval(25) * log ( fix - hlgth )

  do i = 2, 12
    u = hlgth*x(i-1)
    isym = 26-i
    fval(i) = fval(i) * log(u+fix)
    fval(isym) = fval(isym) * log(fix-u)
  end do

  call qcheb ( x, fval, cheb12, cheb24 )
!
!  integr = 2  (or 4)
!
  do i = 1, 13
    res12 = res12+cheb12(i)*rj(i)
    res24 = res24+cheb24(i)*rj(i)
  end do

  do i = 14, 25
    res24 = res24+cheb24(i)*rj(i)
  end do

  if ( integr == 2 ) go to 260

  dc = log(br-bl)
  result = res24*dc
  abserr = abs((res24-res12)*dc)
  res12 = 0.0e+00
  res24 = 0.0e+00
!
!  integr = 4
!
  do i = 1, 13
    res12 = res12+cheb12(i)*rh(i)
    res24 = res24+cheb24(i)*rh(i)
  end do

  do i = 14, 25
    res24 = res24+cheb24(i)*rh(i)
  end do

260 continue

  result = (result+res24)*factor
  abserr = (abserr+abs(res24-res12))*factor

270 continue

  return
end subroutine
subroutine qcheb ( x, fval, cheb12, cheb24 )

!*****************************************************************************80
!
!! QCHEB computes the Chebyshev series expansion.
!
!  Discussion:
!
!    This routine computes the Chebyshev series expansion
!    of degrees 12 and 24 of a function using a fast Fourier transform method
!
!      f(x) = sum(k=1, ...,13) (cheb12(k)*t(k-1,x)),
!      f(x) = sum(k=1, ...,25) (cheb24(k)*t(k-1,x)),
!
!    where T(K,X) is the Chebyshev polynomial of degree K.
!
!  Author:
!
!    Robert Piessens, Elise de Doncker-Kapenger, 
!    Christian Ueberhuber, David Kahaner
!
!  Reference:
!
!    Robert Piessens, Elise de Doncker-Kapenger, 
!    Christian Ueberhuber, David Kahaner,
!    QUADPACK, a Subroutine Package for Automatic Integration,
!    Springer Verlag, 1983
!
!  Parameters:
!
!    Input, real X(11), contains the values of COS(K*PI/24), for K = 1 to 11.
!
!    Input/output, real FVAL(25), the function values at the points
!    (b+a+(b-a)*cos(k*pi/24))/2, k = 0, ...,24, where (a,b) is the 
!    approximation interval.  FVAL(1) and FVAL(25) are divided by two
!    These values are destroyed at output.
!
!    Output, real CHEB12(13), the Chebyshev coefficients for degree 12.
!
!    Output, real CHEB24(25), the Chebyshev coefficients for degree 24.
!
  implicit none

  real alam
  real alam1
  real alam2
  real cheb12(13)
  real cheb24(25)
  real fval(25)
  integer i
  integer j
  real part1
  real part2
  real part3
  real v(12)
  real x(11)

  do i = 1, 12
    j = 26-i
    v(i) = fval(i)-fval(j)
    fval(i) = fval(i)+fval(j)
  end do

  alam1 = v(1)-v(9)
  alam2 = x(6)*(v(3)-v(7)-v(11))
  cheb12(4) = alam1+alam2
  cheb12(10) = alam1-alam2
  alam1 = v(2)-v(8)-v(10)
  alam2 = v(4)-v(6)-v(12)
  alam = x(3)*alam1+x(9)*alam2
  cheb24(4) = cheb12(4)+alam
  cheb24(22) = cheb12(4)-alam
  alam = x(9)*alam1-x(3)*alam2
  cheb24(10) = cheb12(10)+alam
  cheb24(16) = cheb12(10)-alam
  part1 = x(4)*v(5)
  part2 = x(8)*v(9)
  part3 = x(6)*v(7)
  alam1 = v(1)+part1+part2
  alam2 = x(2)*v(3)+part3+x(10)*v(11)
  cheb12(2) = alam1+alam2
  cheb12(12) = alam1-alam2
  alam = x(1)*v(2)+x(3)*v(4)+x(5)*v(6)+x(7)*v(8) &
    +x(9)*v(10)+x(11)*v(12)
  cheb24(2) = cheb12(2)+alam
  cheb24(24) = cheb12(2)-alam
  alam = x(11)*v(2)-x(9)*v(4)+x(7)*v(6)-x(5)*v(8) &
    +x(3)*v(10)-x(1)*v(12)
  cheb24(12) = cheb12(12)+alam
  cheb24(14) = cheb12(12)-alam
  alam1 = v(1)-part1+part2
  alam2 = x(10)*v(3)-part3+x(2)*v(11)
  cheb12(6) = alam1+alam2
  cheb12(8) = alam1-alam2
  alam = x(5)*v(2)-x(9)*v(4)-x(1)*v(6) &
    -x(11)*v(8)+x(3)*v(10)+x(7)*v(12)
  cheb24(6) = cheb12(6)+alam
  cheb24(20) = cheb12(6)-alam
  alam = x(7)*v(2)-x(3)*v(4)-x(11)*v(6)+x(1)*v(8) &
    -x(9)*v(10)-x(5)*v(12)
  cheb24(8) = cheb12(8)+alam
  cheb24(18) = cheb12(8)-alam

  do i = 1, 6
    j = 14-i
    v(i) = fval(i)-fval(j)
    fval(i) = fval(i)+fval(j)
  end do

  alam1 = v(1)+x(8)*v(5)
  alam2 = x(4)*v(3)
  cheb12(3) = alam1+alam2
  cheb12(11) = alam1-alam2
  cheb12(7) = v(1)-v(5)
  alam = x(2)*v(2)+x(6)*v(4)+x(10)*v(6)
  cheb24(3) = cheb12(3)+alam
  cheb24(23) = cheb12(3)-alam
  alam = x(6)*(v(2)-v(4)-v(6))
  cheb24(7) = cheb12(7)+alam
  cheb24(19) = cheb12(7)-alam
  alam = x(10)*v(2)-x(6)*v(4)+x(2)*v(6)
  cheb24(11) = cheb12(11)+alam
  cheb24(15) = cheb12(11)-alam

  do i = 1, 3
    j = 8-i
    v(i) = fval(i)-fval(j)
    fval(i) = fval(i)+fval(j)
  end do

  cheb12(5) = v(1)+x(8)*v(3)
  cheb12(9) = fval(1)-x(8)*fval(3)
  alam = x(4)*v(2)
  cheb24(5) = cheb12(5)+alam
  cheb24(21) = cheb12(5)-alam
  alam = x(8)*fval(2)-fval(4)
  cheb24(9) = cheb12(9)+alam