#! /usr/bin/env python
#
def cauchy_principal_value ( f, a, b, n ):

#*****************************************************************************80
#
## CAUCHY_PRINCIPAL_VALUE estimates the Cauchy Principal Value of an integral.
#
#  Location:
#
#    http://people.sc.fsu.edu/~jburkardt/f_src/cauchy_principal_value/cauchy_principal_value.f90
#
#  Discussion:
#
#    This function can be used to estimate the Cauchy Principal Value of
#    a singular integral of the form 
#      Integral f(t)/(t-x) dt 
#    over an interval which includes the singularity point t=x.
#
#    Isolate the singularity at x in a symmetric interval of finite size delta:
#
#        CAUCHY_PRINCIPAL_VALUE ( Integral ( a         <= t <= b         ) p(t) / ( t - x ) dt )
#      =       Integral ( a         <= t <= x - delta ) p(t) / ( t - x ) dt
#      + CAUCHY_PRINCIPAL_VALUE ( Integral ( x - delta <= t <= x + delta ) p(t) / ( t - x ) dt )
#      +       Integral ( x + delta <= t <= b         ) p(t) / ( t - x ) dt.
#
#    We assume the first and third integrals can be handled in the usual way.
#    The second integral can be rewritten as
#      Integral ( -1 <= s <= +1 ) ( p(s*delta+x) - p(x) ) / s ds
#    and approximated by
#      Sum ( 1 <= i <= N ) w(i) * ( p(xi*delta+x) - p(x) ) / xi(i)
#    = Sum ( 1 <= i <= N ) w(i) * ( p(xi*delta+x) ) / xi(i)
#    if we assume that N is even, so that coefficients of p(x) sum to zero.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    01 April 2015
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Julian Noble,
#    Gauss-Legendre Principal Value Integration,
#    Computing in Science and Engineering,
#    Volume 2, Number 1, January-February 2000, pages 92-95.
#
#  Parameters:
#
#    Input, real F ( X ), the function that evaluates the
#    integrand.
#
#    Input, real A, B, the endpoints of the symmetric interval, 
#    which contains a singularity of the form 1/(X-(A+B)/2).
#
#    Input, integer N, the number of Gauss points to use.
#    N must be even.
#
#    Output, real VALUE, the estimate for the Cauchy Principal Value.
#
  from legendre_set import legendre_set
  from sys import exit
#
#  N must be even.
#
  if ( ( n % 2 ) != 0 ):
    print ( '' )
    print ( 'CAUCHY_PRINCIPAL_VALUE - Fatal error!' )
    print ( '  N must be even.' )
    exit ( 'CAUCHY_PRINCIPAL_VALUE - Fatal error.' )
#
#  Get the Gauss-Legendre rule.
#
  [ x, w ] = legendre_set ( n );
#
#  Estimate the integral.
#
  value = 0.0;
  for i in range ( 0, n ):
    x2 = ( ( 1.0 - x[i] ) * a   \
         + ( 1.0 + x[i] ) * b ) \
         /   2.0
    value = value + w[i] * ( f ( x2 ) ) / x[i]

  return value

def cauchy_principal_value_test01 ( ):

#*****************************************************************************80
#
## CAUCHY_PRINCIPAL_VALUE_TEST01 seeks the CAUCHY_PRINCIPAL_VALUE of Integral ( -1 <= t <= 1 ) exp(t) / t dt
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    02 April 2015
#
#  Author:
#
#    John Burkardt
#
  import platform

  print ( '' )
  print ( 'CAUCHY_PRINCIPAL_VALUE_TEST01:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  CAUCHY_PRINCIPAL_VALUE of Integral ( -1 <= t <= 1 ) exp(t) / t dt' )

  print ( '' )
  print ( '   N           Estimate             Error' )
  print ( '' )

  exact = 2.11450175075
  a = -1.0
  b = +1.0
  for n in range ( 2, 10, 2 ):
    value = cauchy_principal_value ( f01, a, b, n )
    print ( '  %2d  %24.16g  %14.6g' % ( n, value, abs ( value - exact ) ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'CAUCHY_PRINCIPAL_VALUE_TEST01:' )
  print ( '  Normal end of execution.' )
  return

def 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.
#
  import numpy as np

  value = np.exp ( t )

  return value

def cauchy_principal_value_test02 ( ):

#*****************************************************************************80
#
## CAUCHY_PRINCIPAL_VALUE_TEST02 is another test.
#
#  Discussion:
#
#    We seek
#      CAUCHY_PRINCIPAL_VALUE ( Integral ( 1-delta <= t <= 1+delta ) 1/(1-t)^3 dt )
#    which we must rewrite as
#      CAUCHY_PRINCIPAL_VALUE ( 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
#
  import numpy as np
  import platform

  print ( '' )
  print ( 'CAUCHY_PRINCIPAL_VALUE_TEST02:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  Compute CAUCHY_PRINCIPAL_VALUE ( Integral ( 1-delta <= t <= 1+delta ) 1/(1-t)^3 dt )' )
  print ( '  Try this for delta = 1, 1/2, 1/4.' )
  print ( '' )
  print ( '   N          Estimate                  Exact                  Error' )

  delta = 1.0

  for k in range ( 0, 3 ):
    print ( '' )
    r1 = (   delta + 1.5 ) ** 2 + 0.75
    r2 = ( - delta + 1.5 ) ** 2 + 0.75
    r3 = np.arctan ( np.sqrt ( 0.75 ) / (   delta + 1.5 ) )
    r4 = np.arctan ( np.sqrt ( 0.75 ) / ( - delta + 1.5 ) )
    exact = - np.log ( r1 / r2 ) / 6.0 + ( r3 - r4 ) / np.sqrt ( 3.0 )
    for n in range ( 2, 10, 2 ):
      a = 1.0 - delta
      b = 1.0 + delta
      value = cauchy_principal_value ( f02, a, b, n )
      print ( '  %2d  %24.16g  %24.16g %14.6g' \
        % ( n, value, exact, abs ( value - exact ) ) )

    delta = delta / 2.0
#
#  Terminate.
#
  print ( '' )
  print ( 'CAUCHY_PRINCIPAL_VALUE_TEST02:' )
  print ( '  Normal end of execution.' )
  return

def 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:
#
  value = 1.0 / ( 1.0 + t + t * t )

  return value

def cauchy_principal_value_test ( ):

#*****************************************************************************80
#
## CAUCHY_PRINCIPAL_VALUE_TEST tests the CAUCHY_PRINCIPAL_VALUE library.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    02 April 2015
#
#  Author:
#
#    John Burkardt
#
  import platform

  from legendre_set  import legendre_set_test01

  print ( '' )
  print ( 'CAUCHY_PRINCIPAL_VALUE_TEST' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  Test the CAUCHY_PRINCIPAL_VALUE library.' )

  cauchy_principal_value_test01 ( )
  cauchy_principal_value_test02 ( )

  legendre_set_test01 ( )
#
#  Terminate.
#
  print ( '' )
  print ( 'CAUCHY_PRINCIPAL_VALUE_TEST:' )
  print ( '  Normal end of execution.' )
  return

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( )
  cauchy_principal_value_test ( )
  timestamp ( )