#! /usr/bin/env python
#
def a_log ( rule ):

#*****************************************************************************80
#
## A_LOG returns the value of A for an Alpert rule for log singular functions.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    06 December 2015
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer RULE, the index of the rule, between 1 and 10.
#
#    Output, integer A, the value of A.
#
  import numpy as np
  from sys import exit

  a_vec = np.array ( [ \
    1, 2, 2, 3, 3, \
    5, 6, 7, 9, 10 ] )

  if ( rule < 1 or 10 < rule ):
    print ( '' )
    print ( 'A_LOG - Fatal error!' )
    print ( '  Input value of RULE is not between 1 and 10.' )
    exit ( 'A_LOG - Fatal error!' )

  a = a_vec[rule-1]
  
  return a

def j_log ( rule ):

#*****************************************************************************80
#
## J_LOG returns the value of J for an Alpert rule for log singular functions.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    06 December 2015
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer RULE, the index of the rule, between 1 and 10.
#
#    Output, integer J, the value of J.
#
  import numpy as np
  from sys import exit

  j_vec = np.array ( [ \
     1,  2,  3,  4,  5, \
     7, 10, 11, 14, 15 ] )

  if ( rule < 1 or 10 < rule ):
    print ( '' )
    print ( 'J_LOG - Fatal error!' )
    print ( '  Input value of RULE is not between 1 and 10.' )
    exit ( 'J_LOG - Fatal error!' )

  j = j_vec[rule-1]
  
  return j

def num_log ( ):

#*****************************************************************************80
#
## NUM_LOG returns the number of Alpert rules for log singular functions.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    06 December 2015
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Output, integer NUM, the number of rules.
#
  num = 10

  return num

def order_log ( rule ):

#*****************************************************************************80
#
## ORDER_LOG returns the order of an Alpert rule for log singular functions.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    06 December 2015
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer RULE, the index of the rule, between 1 and 10.
#
#    Output, integer ORDER, the order of the rule.
#
  import numpy as np
  from sys import exit

  order_vec = np.array ( [ \
     2,  3,  4,  5,  6, \
     8, 10, 12, 14, 16 ] )

  if ( rule < 1 or 10 < rule ):
    print ( '' )
    print ( 'ORDER_LOG - Fatal error!' )
    print ( '  Input value of RULE is not between 1 and 10.' )
    exit ( 'ORDER_LOG - Fatal error!' )

  order = order_vec[rule-1]
  
  return order

def rule_log ( rule, j ):

