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

#*****************************************************************************80
#
## A_POWER returns the value of A for an Alpert rule for power 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 12.
#
#    Output, integer A, the value of A.
#
  import numpy as np
  from sys import exit

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

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

  a = a_vec[rule-1]
  
  return a

def j_power ( rule ):

#*****************************************************************************80
#
## J_POWER returns the value of J for an Alpert rule for power 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 12.
#
#    Output, integer J, the value of J.
#
  import numpy as np
  from sys import exit

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

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

  j = j_vec[rule-1]
  
  return j

def num_power ( ):

#*****************************************************************************80
#
## NUM_POWER returns the number of Alpert rules for power 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 = 12

  return num

def order_power ( rule ):

#*****************************************************************************80
#
## ORDER_POWER returns the order of an Alpert rule for power 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, real ORDER, the order of the rule.
#
  import numpy as np
  from sys import exit

  order_vec = np.array ( [ \
     1.5,  2.0,  2.5,  3.0,  3.5, \
     4.0,  6.0,  8.0, 10.0, 12.0, \
    14.0, 16.0 ] )

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

  order = order_vec[rule-1]
  
  return order

def rule_power ( rule, j ):

#*****************************************************************************80
#
## RULE_POWER returns an Alpert rule for power 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 12.
#
#    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 12 < rule ):
    print ( '' )
    print ( 'RULE_POWER - Fatal error!' )
    print ( '  Input value of RULE is not between 1 and 12.' )
    exit ( 'RULE_POWER - Fatal error!' )

  if ( rule == 1 ):

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

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

  elif ( rule == 2 ):

    x = np.array ( [ \
      9.252112715421378E-02, \
      1.000000000000000E-00 ] )

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

  elif ( rule == 3 ):

    x = np.array ( [ \
      6.023873796408450E-02, \
      8.780704050676215E-01 ] )

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

  elif ( rule == 4 ):

    x = np.array ( [ \
      7.262978413470474E-03, \
      2.246325512521893E-01, \
      1.000000000000000E+00 ] )

    w = np.array ( [ \
      3.907638767531813E-02, \
      4.873484056646474E-01, \
      9.735752066600344E-01 ] )

  elif ( rule == 5 ):

    x = np.array ( [ \
      1.282368909458828E-02, \
      2.694286346792474E-01, \
      1.018414523786358E+00 ] )

    w = np.array ( [ \
      6.363996663105925E-02, \
      5.077434578043636E-01, \
      9.286165755645772E-01 ] )

  elif ( rule == 6 ):

    x = np.array ( [ \
      1.189242434021285E-02, \
      2.578220434738662E-01, \
      1.007750064585281E+00, \
      2.000000000000000E+00 ] )

    w = np.array ( [ \
      5.927215035616424E-02, \
      4.955981740306228E-01, \
      9.427131290628058E-01, \
      1.002416546550407E+00 ] )

  elif ( rule == 7 ):

    x = np.array ( [ \
      3.317925942699451E-03, \
      8.283019705296352E-02, \
      4.136094925726231E-01, \
      1.088744373688402E+00, \
      2.006482101852379E+00, \
      3.000000000000000E+00 ] )

    w = np.array ( [ \
      1.681780929883469E-02, \
      1.755244404544475E-01, \
      5.039350503858001E-01, \
      8.266241339680867E-01, \
      9.773065848981277E-01, \
      9.997919809947032E-01 ] )

  elif ( rule == 8 ):

    x = np.array ( [ \
      1.214130606523435E-03, \
      3.223952700027058E-02, \
      1.790935383649920E-01, \
      5.437663805244631E-01, \
      1.176116628396759E+00, \
      2.031848210716014E+00, \
      3.001961225690812E+00, \
      4.000000000000000E+00 ] )

    w = np.array ( [ \
      6.199844884297793E-03, \
      7.106286791720044E-02, \
      2.408930104410471E-01, \
      4.975929263668960E-01, \
      7.592446540441226E-01, \
      9.322446399614420E-01, \
      9.928171438160095E-01, \
      9.999449125689846E-01 ] )

  elif ( rule == 9 ):

    x = np.array ( [ \
      1.745862989163252E-04, \
      8.613670540457314E-03, \
      6.733385088703690E-02, \
      2.514488774733840E-01, \
      6.341845573737690E-01, \
      1.248404055083152E+00, \
      2.065688031953401E+00, \
      3.009199358662542E+00, \
      4.000416269690208E+00, \
      5.000000000000000E+00 ] )

    w = np.array ( [ \
      1.016950985948944E-03, \
      2.294670686517670E-02, \
      1.076657968022888E-01, \
      2.734577662465576E-01, \
      4.978815591924992E-01, \
      7.256208919565360E-01, \
      8.952638690320078E-01, \
      9.778157465381624E-01, \
      9.983390781399277E-01, \
      9.999916342408948E-01 ] )

  elif ( rule == 10 ):

    x = np.array ( [ \
      5.710218427206990E-04, \
      1.540424351115548E-02, \
      8.834248407196555E-02, \
      2.824462054509770E-01, \
      6.574869892305580E-01, \
      1.246541060977993E+00, \
      2.039218495130811E+00, \
      2.979333487049800E+00, \
      3.985772595393049E+00, \
      4.997240804311428E+00, \
      5.999868793951190E+00, \
      7.000000000000000E+00 ] )

    w = np.array ( [ \
      2.921018926912141E-03, \
      3.431130611256885E-02, \
      1.224669495638615E-01, \
      2.761108242022520E-01, \
      4.797809643010337E-01, \
      6.966555677271379E-01, \
      8.790077941972658E-01, \
      9.868622449294327E-01, \
      1.015142389688201E+00, \
      1.006209712632210E+00, \
      1.000528829922287E+00, \
      1.000002397796838E+00 ] )

  elif ( rule == 11 ):

    x = np.array ( [ \
      3.419821460249725E-04, \
      9.296593430187960E-03, \
      5.406214771755252E-02, \
      1.763945096508648E-01, \
      4.218486605653738E-01, \
      8.274022895884040E-01, \
      1.410287585637014E+00, \
      2.160997505238153E+00, \
      3.043504749358223E+00, \
      4.005692579069439E+00, \
      4.999732707905968E+00, \
      5.999875191971098E+00, \
      6.999994560568667E+00, \
      8.000000000000000E+00 ] )

    w = np.array ( [ \
      1.750957243202047E-03, \
      2.080726584287380E-02, \
      7.586830616433430E-02, \
      1.766020526671851E-01, \
      3.206624362072232E-01, \
      4.934405290553812E-01, \
      6.707497030698472E-01, \
      8.244959025366557E-01, \
      9.314646742162802E-01, \
      9.845768443163154E-01, \
      9.992852769154770E-01, \
      1.000273112957723E+00, \
      1.000022857402321E+00, \
      1.000000081405180E+00 ] )

  elif ( rule == 12 ):

    x = np.array ( [ \
      2.158438988280793E-04, \
      5.898432743709196E-03, \
      3.462795956896131E-02, \
      1.145586495070213E-01, \
      2.790344218856415E-01, \
      5.600113798653321E-01, \
      9.814091242883119E-01, \
      1.553594853974655E+00, \
      2.270179114036658E+00, \
      3.108234601715371E+00, \
      4.032930893996553E+00, \
      5.006803270228157E+00, \
      6.000815466735179E+00, \
      7.000045035079542E+00, \
      8.000000738923901E+00, \
      9.000000000000000E+00 ] )

    w = np.array ( [ \
      1.105804873501181E-03, \
      1.324499944707956E-02, \
      4.899842307592144E-02, \
      1.165326192868815E-01, \
      2.178586693194957E-01, \
      3.481766016945031E-01, \
      4.964027915911545E-01, \
      6.469026189623831E-01, \
      7.823688971783889E-01, \
      8.877772445893361E-01, \
      9.551665077035583E-01, \
      9.876285579741800E-01, \
      9.979929183863017E-01, \
      9.998470620634641E-01, \
      9.999962891645340E-01, \
      9.999999946893169E-01 ] )

  return x, w

def alpert_power_test ( ):

#*****************************************************************************80
#
## ALPERT_POWER_TEST tests the Alpert rule on the power 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_power
  from test_functions import integrand_power

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

  v2 = integral_power ( )

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

    a_p = a_power ( rule )
    j_p = j_power ( rule )
    order_p = order_power ( rule )
    x_p, w_p = rule_power ( rule, j_p )

    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, 7 ):

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

      x1 = h * x_p 
      f1 = integrand_power ( j_p, x1 )
      s1 = np.dot ( w_p, f1 )

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

      x3 = 1.0 - h * x_r
      f3 = integrand_power ( 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_p, j_p, a_p, n, n + j_p + j_r, h, v1, abs ( v1 - v2 ) ) )

    print ( '' )

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

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