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

#*****************************************************************************80
#
## R8_B0MP evaluates the modulus and phase for the Bessel J0 and Y0 functions.
#
#  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 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

  bm0cs = np.array ( [ \
      +0.9211656246827742712573767730182E-01, \
      -0.1050590997271905102480716371755E-02, \
      +0.1470159840768759754056392850952E-04, \
      -0.5058557606038554223347929327702E-06, \
      +0.2787254538632444176630356137881E-07, \
      -0.2062363611780914802618841018973E-08, \
      +0.1870214313138879675138172596261E-09, \
      -0.1969330971135636200241730777825E-10, \
      +0.2325973793999275444012508818052E-11, \
      -0.3009520344938250272851224734482E-12, \
      +0.4194521333850669181471206768646E-13, \
      -0.6219449312188445825973267429564E-14, \
      +0.9718260411336068469601765885269E-15, \
      -0.1588478585701075207366635966937E-15, \
      +0.2700072193671308890086217324458E-16, \
      -0.4750092365234008992477504786773E-17, \
      +0.8615128162604370873191703746560E-18, \
      -0.1605608686956144815745602703359E-18, \
      +0.3066513987314482975188539801599E-19, \
      -0.5987764223193956430696505617066E-20, \
      +0.1192971253748248306489069841066E-20, \
      -0.2420969142044805489484682581333E-21, \
      +0.4996751760510616453371002879999E-22, \
      -0.1047493639351158510095040511999E-22, \
      +0.2227786843797468101048183466666E-23, \
      -0.4801813239398162862370542933333E-24, \
      +0.1047962723470959956476996266666E-24, \
      -0.2313858165678615325101260800000E-25, \
      +0.5164823088462674211635199999999E-26, \
      -0.1164691191850065389525401599999E-26, \
      +0.2651788486043319282958336000000E-27, \
      -0.6092559503825728497691306666666E-28, \
      +0.1411804686144259308038826666666E-28, \
      -0.3298094961231737245750613333333E-29, \
      +0.7763931143074065031714133333333E-30, \
      -0.1841031343661458478421333333333E-30, \
      +0.4395880138594310737100799999999E-31 ] )

  bm02cs = np.array ( [ \
      +0.9500415145228381369330861335560E-01, \
      -0.3801864682365670991748081566851E-03, \
      +0.2258339301031481192951829927224E-05, \
      -0.3895725802372228764730621412605E-07, \
      +0.1246886416512081697930990529725E-08, \
      -0.6065949022102503779803835058387E-10, \
      +0.4008461651421746991015275971045E-11, \
      -0.3350998183398094218467298794574E-12, \
      +0.3377119716517417367063264341996E-13, \
      -0.3964585901635012700569356295823E-14, \
      +0.5286111503883857217387939744735E-15, \
      -0.7852519083450852313654640243493E-16, \
      +0.1280300573386682201011634073449E-16, \
      -0.2263996296391429776287099244884E-17, \
      +0.4300496929656790388646410290477E-18, \
      -0.8705749805132587079747535451455E-19, \
      +0.1865862713962095141181442772050E-19, \
      -0.4210482486093065457345086972301E-20, \
      +0.9956676964228400991581627417842E-21, \
      -0.2457357442805313359605921478547E-21, \
      +0.6307692160762031568087353707059E-22, \
      -0.1678773691440740142693331172388E-22, \
      +0.4620259064673904433770878136087E-23, \
      -0.1311782266860308732237693402496E-23, \
      +0.3834087564116302827747922440276E-24, \
      -0.1151459324077741271072613293576E-24, \
      +0.3547210007523338523076971345213E-25, \
      -0.1119218385815004646264355942176E-25, \
      +0.3611879427629837831698404994257E-26, \
      -0.1190687765913333150092641762463E-26, \
      +0.4005094059403968131802476449536E-27, \
      -0.1373169422452212390595193916017E-27, \
      +0.4794199088742531585996491526437E-28, \
      -0.1702965627624109584006994476452E-28, \
      +0.6149512428936330071503575161324E-29, \
      -0.2255766896581828349944300237242E-29, \
      +0.8399707509294299486061658353200E-30, \
      -0.3172997595562602355567423936152E-30, \
      +0.1215205298881298554583333026514E-30, \
      -0.4715852749754438693013210568045E-31 ] )

  bt02cs = np.array ( [ \
      -0.24548295213424597462050467249324, \
      +0.12544121039084615780785331778299E-02, \
      -0.31253950414871522854973446709571E-04, \
      +0.14709778249940831164453426969314E-05, \
      -0.99543488937950033643468850351158E-07, \
      +0.85493166733203041247578711397751E-08, \
      -0.86989759526554334557985512179192E-09, \
      +0.10052099533559791084540101082153E-09, \
      -0.12828230601708892903483623685544E-10, \
      +0.17731700781805131705655750451023E-11, \
      -0.26174574569485577488636284180925E-12, \
      +0.40828351389972059621966481221103E-13, \
      -0.66751668239742720054606749554261E-14, \
      +0.11365761393071629448392469549951E-14, \
      -0.20051189620647160250559266412117E-15, \
      +0.36497978794766269635720591464106E-16, \
      -0.68309637564582303169355843788800E-17, \
      +0.13107583145670756620057104267946E-17, \
      -0.25723363101850607778757130649599E-18, \
      +0.51521657441863959925267780949333E-19, \
      -0.10513017563758802637940741461333E-19, \
      +0.21820381991194813847301084501333E-20, \
      -0.46004701210362160577225905493333E-21, \
      +0.98407006925466818520953651199999E-22, \
      -0.21334038035728375844735986346666E-22, \
      +0.46831036423973365296066286933333E-23, \
      -0.10400213691985747236513382399999E-23, \
      +0.23349105677301510051777740800000E-24, \
      -0.52956825323318615788049749333333E-25, \
      +0.12126341952959756829196287999999E-25, \
      -0.28018897082289428760275626666666E-26, \
      +0.65292678987012873342593706666666E-27, \
      -0.15337980061873346427835733333333E-27, \
      +0.36305884306364536682359466666666E-28, \
      -0.86560755713629122479172266666666E-29, \
      +0.20779909972536284571238399999999E-29, \
      -0.50211170221417221674325333333333E-30, \
      +0.12208360279441714184191999999999E-30, \
      -0.29860056267039913454250666666666E-31 ] )

  bth0cs = np.array ( [ \
      -0.24901780862128936717709793789967, \
      +0.48550299609623749241048615535485E-03, \
      -0.54511837345017204950656273563505E-05, \
      +0.13558673059405964054377445929903E-06, \
      -0.55691398902227626227583218414920E-08, \
      +0.32609031824994335304004205719468E-09, \
      -0.24918807862461341125237903877993E-10, \
      +0.23449377420882520554352413564891E-11, \
      -0.26096534444310387762177574766136E-12, \
      +0.33353140420097395105869955014923E-13, \
      -0.47890000440572684646750770557409E-14, \
      +0.75956178436192215972642568545248E-15, \
      -0.13131556016891440382773397487633E-15, \
      +0.24483618345240857495426820738355E-16, \
      -0.48805729810618777683256761918331E-17, \
      +0.10327285029786316149223756361204E-17, \
      -0.23057633815057217157004744527025E-18, \
      +0.54044443001892693993017108483765E-19, \
      -0.13240695194366572724155032882385E-19, \
      +0.33780795621371970203424792124722E-20, \
      -0.89457629157111779003026926292299E-21, \
      +0.24519906889219317090899908651405E-21, \
      -0.69388422876866318680139933157657E-22, \
      +0.20228278714890138392946303337791E-22, \
      -0.60628500002335483105794195371764E-23, \
      +0.18649748964037635381823788396270E-23, \
      -0.58783732384849894560245036530867E-24, \
      +0.18958591447999563485531179503513E-24, \
      -0.62481979372258858959291620728565E-25, \
      +0.21017901684551024686638633529074E-25, \
      -0.72084300935209253690813933992446E-26, \
      +0.25181363892474240867156405976746E-26, \
      -0.89518042258785778806143945953643E-27, \
      +0.32357237479762298533256235868587E-27, \
      -0.11883010519855353657047144113796E-27, \
      +0.44306286907358104820579231941731E-28, \
      -0.16761009648834829495792010135681E-28, \
      +0.64292946921207466972532393966088E-29, \
      -0.24992261166978652421207213682763E-29, \
      +0.98399794299521955672828260355318E-30, \
      -0.39220375242408016397989131626158E-30, \
      +0.15818107030056522138590618845692E-30, \
      -0.64525506144890715944344098365426E-31, \
      +0.26611111369199356137177018346367E-31 ] )

  eta = 0.1 * r8_mach ( 3 )
  nbm0 = r8_inits ( bm0cs, 37, eta )
  nbt02 = r8_inits ( bt02cs, 39, eta )
  nbm02 = r8_inits ( bm02cs, 40, eta )
  nbth0 = r8_inits ( bth0cs, 44, eta )
  xmax = 1.0 / r8_mach ( 4 )

  if ( x < 4.0 ):
    print ( '' )
    print ( 'R8_B0MP - Fatal error!' )
    print ( '  X < 4.' )
    exit ( 'R8_B0MP - Fatal error!' )
  elif ( x <= 8.0 ):
    z = ( 128.0 / x / x - 5.0 ) / 3.0
    ampl = ( 0.75 + r8_csevl ( z, bm0cs, nbm0 ) ) / np.sqrt ( x )
    theta = x - 0.25 * np.pi + r8_csevl ( z, bt02cs, nbt02 ) / x
  else:
    z = 128.0 / x / x - 1.0
    ampl = ( 0.75 + r8_csevl ( z, bm02cs, nbm02) ) / np.sqrt ( x )
    theta = x - 0.25 * np.pi + r8_csevl ( z, bth0cs, nbth0 ) / x

  return ampl, theta