#*****************************************************************************80
#
## RULE_LOG returns an Alpert rule for log singular functions.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    06 December 2015
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    Bradley Alpert,
#    Hybrid Gauss-Trapezoidal Quadrature Rules,
#    SIAM Journal on Scientific Computing,
#    Volume 20, Number 5, pages 1551-1584, 1999.
#
#  Parameters:
#
#    Input, integer RULE, the index of the rule, between 1 and 10.
#
#    Input, integer J, the number of points in the rule.
#
#    Output, real X[J], W[J], the points and weights needed for the rule.
#
  import numpy as np
  from sys import exit

  if ( rule < 1 or 10 < rule ):
    print ( '' )
    print ( 'RULE_LOG - Fatal error!' )
    print ( '  Input value of RULE is not between 1 and 10.' )
    exit ( 'RULE_LOG - Fatal error!' )

  if ( rule == 1 ):

    x = np.array ( [ \
      1.591549430918953E-01 ] )

    w = np.array ( [ \
      5.0E-01 ] )

  elif ( rule == 2 ):

    x = np.array ( [ \
      1.150395811972836E-01, \
      9.365464527949632E-01 ] )

    w = np.array ( [ \
      3.913373788753340E-01, \
      1.108662621124666E+00 ] )

  elif ( rule == 3 ):

    x = np.array ( [ \
      2.379647284118974E-02, \
      2.935370741501914E-01, \
      1.023715124251890E+00 ] )

    w = np.array ( [ \
      8.795942675593887E-02, \
      4.989017152913699E-01, \
      9.131388579526912E-01 ] )

  elif ( rule == 4 ):

    x = np.array ( [ \
      2.339013027203800E-02, \
      2.854764931311984E-01, \
      1.005403327220700E+00, \
      1.994970303994294E+00 ] )

    w = np.array ( [ \
      8.609736556158105E-02, \
      4.847019685417959E-01, \
      9.152988869123725E-01, \
      1.013901778984250E+00 ] )

  elif ( rule == 5 ):

    x = np.array ( [ \
      4.004884194926570E-03, \
      7.745655373336686E-02, \
      3.972849993523248E-01, \
      1.075673352915104E+00, \
      2.003796927111872E+00 ] )

    w = np.array ( [ \
      1.671879691147102E-02, \
      1.636958371447360E-01, \
      4.981856569770637E-01, \
      8.372266245578912E-01, \
      9.841730844088381E-01 ] )

  elif ( rule == 6 ):

    x = np.array ( [ \
      6.531815708567918E-03, \
      9.086744584657729E-02, \
      3.967966533375878E-01, \
      1.027856640525646E+00, \
      1.945288592909266E+00, \
      2.980147933889640E+00, \
      3.998861349951123E+00 ] )

    w = np.array ( [ \
      2.462194198995203E-02, \
      1.701315866854178E-01, \
      4.609256358650077E-01, \
      7.947291148621894E-01, \
      1.008710414337933E+00, \
      1.036093649726216E+00, \
      1.004787656533285E+00 ] )

  elif ( rule == 7 ):

    x = np.array ( [ \
      1.175089381227308E-03, \
      1.877034129831289E-02, \
      9.686468391426860E-02, \
      3.004818668002884E-01, \
      6.901331557173356E-01, \
      1.293695738083659E+00, \
      2.090187729798780E+00, \
      3.016719313149212E+00, \
      4.001369747872486E+00, \
      5.000025661793423E+00 ] )

    w = np.array ( [ \
      4.560746882084207E-03, \
      3.810606322384757E-02, \
      1.293864997289512E-01, \
      2.884360381408835E-01, \
      4.958111914344961E-01, \
      7.077154600594529E-01, \
      8.741924365285083E-01, \
      9.661361986515218E-01, \
      9.957887866078700E-01, \
      9.998665787423845E-01 ] )

  elif ( rule == 8 ):

    x = np.array ( [ \
      1.674223682668368E-03, \
      2.441110095009738E-02, \
      1.153851297429517E-01, \
      3.345898490480388E-01, \
      7.329740531807683E-01, \
      1.332305048525433E+00, \
      2.114358752325948E+00, \
      3.026084549655318E+00, \
      4.003166301292590E+00, \
      5.000141170055870E+00, \
      6.000001002441859E+00 ] )

    w = np.array ( [ \
      6.364190780720557E-03, \
      4.723964143287529E-02, \
      1.450891158385963E-01, \
      3.021659470785897E-01, \
      4.984270739715340E-01, \
      6.971213795176096E-01, \
      8.577295622757315E-01, \
      9.544136554351155E-01, \
      9.919938052776484E-01, \
      9.994621875822987E-01, \
      9.999934408092805E-01 ] )

  elif ( rule == 9 ):

    x = np.array ( [ \
      9.305182368545380E-04, \
      1.373832458434617E-02, \
      6.630752760779359E-02, \
      1.979971397622003E-01, \
      4.504313503816532E-01, \
      8.571888631101634E-01, \
      1.434505229617112E+00, \
      2.175177834137754E+00, \
      3.047955068386372E+00, \
      4.004974906813428E+00, \
      4.998525901820967E+00, \
      5.999523015116678E+00, \
      6.999963617883990E+00, \
      7.999999488130134E+00 ] )

    w = np.array ( [ \
      3.545060644780164E-03, \
      2.681514031576498E-02, \
      8.504092035093420E-02, \
      1.854526216643691E-01, \
      3.251724374883192E-01, \
      4.911553747260108E-01, \
      6.622933417369036E-01, \
      8.137254578840510E-01, \
      9.235595514944174E-01, \
      9.821609923744658E-01, \
      1.000047394596121E+00, \
      1.000909336693954E+00, \
      1.000119534283784E+00, \
      1.000002835746089E+00 ] )

  elif ( rule == 10 ):

    x = np.array ( [ \
      8.371529832014113E-04, \
      1.239382725542637E-02, \
      6.009290785739468E-02, \
      1.805991249601928E-01, \
      4.142832599028031E-01, \
      7.964747731112430E-01, \
      1.348993882467059E+00, \
      2.073471660264395E+00, \
      2.947904939031494E+00, \
      3.928129252248612E+00, \
      4.957203086563112E+00, \
      5.986360113977494E+00, \
      6.997957704791519E+00, \
      7.999888757524622E+00, \
      8.999998754306120E+00 ] )

    w = np.array ( [ \
      3.190919086626234E-03, \
      2.423621380426338E-02, \
      7.740135521653088E-02, \
      1.704889420286369E-01, \
      3.029123478511309E-01, \
      4.652220834914617E-01, \
      6.401489637096768E-01, \
      8.051212946181061E-01, \
      9.362411945698647E-01, \
      1.014359775369075E+00, \
      1.035167721053657E+00, \
      1.020308624984610E+00, \
      1.004798397441514E+00, \
      1.000395017352309E+00, \
      1.000007149422537E+00 ] )

  return x, w

