#! /usr/bin/env python
#
def r8_b1mp ( x ):

#*****************************************************************************80
#
## R8_B1MP evaluates the modulus and phase for the Bessel J1 and Y1 functions.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    26 April 2016
#
#  Author:
#
#    Original FORTRAN77version by Wayne Fullerton.
#    Python version by John Burkardt.
#
#  Reference:
#
#    Wayne Fullerton,
#    Portable Special Function Routines,
#    in Portability of Numerical Software,
#    edited by Wayne Cowell,
#    Lecture Notes in Computer Science, Volume 57,
#    Springer 1977,
#    ISBN: 978-3-540-08446-4,
#    LC: QA297.W65.
#
#  Parameters:
#
#    Input, real X, the argument.
#
#    Output, real AMPL, THETA, the modulus and phase.
#
  import numpy as np
  from r8_csevl import r8_csevl
  from r8_inits import r8_inits
  from machine import r8_mach
  from sys import exit

  bm12cs = np.array ( [ \
      +0.9807979156233050027272093546937E-01, \
      +0.1150961189504685306175483484602E-02, \
      -0.4312482164338205409889358097732E-05, \
      +0.5951839610088816307813029801832E-07, \
      -0.1704844019826909857400701586478E-08, \
      +0.7798265413611109508658173827401E-10, \
      -0.4958986126766415809491754951865E-11, \
      +0.4038432416421141516838202265144E-12, \
      -0.3993046163725175445765483846645E-13, \
      +0.4619886183118966494313342432775E-14, \
      -0.6089208019095383301345472619333E-15, \
      +0.8960930916433876482157048041249E-16, \
      -0.1449629423942023122916518918925E-16, \
      +0.2546463158537776056165149648068E-17, \
      -0.4809472874647836444259263718620E-18, \
      +0.9687684668292599049087275839124E-19, \
      -0.2067213372277966023245038117551E-19, \
      +0.4646651559150384731802767809590E-20, \
      -0.1094966128848334138241351328339E-20, \
      +0.2693892797288682860905707612785E-21, \
      -0.6894992910930374477818970026857E-22, \
      +0.1830268262752062909890668554740E-22, \
      -0.5025064246351916428156113553224E-23, \
      +0.1423545194454806039631693634194E-23, \
      -0.4152191203616450388068886769801E-24, \
      +0.1244609201503979325882330076547E-24, \
      -0.3827336370569304299431918661286E-25, \
      +0.1205591357815617535374723981835E-25, \
      -0.3884536246376488076431859361124E-26, \
      +0.1278689528720409721904895283461E-26, \
      -0.4295146689447946272061936915912E-27, \
      +0.1470689117829070886456802707983E-27, \
      -0.5128315665106073128180374017796E-28, \
      +0.1819509585471169385481437373286E-28, \
      -0.6563031314841980867618635050373E-29, \
      +0.2404898976919960653198914875834E-29, \
      -0.8945966744690612473234958242979E-30, \
      +0.3376085160657231026637148978240E-30, \
      -0.1291791454620656360913099916966E-30, \
      +0.5008634462958810520684951501254E-31 ] )

  bm1cs = np.array ( [ \
      +0.1069845452618063014969985308538, \
      +0.3274915039715964900729055143445E-02, \
      -0.2987783266831698592030445777938E-04, \
      +0.8331237177991974531393222669023E-06, \
      -0.4112665690302007304896381725498E-07, \
      +0.2855344228789215220719757663161E-08, \
      -0.2485408305415623878060026596055E-09, \
      +0.2543393338072582442742484397174E-10, \
      -0.2941045772822967523489750827909E-11, \
      +0.3743392025493903309265056153626E-12, \
      -0.5149118293821167218720548243527E-13, \
      +0.7552535949865143908034040764199E-14, \
      -0.1169409706828846444166290622464E-14, \
      +0.1896562449434791571721824605060E-15, \
      -0.3201955368693286420664775316394E-16, \
      +0.5599548399316204114484169905493E-17, \
      -0.1010215894730432443119390444544E-17, \
      +0.1873844985727562983302042719573E-18, \
      -0.3563537470328580219274301439999E-19, \
      +0.6931283819971238330422763519999E-20, \
      -0.1376059453406500152251408930133E-20, \
      +0.2783430784107080220599779327999E-21, \
      -0.5727595364320561689348669439999E-22, \
      +0.1197361445918892672535756799999E-22, \
      -0.2539928509891871976641440426666E-23, \
      +0.5461378289657295973069619199999E-24, \
      -0.1189211341773320288986289493333E-24, \
      +0.2620150977340081594957824000000E-25, \
      -0.5836810774255685901920938666666E-26, \
      +0.1313743500080595773423615999999E-26, \
      -0.2985814622510380355332778666666E-27, \
      +0.6848390471334604937625599999999E-28, \
      -0.1584401568222476721192960000000E-28, \
      +0.3695641006570938054301013333333E-29, \
      -0.8687115921144668243012266666666E-30, \
      +0.2057080846158763462929066666666E-30, \
      -0.4905225761116225518523733333333E-31 ] )

  bt12cs = np.array ( [ \
      +0.73823860128742974662620839792764, \
      -0.33361113174483906384470147681189E-02, \
      +0.61463454888046964698514899420186E-04, \
      -0.24024585161602374264977635469568E-05, \
      +0.14663555577509746153210591997204E-06, \
      -0.11841917305589180567005147504983E-07, \
      +0.11574198963919197052125466303055E-08, \
      -0.13001161129439187449366007794571E-09, \
      +0.16245391141361731937742166273667E-10, \
      -0.22089636821403188752155441770128E-11, \
      +0.32180304258553177090474358653778E-12, \
      -0.49653147932768480785552021135381E-13, \
      +0.80438900432847825985558882639317E-14, \
      -0.13589121310161291384694712682282E-14, \
      +0.23810504397147214869676529605973E-15, \
      -0.43081466363849106724471241420799E-16, \
      +0.80202544032771002434993512550400E-17, \
      -0.15316310642462311864230027468799E-17, \
      +0.29928606352715568924073040554666E-18, \
      -0.59709964658085443393815636650666E-19, \
      +0.12140289669415185024160852650666E-19, \
      -0.25115114696612948901006977706666E-20, \
      +0.52790567170328744850738380799999E-21, \
      -0.11260509227550498324361161386666E-21, \
      +0.24348277359576326659663462400000E-22, \
      -0.53317261236931800130038442666666E-23, \
      +0.11813615059707121039205990399999E-23, \
      -0.26465368283353523514856789333333E-24, \
      +0.59903394041361503945577813333333E-25, \
      -0.13690854630829503109136383999999E-25, \
      +0.31576790154380228326413653333333E-26, \
      -0.73457915082084356491400533333333E-27, \
      +0.17228081480722747930705920000000E-27, \
      -0.40716907961286507941068800000000E-28, \
      +0.96934745136779622700373333333333E-29, \
      -0.23237636337765716765354666666666E-29, \
      +0.56074510673522029406890666666666E-30, \
      -0.13616465391539005860522666666666E-30, \
      +0.33263109233894654388906666666666E-31 ] )

  bth1cs = np.array ( [ \
      +0.74749957203587276055443483969695, \
      -0.12400777144651711252545777541384E-02, \
      +0.99252442404424527376641497689592E-05, \
      -0.20303690737159711052419375375608E-06, \
      +0.75359617705690885712184017583629E-08, \
      -0.41661612715343550107630023856228E-09, \
      +0.30701618070834890481245102091216E-10, \
      -0.28178499637605213992324008883924E-11, \
      +0.30790696739040295476028146821647E-12, \
      -0.38803300262803434112787347554781E-13, \
      +0.55096039608630904934561726208562E-14, \
      -0.86590060768383779940103398953994E-15, \
      +0.14856049141536749003423689060683E-15, \
      -0.27519529815904085805371212125009E-16, \
      +0.54550796090481089625036223640923E-17, \
      -0.11486534501983642749543631027177E-17, \
      +0.25535213377973900223199052533522E-18, \
      -0.59621490197413450395768287907849E-19, \
      +0.14556622902372718620288302005833E-19, \
      -0.37022185422450538201579776019593E-20, \
      +0.97763074125345357664168434517924E-21, \
      -0.26726821639668488468723775393052E-21, \
      +0.75453300384983271794038190655764E-22, \
      -0.21947899919802744897892383371647E-22, \
      +0.65648394623955262178906999817493E-23, \
      -0.20155604298370207570784076869519E-23, \
      +0.63417768556776143492144667185670E-24, \
      -0.20419277885337895634813769955591E-24, \
      +0.67191464220720567486658980018551E-25, \
      -0.22569079110207573595709003687336E-25, \
      +0.77297719892989706370926959871929E-26, \
      -0.26967444512294640913211424080920E-26, \
      +0.95749344518502698072295521933627E-27, \
      -0.34569168448890113000175680827627E-27, \
      +0.12681234817398436504211986238374E-27, \
      -0.47232536630722639860464993713445E-28, \
      +0.17850008478186376177858619796417E-28, \
      -0.68404361004510395406215223566746E-29, \
      +0.26566028671720419358293422672212E-29, \
      -0.10450402527914452917714161484670E-29, \
      +0.41618290825377144306861917197064E-30, \
      -0.16771639203643714856501347882887E-30, \
      +0.68361997776664389173535928028528E-31, \
      -0.28172247861233641166739574622810E-31 ] )

  eta = 0.1 * r8_mach ( 3 )
  nbm1 = r8_inits ( bm1cs, 37, eta )
  nbt12 = r8_inits ( bt12cs, 39, eta )
  nbm12 = r8_inits ( bm12cs, 40, eta )
  nbth1 = r8_inits ( bth1cs, 44, eta )
  xmax = 1.0 / r8_mach ( 4 )

  if ( x < 4.0 ):
    print ( '' )
    print ( 'R8_B1MP - Fatal error!' )
    print ( '  X < 4.' )
    exit ( 'R8_B1MP - Fatal error!' )
  elif ( x <= 8.0 ):
    z = ( 128.0 / x / x - 5.0 ) / 3.0
    ampl = ( 0.75 + r8_csevl ( z, bm1cs, nbm1 ) ) / np.sqrt ( x )
    theta = x - 0.75 * np.pi + r8_csevl ( z, bt12cs, nbt12 ) / x
  else:
    z = 128.0 / x / x - 1.0
    ampl = ( 0.75 + r8_csevl ( z, bm12cs, nbm12 ) ) / np.sqrt ( x )
    theta = x - 0.75 * np.pi + r8_csevl ( z, bth1cs, nbth1 ) / x

  return ampl, theta