def r8_besj0 ( x ):

#*****************************************************************************80
#
## R8_BESJ0 evaluates the Bessel function J of order 0 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 0 of X.
#
  import numpy as np
  from r8_csevl import r8_csevl
  from r8_inits import r8_inits
  from machine import r8_mach

  bj0cs = np.array ( [ \
      +0.10025416196893913701073127264074, \
      -0.66522300776440513177678757831124, \
      +0.24898370349828131370460468726680, \
      -0.33252723170035769653884341503854E-01, \
      +0.23114179304694015462904924117729E-02, \
      -0.99112774199508092339048519336549E-04, \
      +0.28916708643998808884733903747078E-05, \
      -0.61210858663032635057818407481516E-07, \
      +0.98386507938567841324768748636415E-09, \
      -0.12423551597301765145515897006836E-10, \
      +0.12654336302559045797915827210363E-12, \
      -0.10619456495287244546914817512959E-14, \
      +0.74706210758024567437098915584000E-17, \
      -0.44697032274412780547627007999999E-19, \
      +0.23024281584337436200523093333333E-21, \
      -0.10319144794166698148522666666666E-23, \
      +0.40608178274873322700800000000000E-26, \
      -0.14143836005240913919999999999999E-28, \
      +0.43910905496698880000000000000000E-31 ] )

  ntj0 = r8_inits ( bj0cs, 19, 0.1 * r8_mach ( 3 ) )
  xsml = np.sqrt ( 4.0 * r8_mach ( 3 ) )

  y = abs ( x )

  if ( y <= xsml ):
    value = 1.0
  elif ( y <= 4.0 ):
    value = r8_csevl ( 0.125 * y * y - 1.0, bj0cs, ntj0 )
  else:
    ampl, theta = r8_b0mp ( y )
    value = ampl * np.cos ( theta )

  return value