def alpert_log_test ( ):

#*****************************************************************************80
#
## ALPERT_LOG_TEST tests the Alpert rule on the log integrand.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    06 December 2015
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  from alpert_regular import a_regular
  from alpert_regular import j_regular
  from alpert_regular import order_regular
  from alpert_regular import rule_regular
  from test_functions import integral_log
  from test_functions import integrand_log

  print ( '' )
  print ( 'ALPERT_LOG_TEST' )
  print ( '  Test the Alpert rule on the log integrand.' )
  print ( '' )
  print ( '  Rule    Order   J   A        N     N+2J               H        Estimate           Error' )
  print ( '' )

  v2 = integral_log ( )

  num_l = num_log ( )
#
#  For the righthand interval, use the regular rule of the same index.
#
  for rule in range ( 1, num_l + 1 ):

    a_l = a_log ( rule )
    j_l = j_log ( rule )
    order_l = order_log ( rule )
    x_l, w_l = rule_log ( rule, j_l )

    a_r = a_regular ( rule )
    j_r = j_regular ( rule )
    order_r = order_regular ( rule )
    x_r, w_r = rule_regular ( rule, j_r )

    n = 8

    for nlog in range ( 4, 8 ):

      n = n * 2
      h = 1.0 / float ( n + a_l + a_r - 1 )

      x1 = h * x_l
      f1 = integrand_log ( j_l, x1 )
      s1 = np.dot ( w_l, f1 )

      x2 = np.linspace ( a_l * h, ( a_l + n - 1 ) * h, n )
      f2 = integrand_log ( n, x2 )
      s2 = np.sum ( f2 )

      x3 = 1.0 - h * x_r
      f3 = integrand_log ( j_r, x3 )
      s3 = np.dot ( w_r, f3 )

      v1 = h * ( s1 + s2 + s3 )

      print ( '    %2d     %4.1f  %2d  %2d  %7d  %7d  %14.6g  %14.6g  %14.6g' \
        % ( rule, order_l, j_l, a_l, n, n + j_l + j_r, h, v1, abs ( v1 - v2 ) ) )

    print ( '' )

  print ( '' )
  print ( '                                                    Exact: %14.6g' % ( v2 ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'ALPERT_LOG_TEST' )
  print ( '  Normal end of execution.' )
  return

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