program main

!*****************************************************************************80
!
!! CPV_PRB tests the CPV library.
!
!  Location:
!
!    http://people.sc.fsu.edu/~jburkardt/f_src/cauchy_principal_value/cpv_prb.f90
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    01 April 2015
!
!  Author:
!
!    John Burkardt
!
  implicit none

  call timestamp ( )
  write ( *, '(a)' ) ''
  write ( *, '(a)' ) 'CPV_PRB'
  write ( *, '(a)' ) '  FORTRAN90 version'
  write ( *, '(a)' ) '  Test the CPV library.'

  call cpv_test01 ( )
  call cpv_test02 ( )
!
!  Terminate.
!
  write ( *, '(a)' ) ''
  write ( *, '(a)' ) 'CPV_PRB'
  write ( *, '(a)' ) '  Normal end of execution.'
  write ( *, '(a)' ) ''
  call timestamp ( )

  return
end
subroutine cpv_test01 ( )

!*****************************************************************************80
!
!! CPV_TEST01 seeks the CPV of Integral ( -1 <= t <= 1 ) exp(t) / t dt
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    01 April 2015
!
!  Author:
!
!    John Burkardt
!
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) b
  real ( kind = 8 ) cpv
  real ( kind = 8 ) exact
  real ( kind = 8 ), external :: f01
  integer ( kind = 4 ) n
  real ( kind = 8 ) value

  write ( *, '(a)' ) ''
  write ( *, '(a)' ) 'CPV_TEST01:'
  write ( *, '(a)' ) '  CPV of Integral ( -1 <= t <= 1 ) exp(t) / t dt'

  write ( *, '(a)' ) ''
  write ( *, '(a)' ) '   N           Estimate             Error'
  write ( *, '(a)' ) ''

  exact = 2.11450175075D+00
  a = -1.0D+00
  b = +1.0D+00
  do n = 2, 8, 2
    value = cpv ( f01, a, b, n )
    write ( *, '(2x,i2,2x,g24.16,2x,g14.6)' ) n, value, abs ( value - exact )
  end do

  return
end
function f01 ( t )

!*****************************************************************************80
!
!! F01 evaluates the integrand of Integral ( -1 <= t <= 1 ) exp(t) / t dt
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    01 April 2015
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real T, the argument.
!
!    Output, real VALUE, the value of the integrand.
!
  implicit none

  real ( kind = 8 ) f01
  real ( kind = 8 ) t

  f01 = exp ( t )

  return
end
subroutine cpv_test02 ( )

!*****************************************************************************80
!
!! CPV_TEST02 is another test.
!
!  Discussion:
!
!    We seek
!      CPV ( Integral ( 1-delta <= t <= 1+delta ) 1/(1-t)^3 dt )
!    which we must rewrite as
!      CPV ( Integral ( 1-delta <= t <= 1+delta ) 1/(1+t+t^2) 1/(1-t) dt )
!    so that our "integrand" is 1/(1+t+t^2).
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    01 April 2015
!
!  Author:
!
!    John Burkardt
!
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) b
  real ( kind = 8 ) cpv
  real ( kind = 8 ) delta
  real ( kind = 8 ) exact
  real ( kind = 8 ), external :: f02
  integer ( kind = 4 ) k
  integer ( kind = 4 ) n
  real ( kind = 8 ) r1
  real ( kind = 8 ) r2
  real ( kind = 8 ) r3
  real ( kind = 8 ) r4
  real ( kind = 8 ) value

  write ( *, '(a)' ) ''
  write ( *, '(a)' ) 'CPV_TEST02:'
  write ( *, '(a)' ) &
    '  Compute CPV ( Integral ( 1-delta <= t <= 1+delta ) 1/(1-t)^3 dt )'
  write ( *, '(a)' ) '  Try this for delta = 1, 1/2, 1/4.'
  write ( *, '(a)' ) ''
  write ( *, '(a)' ) &
    '   N          Estimate                  Exact                  Error'
  delta = 1.0D+00
  do k = 1, 3
    write ( *, '(a)' ) ''
    r1 = (   delta + 1.5D+00 ) ** 2 + 0.75D+00
    r2 = ( - delta + 1.5D+00 ) ** 2 + 0.75D+00
    r3 = atan ( sqrt ( 0.75D+00 ) / (   delta + 1.5D+00 ) )
    r4 = atan ( sqrt ( 0.75D+00 ) / ( - delta + 1.5D+00 ) )
    exact = - log ( r1 / r2 ) / 6.0D+00 + ( r3 - r4 ) / sqrt ( 3.0D+00 )
    do n = 2, 8, 2
      a = 1.0D+00 - delta
      b = 1.0D+00 + delta
      value = cpv ( f02, a, b, n )
      write ( *, '(2x,i2,g24.16,2x,g24.16,2x,g14.6)' ) &
        n, value, exact, abs ( exact - value )
    end do
    delta = delta / 2.0D+00
  end do

  return
end
function f02 ( t )

!*****************************************************************************80
!
!! F02: integrand of Integral ( 1-delta <= t <= 1+delta ) 1/(1-t^3) dt
!
!  Discussion:
!
!    1/(1-t^3) = 1/(1+t+t^2) * 1/(1-t)
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    01 April 2015
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
  implicit none

  real ( kind = 8 ) f02
  real ( kind = 8 ) t

  f02 = 1.0D+00 / ( 1.0D+00 + t + t * t )

  return
end