clean up besselj structure

Author: heltonmc

src/U_polynomials.jl

@@ -43,8 +43,12 @@ function Uk_poly_In(p, v, p2, ::Type{Float32})
     u13 = -1 / 955134445051714068660879360000 * evalpoly(p2, (17438611142828905996129258798828125, -3698121486504259988094897605296209375, 136735019134677724428035696765082813750, -2069933923586966756183324291232117362250, 16843538631795795357786827345307534156063, -83924867223075156862785921508524155665245, 274983827478138958623041508409195988431140, -616410216242554698436702353237166008656700, 962926533925253374359704288384340809260875, -1049095945162229046324321461816274931715625, 782463969315283937856703223178540650343750, -381190503845282445953419057314665534156250, 109361210755577700442544717509565392265625, -14020668045586884672121117629431460546875))
     u14 = 1 / 45846453362482275295722209280000 * evalpoly(p2, (5448320367052402487647812713291015625, -1338184074771428116079233795614103631250, 57170953417612444837142230812990944671875, -1000503839668383458999731491914516564625300, 9440449669103391509091075981237243128469201, -54857705817658080981995319669299096598482382, 211477117385619365164298957821904115470535555, -564850830044980230366140582063618983657685400, 1070439683260179398514645952824098811025619475, -1451823699927947602004297385351260623500372750, 1401302601668131482630653233972052729323190625, -940627071986145750405920450097257805227812500, 417630985812040040477569028872405152769921875, -110320224449895843354117801955418504167031250, 13133360053559028970728309756597440972265625))   
     u15 = -1 / 9820310310243703368343697227776000000 * evalpoly(p2, (8178936810213560828419581728001773291015625, -2303431987527333128955769182299845911305390625, 112597271053778779753048514469995937998172890625, -2254933495791765108580529087615802250458013685625, 24403480234538299231733883413666768614198435948125, -163359140754958502896104062604202934925448173291477, 730367145705123976114617970888707594104468177381925, -2284251621937242886581917667066122422330060024456125, 5136561256208409671660362298619778869859994724706875, -8420533422834140468835467666391400380550043688871875, 10085018700249896522602267572484630409640764997271875, -8735135969643867540297524795790262235822823374296875, 5329871927856528282113994744458999865006055974609375, -2173722139119126857509156976742601742357422740234375, 532039967451707060045861997017872377315039521484375, -59115551939078562227317999668652486368337724609375))
-
-    return evalpoly(-p/v, (u0, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10, u11, u12, u13, u14, u15))
+    u16 = evalpoly(p2, (6.252951493434797e6, -2.0016469281917763e9, 1.1099740513917902e11, -2.5215584749128545e12, 3.100743647289646e13, -2.3665253045164925e14, 1.2126758042503475e15, -4.3793258383640155e15, 1.1486706978449752e16, -2.2268225133911144e16, 3.213827526858624e16, -3.4447226006485144e16, 2.705471130619708e16, -1.5129826322457682e16, 5.705782159023671e15, -1.3010127235496995e15, 1.3552215870309369e14))
+    u17 = -evalpoly(p2, (5.0069589531988926e7, -1.8078220384658062e10, 1.128709145410874e12, -2.886383763141476e13, 4.0004445704303625e14, -3.4503855118462725e15, 2.0064271476309532e16, -8.270945651585064e16, 2.4960365126160426e17, -5.62631788074636e17, 9.575335098169139e17, -1.2336116931960694e18, 1.1961991142756308e18, -8.592577980317548e17, 4.4347954614171904e17, -1.5552983504313904e17, 3.3192764720355224e16, -3.254192619642669e15))
+    u18 = evalpoly(p2, (4.259392165047669e8, -1.722832387173505e11, 1.2030115826419191e13, -3.4396530474307594e14, 5.335106978708839e15, -5.1605093193485224e16, 3.37667624979061e17, -1.5736434765189599e18, 5.402894876715982e18, -1.3970803516443374e19, 2.757282981650519e19, -4.178861444656839e19, 4.859942729324836e19, -4.301555703831444e19, 2.846521225167657e19, -1.3639420410571592e19, 4.47020096401231e18, -8.966114215270463e17, 8.30195760673191e16))
+    u19 = -evalpoly(p2, (3.8362551802304335e9, -1.7277040123529995e12, 1.3412416915180639e14, -4.2619355104268985e15, 7.351663610930971e16, -7.921651119323832e17, 5.789887667664653e18, -3.025566598990372e19, 1.1707490535797259e20, -3.434621399768417e20, 7.756704953461136e20, -1.360203777284994e21, 1.8571089321463453e21, -1.9677247077053125e21, 1.6016898573693598e21, -9.824438427689858e20, 4.392792200888712e20, -1.351217503435996e20, 2.5563802960529236e19, -2.242438856186775e18))
+    u20 = evalpoly(p2, (3.646840080706556e10, -1.818726203851104e13, 1.5613123930484672e15, -5.48403360388329e16, 1.0461721131134344e18, -1.2483700995047234e19, 1.0126774169536592e20, -5.8917941350694964e20, 2.548961114664972e21, -8.405915817108351e21, 2.1487414815055883e22, -4.302534303482379e22, 6.783661642951883e22, -8.423222750084323e22, 8.19433100543513e22, -6.173206302884415e22, 3.528435843903409e22, -1.4787743528433614e22, 4.285296082829494e21, -7.671943936729004e20, 6.393286613940837e19))
+    return evalpoly(-p/v, (u0, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10, u11, u12, u13, u14, u15, u16,  u17, u18, u19, u20))
 end
 
 function Uk_poly_Jn(p, v, p2, ::Type{T}) where T <: Float64