def r8_besj0_test ( ):

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

  print ( '' )
  print ( 'R8_BESJ0_TEST:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  R8_BESJ0 evaluates the Bessel J0(x) function' )
  print ( '' )
  print ( '             X       BESJ0(X)  R8_BESJ0(X)         Diff' )
  print ( '' )

  n_data = 0

  while ( True ):

    n_data, x, fx1 = bessel_j0_values ( n_data )

    if ( n_data == 0 ):
      break

    fx2 = r8_besj0 ( x )

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

def r8_besy0 ( x ):

#*****************************************************************************80
#
## R8_BESY0 evaluates the Bessel function Y of order 0 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 0 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

  alnhaf = -0.69314718055994530941723212145818
  twodpi = 0.636619772367581343075535053490057

  by0cs = np.array ( [ \
      -0.1127783939286557321793980546028E-01, \
      -0.1283452375604203460480884531838, \
      -0.1043788479979424936581762276618, \
      +0.2366274918396969540924159264613E-01, \
      -0.2090391647700486239196223950342E-02, \
      +0.1039754539390572520999246576381E-03, \
      -0.3369747162423972096718775345037E-05, \
      +0.7729384267670667158521367216371E-07, \
      -0.1324976772664259591443476068964E-08, \
      +0.1764823261540452792100389363158E-10, \
      -0.1881055071580196200602823012069E-12, \
      +0.1641865485366149502792237185749E-14, \
      -0.1195659438604606085745991006720E-16, \
      +0.7377296297440185842494112426666E-19, \
      -0.3906843476710437330740906666666E-21, \
      +0.1795503664436157949829120000000E-23, \
      -0.7229627125448010478933333333333E-26, \
      +0.2571727931635168597333333333333E-28, \
      -0.8141268814163694933333333333333E-31 ] )

  nty0 = r8_inits ( by0cs, 19, 0.1 * r8_mach ( 3 ) )
  xsml = np.sqrt ( 4.0 * r8_mach ( 3 ) )

  if ( x <= 0.0 ):
    print ( '' )
    print ( 'R8_BESY0 - Fatal error!' )
    print ( '  X <= 0.' )
    exit ( 'R8_BESY0 - Fatal error!' )
  elif ( x <= xsml ):
    y = 0.0
    value = twodpi * ( alnhaf + np.log ( x ) ) * r8_besj0 ( x ) \
      + 0.375 + r8_csevl ( 0.125 * y - 1.0, by0cs, nty0 )
  elif ( x <= 4.0 ):
    y = x * x
    value = twodpi * ( alnhaf + np.log ( x ) ) * r8_besj0 ( x ) \
      + 0.375 + r8_csevl ( 0.125 * y - 1.0, by0cs, nty0 )
  else:
    ampl, theta = r8_b0mp ( x )
    value = ampl * np.sin ( theta )

  return value

def r8_besy0_test ( ):

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

  print ( '' )
  print ( 'R8_BESY0_TEST:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  R8_BESY0 evaluates the Bessel Y0(X) function.' )
  print ( '' )
  print ( '             X       BESY0(X)  R8_BESY0(X)         Diff' )
  print ( '' )

  n_data = 0

  while ( True ):

    n_data, x, fx1 = bessel_y0_values ( n_data )

    if ( n_data == 0 ):
      break

    fx2 = r8_besy0 ( x )

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

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( )
  r8_besj0_test ( )
  r8_besy0_test ( )
  timestamp ( )