def r8_besj1 ( x ):

#*****************************************************************************80
#
## R8_BESJ1 evaluates the Bessel function J of order 1 of an R8 argument.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    26 April 2016
#
#  Author:
#
#    Original FORTRAN77 version by Wayne Fullerton.
#    Python version by John Burkardt.
#
#  Reference:
#
#    Wayne Fullerton,
#    Portable Special Function Routines,
#    in Portability of Numerical Software,
#    edited by Wayne Cowell,
#    Lecture Notes in Computer Science, Volume 57,
#    Springer 1977,
#    ISBN: 978-3-540-08446-4,
#    LC: QA297.W65.
#
#  Parameters:
#
#    Input, real X, the argument.
#
#    Output, real VALUE, the Bessel function J of order 1 of X.
#
  import numpy as np
  from r8_csevl import r8_csevl
  from r8_inits import r8_inits
  from machine import r8_mach

  bj1cs = np.array ( [ \
      -0.117261415133327865606240574524003, \
      -0.253615218307906395623030884554698, \
      +0.501270809844695685053656363203743E-01, \
      -0.463151480962508191842619728789772E-02, \
      +0.247996229415914024539124064592364E-03, \
      -0.867894868627882584521246435176416E-05, \
      +0.214293917143793691502766250991292E-06, \
      -0.393609307918317979229322764073061E-08, \
      +0.559118231794688004018248059864032E-10, \
      -0.632761640466139302477695274014880E-12, \
      +0.584099161085724700326945563268266E-14, \
      -0.448253381870125819039135059199999E-16, \
      +0.290538449262502466306018688000000E-18, \
      -0.161173219784144165412118186666666E-20, \
      +0.773947881939274637298346666666666E-23, \
      -0.324869378211199841143466666666666E-25, \
      +0.120223767722741022720000000000000E-27, \
      -0.395201221265134933333333333333333E-30, \
      +0.116167808226645333333333333333333E-32 ] )

  ntj1 = r8_inits ( bj1cs, 19, 0.1 * r8_mach ( 3 ) )
  xsml = np.sqrt ( 4.0 * r8_mach ( 3 ) )
  xmin = 2.0 * r8_mach ( 1 )

  y = abs ( x )

  if ( y <= xmin ):
    value = 0.0
  elif ( y <= xsml ):
    value = 0.5 * x
  elif ( y <= 4.0 ):
    value = x * ( 0.25 + r8_csevl ( 0.125 * y * y - 1.0, bj1cs, ntj1 ) )
  else:
    ampl, theta = r8_b1mp ( y )
    if ( x < 0.0 ):
      value = - ampl * np.cos ( theta )
    else:
      value = + ampl * np.cos ( theta )

  return value