@@ -66,11 +70,12 @@ function Uk_poly_Jn(p, v, p2, ::Type{T}) where T <: Float64
     u15 = -evalpoly(p2, (832859.3040162893, -2.3455796352225152e8, 1.1465754899448236e10, -2.2961937296824646e11, 2.4850009280340854e12, -1.663482472489248e13, 7.437312290867914e13, -2.3260483118893994e14, 5.230548825784446e14, -8.57461032982895e14, 1.0269551960827625e15, -8.894969398810265e14, 5.4273966498765975e14, -2.213496387025252e14, 5.417751075510605e13, -6.019723417234006e12))
     u16 = evalpoly(p2, (6.252951493434797e6, -2.0016469281917763e9, 1.1099740513917902e11, -2.5215584749128545e12, 3.100743647289646e13, -2.3665253045164925e14, 1.2126758042503475e15, -4.3793258383640155e15, 1.1486706978449752e16, -2.2268225133911144e16, 3.213827526858624e16, -3.4447226006485144e16, 2.705471130619708e16, -1.5129826322457682e16, 5.705782159023671e15, -1.3010127235496995e15, 1.3552215870309369e14))
     u17 = -evalpoly(p2, (5.0069589531988926e7, -1.8078220384658062e10, 1.128709145410874e12, -2.886383763141476e13, 4.0004445704303625e14, -3.4503855118462725e15, 2.0064271476309532e16, -8.270945651585064e16, 2.4960365126160426e17, -5.62631788074636e17, 9.575335098169139e17, -1.2336116931960694e18, 1.1961991142756308e18, -8.592577980317548e17, 4.4347954614171904e17, -1.5552983504313904e17, 3.3192764720355224e16, -3.254192619642669e15))
-    
-    return evalpoly(-p/v, (u0, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10, u11, u12, u13, u14, u15, u16,  u17))
+    u18 = evalpoly(p2, (4.259392165047669e8, -1.722832387173505e11, 1.2030115826419191e13, -3.4396530474307594e14, 5.335106978708839e15, -5.1605093193485224e16, 3.37667624979061e17, -1.5736434765189599e18, 5.402894876715982e18, -1.3970803516443374e19, 2.757282981650519e19, -4.178861444656839e19, 4.859942729324836e19, -4.301555703831444e19, 2.846521225167657e19, -1.3639420410571592e19, 4.47020096401231e18, -8.966114215270463e17, 8.30195760673191e16))
+    u19 = -evalpoly(p2, (3.8362551802304335e9, -1.7277040123529995e12, 1.3412416915180639e14, -4.2619355104268985e15, 7.351663610930971e16, -7.921651119323832e17, 5.789887667664653e18, -3.025566598990372e19, 1.1707490535797259e20, -3.434621399768417e20, 7.756704953461136e20, -1.360203777284994e21, 1.8571089321463453e21, -1.9677247077053125e21, 1.6016898573693598e21, -9.824438427689858e20, 4.392792200888712e20, -1.351217503435996e20, 2.5563802960529236e19, -2.242438856186775e18))
+    u20 = evalpoly(p2, (3.646840080706556e10, -1.818726203851104e13, 1.5613123930484672e15, -5.48403360388329e16, 1.0461721131134344e18, -1.2483700995047234e19, 1.0126774169536592e20, -5.8917941350694964e20, 2.548961114664972e21, -8.405915817108351e21, 2.1487414815055883e22, -4.302534303482379e22, 6.783661642951883e22, -8.423222750084323e22, 8.19433100543513e22, -6.173206302884415e22, 3.528435843903409e22, -1.4787743528433614e22, 4.285296082829494e21, -7.671943936729004e20, 6.393286613940837e19))
+    return evalpoly(-p/v, (u0, u1, u2, u3, u4, u5, u6, u7, u8, u9, u10, u11, u12, u13, u14, u15, u16,  u17, u18, u19, u20))
 end
 
