add asymptotic phase

Author: heltonmc

src/besselj.jl

@@ -150,3 +150,26 @@ function besselj1(x::Float32)
         return p * s
     end
 end
+
+function besselj_large_argument(v, x::T) where T
+    μ = 4 * v^2
+    s0 = 1
+    s1 = (1 - μ) / 8
+    s2 = (-μ^2 + 26μ - 25) / 128
+    s3 = (-μ^3 + 115μ^2 - 1187μ + 1073) / 1024
+    s4 = (-5μ^4 + 1540μ^3 - 56238μ^2 + 430436μ - 375733) / 32768
+    s5 = (-7μ^5 + 4515μ^4 - 397190μ^3 + 9716998μ^2 - 64709091μ + 55384775) / 262144
+    s6 = (-21μ^6 + 24486μ^5 - 4238531μ^4 + 235370036μ^3 - 4733751627μ^2 + 29215626566μ - 24713030909) / 4194304
+    s7 = (-33μ^7 + 63063μ^6 - 18939349μ^5 + 1989535379μ^4 - 87480924739μ^3 + 1573356635461μ^2 - 9268603618823μ + 7780757249041) / 33554432
+    s8 = (-429μ^8 + 1252680μ^7 - 598859404μ^6 + 106122595896μ^5 - 8593140373614μ^4 + 329343318168440μ^3 - 5517359285625804μ^2 + 31505470994964360μ - 26308967412122125) / 2147483648
+    s9 = (-715μ^9 + 3026595μ^8 - 2163210764μ^7 + 597489288988μ^6 - 79979851361130μ^5 + 5536596631240042μ^4 - 194026764558396188μ^3 + 3095397360215483916μ^2 - 17285630127691126691μ + 14378802319925055947) / 17179869184
+    s10 = (-2431μ^10 + 14318590μ^9 - 14587179035μ^8 + 5925778483368μ^7 - 1216961874423502μ^6 + 137164402798287604μ^5 - 8556293060689145118μ^4 + 281576004385356401192μ^3 - 4334421752432088249971μ^2 + 23801928703130666089534μ -19740662615375374580231) / 274877906944
+
+    αp = s0 + s1/x^2 + s2/x^4 + s3/x^6 + s4/x^8 + s5/x^10 + s6/x^12 + s7/x^14 + s8/x^16 + s9/x^18 + s10/x^20
+    α = s0*x - s1/x - s2/3x^3 - s3/5x^5 - s4/7x^7 - s5/9x^9 - s6/11x^11 - s7/13x^13 - s8/15x^15 - s9/17x^17 - s10/19x^19
+
+    α = α - T(pi)/4 - T(pi)/2 * v
+
+    b = sqrt(2 / T(pi)) / sqrt(αp * x)
+    return cos(α)*b
+end