def r8_besj1_test ( ):

#*****************************************************************************80
#
## R8_BESJ1_TEST tests R8_BESJ1.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    26 April 2016
#
#  Author:
#
#    John Burkardt
#
  import platform
  from bessel_j1_values import bessel_j1_values

  print ( '' )
  print ( 'R8_BESJ1_TEST:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  R8_BESJ1 evaluates the Bessel J1(x) function' )
  print ( '' )
  print ( '             X       BESJ1(X)  R8_BESJ1(X)         Diff' )
  print ( '' )

  n_data = 0

  while ( True ):

    n_data, x, fx1 = bessel_j1_values ( n_data )

    if ( n_data == 0 ):
      break

    fx2 = r8_besj1 ( x )

    print ( '  %14.4f  %14.6g  %14.6g  %14.6g' % ( x, fx1, fx2, abs ( fx1 - fx2 ) ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'R8_BESJ1_TEST:' )
  print ( '  Normal end of execution.' )
  return

def r8_besy1 ( x ):

#*****************************************************************************80
#
## R8_BESY1 evaluates the Bessel function Y of order 1 of an R8 argument.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    26 April 2016
#
#  Author:
#
#    Original FORTRAN77 version by Wayne Fullerton.
#    Python version by John Burkardt.
#
#  Reference:
#
#    Wayne Fullerton,
#    Portable Special Function Routines,
#    in Portability of Numerical Software,
#    edited by Wayne Cowell,
#    Lecture Notes in Computer Science, Volume 57,
#    Springer 1977,
#    ISBN: 978-3-540-08446-4,
#    LC: QA297.W65.
#
#  Parameters:
#
#    Input, real X, the argument.
#
#    Output, real VALUE, the Bessel function Y of order 1 of X.
#
  import numpy as np
  from r8_csevl import r8_csevl
  from r8_inits import r8_inits
  from machine import r8_mach
  from sys import exit

  twodpi = 0.636619772367581343075535053490057

  by1cs = np.array ( [ \
      +0.320804710061190862932352018628015E-01, \
      +0.126270789743350044953431725999727E+01, \
      +0.649996189992317500097490637314144E-02, \
      -0.893616452886050411653144160009712E-01, \
      +0.132508812217570954512375510370043E-01, \
      -0.897905911964835237753039508298105E-03, \
      +0.364736148795830678242287368165349E-04, \
      -0.100137438166600055549075523845295E-05, \
      +0.199453965739017397031159372421243E-07, \
      -0.302306560180338167284799332520743E-09, \
      +0.360987815694781196116252914242474E-11, \
      -0.348748829728758242414552947409066E-13, \
      +0.278387897155917665813507698517333E-15, \
      -0.186787096861948768766825352533333E-17, \
      +0.106853153391168259757070336000000E-19, \
      -0.527472195668448228943872000000000E-22, \
      +0.227019940315566414370133333333333E-24, \
      -0.859539035394523108693333333333333E-27, \
      +0.288540437983379456000000000000000E-29, \
      -0.864754113893717333333333333333333E-32 ] )

  nty1 = r8_inits ( by1cs, 20, 0.1 * r8_mach ( 3 ) )
  xmin = 1.571 * np.exp ( max ( np.log ( r8_mach ( 1 ) ), \
    - np.log ( r8_mach ( 2 ) ) ) + 0.01 )
  xsml = np.sqrt ( 4.0 * r8_mach ( 3 ) )

  if ( x <= 0.0 ):
    print ( '' )
    print ( 'R8_BESY1 - Fatal error!' )
    print ( '  X <= 0.' )
    error ( 'R8_BESY1 - Fatal error!' )
  elif ( x <= xmin ):
    y = 0.0
    value = twodpi * np.log ( 0.5 * x ) * r8_besj1 ( x ) \
      + ( 0.5 + r8_csevl ( 0.125 * y - 1.0, by1cs, nty1 ) ) / x
  elif ( x <= 4.0 ):
    y = x * x
    value = twodpi * np.log ( 0.5 * x ) * r8_besj1 ( x ) \
      + ( 0.5 + r8_csevl ( 0.125 * y - 1.0, by1cs, nty1 ) ) / x
  else:
    ampl, theta = r8_b1mp ( x )
    value = ampl * np.sin ( theta )

  return value

def r8_besy1_test ( ):

#*****************************************************************************80
#
## R8_BESY1_TEST tests R8_BESY1.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    26 April 2016
#
#  Author:
#
#    John Burkardt
#
  import platform
  from bessel_y1_values import bessel_y1_values

  print ( '' )
  print ( 'R8_BESY1_TEST:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  R8_BESY1 evaluates the Bessel Y1(x) function' )
  print ( '' )
  print ( '             X       BESY1(X)  R8_BESY1(X)         Diff' )
  print ( '' )

  n_data = 0

  while ( True ):

    n_data, x, fx1 = bessel_y1_values ( n_data )

    if ( n_data == 0 ):
      break

    fx2 = r8_besy1 ( x )

    print ( '  %14.4g  %14.6g  %14.6g  %14.6g' % ( x, fx1, fx2, abs ( fx1 - fx2 ) ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'R8_BESY1_TEST:' )
  print ( '  Normal end of execution.' )
  return

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( )
  r8_besj1_test ( )
  r8_besy1_test ( )
  timestamp ( )