-
 #=
 function Uk_poly_Kn(p, v, p2, ::Type{T}) where T <: BigFloat
     u0 = one(T)

src/besselj.jl

@@ -159,55 +159,42 @@ Nu must be real.
 function _besselj(nu, x)
     nu == 0 && return besselj0(x)
     nu == 1 && return besselj1(x)
-    if x < 4.0
-        if nu > 60.0
-            return log_besselj_small_arguments_orders(nu, x)
-        else
-            return besselj_small_arguments_orders(nu, x)
-        end
-    elseif x < 20.0
-        if nu > 45
-            return besselj_debye(nu, x)
-        else
-            v = nu + 45
-            jnu = besselj_debye(v, x)
-            jnup1 = besselj_debye(v+1, x)
-            return besselj_down_recurrence(x, jnu, jnup1, v, nu)[2]
-        end
-    elseif x > 1.65*nu
-        return besselj_large_argument(nu, x)
-    end
-    if nu <= 100
-        if 2.5*x < nu
-            return besselj_debye(nu, x)
-        elseif x <= nu
-            # this can produce inaccurate results when x is close to root of J0 or J1
-            #if isinteger(nu)
-            #    jnu = besselj1(x)
-            #    jnum1 = besselj0(x)
-            #    return besselj_up_recurrence(x, jnu, jnum1, 1.0, nu)[2]
-            #else
-                v = nu + ceil(nu*1.3)
-                jnu = besselj_debye(v, x)
-                jnup1 = besselj_debye(v+1, x)
-                return besselj_down_recurrence(x, jnu, jnup1, v, nu)[2]
-            #end
-        else
-            v = floor(nu - x/3)
-            v2 = nu - v
-            jnu = besselj_large_argument(v2, x)
-            jnum1 = besselj_large_argument(v2 -1, x)
-            return besselj_up_recurrence(x, jnu, jnum1, v2, nu)[2]
-        end
-    elseif 1.5*x < nu
-        return besselj_debye(nu, x)
-    elseif nu < 1000
-        v = nu + ceil(nu*1.3)
+
+    x < 4.0 && return besselj_small_arguments_orders(nu, x)
+
+    large_arg_cutoff = 1.65*nu
+    (x > large_arg_cutoff && x > 20.0) && return besselj_large_argument(nu, x)
+
+
+    debye_cutoff = 2.0 + 1.00035*x + (302.681*x)^(1/3)
+    nu > debye_cutoff && return besselj_debye(nu, x)
+
+    if nu >= x
+        nu_shift = ceil(7.0 + 1.00033*x + (1427.61*x)^(1/3))
+        v = nu + nu_shift
         jnu = besselj_debye(v, x)
         jnup1 = besselj_debye(v+1, x)
         return besselj_down_recurrence(x, jnu, jnup1, v, nu)[2]
+    end
+
+    # at this point x > nu and  x < nu * 1.65
+    # in this region forward recurrence is stable
+    # we must decide if we should do backward recurrence if we are closer to debye accuracy
+    # or if we should do forward recurrence if we are closer to large argument expansion
+    debye_cutoff = 5.0 + 1.00033*x + (1427.61*x)^(1/3)
+
+    debye_diff = debye_cutoff - nu
+    large_arg_diff = nu - x / 2.0
+
+    if (debye_diff > large_arg_diff && x > 20.0)
+        nu_shift = ceil(large_arg_diff)
+        v2 = nu - nu_shift
+        jnu = besselj_large_argument(v2, x)
+        jnum1 = besselj_large_argument(v2 - 1, x)
+        return besselj_up_recurrence(x, jnu, jnum1, v2, nu)[2]
     else
-        v = nu + 500
+        nu_shift = ceil(debye_diff)
+        v = nu + nu_shift
         jnu = besselj_debye(v, x)
         jnup1 = besselj_debye(v+1, x)
         return besselj_down_recurrence(x, jnu, jnup1, v, nu)[2]
@@ -227,8 +214,9 @@ end
 # this needs a better way to sum these as it produces large errors
 # only valid in non-oscillatory regime (v>1/2, 0<t<sqrt(v^2 - 0.25))
 # power series has premature underflow for large orders
-# gives wrong answers for x > 20.0 (might want to fix this)
 function besselj_small_arguments_orders(v, x::T) where T
+    v > 60 && return log_besselj_small_arguments_orders(v, x)
+
     MaxIter = 2000
     out = zero(T)
     a = (x/2)^v / gamma(v + one(T))
@@ -245,7 +233,7 @@ end
 # use when v is large and x is small
 # need for bessely 
 function log_besselj_small_arguments_orders(v, x::T) where T
-    MaxIter = 1000
+    MaxIter = 2000
     out = zero(T)
     a = one(T)
     x2 = (x/2)^2