Skip to content

Commit bef8cb8

Browse files
heltonmcheltonmc
committed
add complex besseli0 and besselj0
1 parent 4c4b013 commit bef8cb8

File tree

1 file changed

+145
-0
lines changed

1 file changed

+145
-0
lines changed

‎src/i0.jl

Lines changed: 145 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,145 @@
1+
"""
2+
besseli0(x::T) where T <: Union{Float32, Float64}
3+
4+
Modified Bessel function of the first kind of order zero, ``I_0(x)``.
5+
"""
6+
function besseli0(z::Complex{T}) where T <: Union{Float64}
7+
check = false
8+
if real(z) < 0.0
9+
z = -z
10+
end
11+
if imag(z) < 0.0
12+
z = conj(z)
13+
check = true
14+
end
15+
16+
if abs(z) > 20.0
17+
zinv = 1 / z
18+
e = exp(z)
19+
s = sqrt(2 * z * π)
20+
p = evalpoly(zinv*zinv, (1.0, 0.0703125, 0.112152099609375, 0.5725014209747314, 6.074042001273483, 110.01714026924674, 3038.090510922384, 118838.42625678325, 6.252951493434797e6, 4.259392165047669e8, 3.646840080706556e10, 3.8335346613939443e12))
21+
p2 = zinv*evalpoly(zinv*zinv, (0.125, 0.0732421875, 0.22710800170898438, 1.7277275025844574, 24.380529699556064, 551.3358961220206, 18257.755474293175, 832859.3040162893, 5.0069589531988926e7, 3.8362551802304335e9, 3.6490108188498334e11))
22+
r = e / s * (p + p2) + im * (p - p2) / (e * s)
23+
elseif abs(z) < 5.0
24+
r = evalpoly(z*z, (1.0, 0.25, 0.015625, 0.00043402777777777775, 6.781684027777777e-6, 6.781684027777778e-8, 4.709502797067901e-10, 2.4028075495244395e-12, 9.385966990329842e-15, 2.896903392077112e-17, 7.242258480192779e-20, 1.4963343967340453e-22, 2.5978027721077174e-25, 3.842903509035085e-28, 4.9016626390753635e-31, 5.4462918211948485e-34, 5.318644356635594e-37))
25+
else
26+
if angle(z) <= pi/5
27+
#r = besseli_power_series(0.0, z)
28+
r = evalpoly(z*z, (1.0, 0.25, 0.015625, 0.00043402777777777775, 6.781684027777777e-6, 6.781684027777778e-8, 4.709502797067901e-10, 2.4028075495244395e-12, 9.385966990329842e-15, 2.896903392077112e-17, 7.242258480192779e-20, 1.4963343967340453e-22, 2.5978027721077174e-25, 3.842903509035085e-28, 4.9016626390753635e-31, 5.4462918211948485e-34, 5.318644356635594e-37, 4.60090342269515e-40, 3.5500798014623073e-43, 2.458504017633177e-46, 1.5365650110207356e-49, 8.71068600351891e-53, 4.499321282809354e-56, 2.1263333094562166e-59, 9.228877211181495e-63, 3.691550884472598e-66, 1.3652185223641266e-69, 4.681819349671216e-73, 1.4929270885431173e-76, 4.4379521062518346e-80, 1.232764473958843e-83, 3.206983543077115e-87, 7.829549665715613e-91, 1.7974172786307652e-94, 3.887148093924665e-98, 7.932955293723807e-102))
29+
else
30+
zz = [
31+
-1.678231235746891147897486007423140108585 + 20.29678265316490737291132973041385412216im,
32+
-1.523025526043160660094599734293296933174 + 18.74439925781758020661982300225645303726im,
33+
-1.602232319163288343233375599083956331015 + 19.57994203993233739424795203376561403275im,
34+
-1.297548314533615521071396869956515729427 + 15.67962995660580283185936423251405358315im,
35+
-1.083284426743612138821504231600556522608 + 13.2089741880721351918737127562053501606im,
36+
-0.8756326508559164611966707525425590574741 + 10.59162940924948692611451406264677643776im,
37+
-0.6632096279193985255417942425992805510759 + 8.104093334608059251422673696652054786682im,
38+
-0.420162057582130654687091464438708499074 + 5.169659732326902457089090603403747081757im,
39+
-0.5248407961081728023700065932644065469503 + 6.352576446680355815033180988393723964691im,
40+
16.81912545211429588221108133438974618912 + 11.61870568962804384227638365700840950012im,
41+
3.956891894434294787430417272844351828098 + 2.733954688042517844337453425396233797073im,
42+
0.1577160301669917608080595528008416295052 + 20.49513202616177309778322523925453424454im,
43+
-0.390590921015831871176970935266581363976 + 4.825951828615894889651372068328782916069im,
44+
5.086291973430167701053505879826843738556 + 19.856103426521318766617696383036673069im,
45+
-1.393155126117511022343364857078995555639 + 16.86197099455173997739620972424745559692im,
46+
12.05655555056696748295053112087771296501 + 16.57204880191947538037311460357159376144im,
47+
2.460344943558517627479886868968605995178 + 4.13151392233745973214809055207297205925im,
48+
-1.582576281040411148026691989798564463854 + 20.43576085524474805765748897101730108261im,
49+
8.550688367994505156843842996750026941299 + 5.903012749617017718151146254967898130417im,
50+
-0.9675432524019281776972434272465761750937 + 11.81279027041785312235333549324423074722im,
51+
4.703843537806423391600674221990630030632 + 3.252480256589709295411694256472401320934im,
52+
16.06349030987121295765973627567291259766 + 12.69219331974861475487159623298794031143im,
53+
1.937834775724029512389279261697083711624 + 20.39705296333508854900173901114612817764im,
54+
-1.594679144973751627745173209405038505793 + 19.93238790107540125973173417150974273682im
55+
]
56+
57+
58+
w = [ -0.006546682597891755918395606528292773873545 + 0.0im,
59+
-0.008478994639378135619867116190562228439376 - 0.006865842132444678960756512253738037543371im,
60+
-0.04084417105131761538405754663472180254757 - 0.06528666329071823593022116938300314359367im,
61+
-0.04712461267953402949126839871496486011893 + 0.01452078721700496635738097950252267764881im,
62+
-0.007118804602080300823752079253381452872418 + 0.06492808208351363852273152588168159127235im,
63+
-0.06116738548641747347245356536404869984835 - 0.02725053730346357894198661142581840977073im,
64+
-0.02494387696913849192248413544348295545205 + 0.02277714123083315542195315117623977130279im,
65+
-0.04245963986977926291066509634219983126968 + 0.004116793595813594310028893374919789494015im,
66+
-0.02406420353735235287406801774068298982456 + 0.02715441509989760873744479852121003204957im,
67+
-0.4904880414270129662668296077754348516464 + 0.1081159025127441941638295475058839656413im,
68+
-0.03348344660847365344968906697431521024555 + 0.05785977178570968909587790562909503933042im,
69+
0.003536996149578506415389611561295168939978 + 0.07752104476498027085806796776523697189987im,
70+
0.01586745997824464543546341133151145186275 - 0.01721172123513132687366855577693058876321im,
71+
0.1214650011644559657320030510163633152843 + 0.0199244295229646650735588764291605912149im,
72+
-0.01413089963393250950152157940920005785301 + 0.03915878901517700488854600848753761965781im,
73+
0.1940642102318209105682456083741271868348 + 0.189125716371402241566812563178245909512im,
74+
-0.0629596993796838894086320692622393835336 + 0.03943410994216315496041502797197608742863im,
75+
0.02664896602239549827650932911637937650084 + 0.007480078058428093340515019349368230905384im,
76+
-0.01891162412416614799215430764434131560847 - 0.224558251516702916950052326683362480253im,
77+
-0.07674844370101535639960843582230154424906 + 0.06062885086104936177564539434570178855211im,
78+
-0.05865831566992357748446806908759754151106 - 0.1487332818764729724936302091009565629065im,
79+
0.5977968397140664968958390090847387909889 - 0.3921463933084315400812158713961252942681im,
80+
0.09468179551658635617616965873821754939854 + 0.06801353116428925094094637415764736942947im,
81+
-0.03593292611744278858276757659950817469507 + 0.08129780087330089333175209276305395178497im]
82+
83+
84+
f = [ 3.267431673967855942919413791969418525696 - 11.08592268019023485692287067649886012077im,
85+
-1.408157475121772916892837201885413378477 + 8.194246359208529284501310030464082956314im,
86+
10.16651031219827316931514360476285219193 + 1.139462837357397617665810685139149427414im,
87+
0.05333588002795774246633797588401648681611 + 5.332849960825335244862799299880862236023im,
88+
3.731512177678290687055095986579544842243 + 1.009153918438497665732711539021693170071im,
89+
2.079126843703456906098381296033039689064 - 1.574054723602581251640231130295433104038im,
90+
-0.3027225201326076420293986757314996793866 - 1.336196815913121005436892119178082793951im,
91+
-0.321207149871566732812766531424131244421 - 0.5904585480867294844387060948065482079983im,
92+
0.5332409552090798809942384650639723986387 + 1.12410228905475118033052694954676553607im,
93+
0.4009721205547844280481228906864998862147 - 0.001452802358581500212844628272534919233294im,
94+
0.4076248719737363135351415621698834002018 - 0.007221143389127381538583616560345035395585im,
95+
0.3572930776018199416910192667273804545403 - 0.2904023883125960159290457340830471366644im,
96+
0.2227466624325098176750969969361904077232 - 0.86310886714297674338070009980583563447im,
97+
0.3994990829482118477322671878937399014831 - 0.002392663621286134686266811044674795994069im,
98+
5.223429579451738469231258932268247008324 - 4.319291202018760600367386359721422195435im,
99+
0.4003498255225810820157050784473540261388 - 0.002032103233714888101957285471144132316113im,
100+
0.4060493902079563288687324984493898227811 - 0.01098448458664161297981820553104626014829im,
101+
0.164972906678356817655739519068447407335 - 9.454457188929774602570432762149721384048im,
102+
0.4029739937090559553922730628983117640018 - 0.002997378292125214404445499027929145086091im,
103+
-2.361384545672124968263005939661525189877 + 0.1425123166686921016843569987031514756382im,
104+
0.406335503077652593351132281895843334496 - 0.005903178413394780457701394027481001103297im,
105+
0.4008673043752389864025076349207665771246 - 0.001578422845174482381375158368541633535642im,
106+
0.3995411619962026539276678249734686687589 - 0.01068906314642052712837738681628252379596im,
107+
8.455251786581097661610328941605985164642 - 5.380704447124407430180781375383958220482im
108+
]
109+
110+
111+
C = 1 ./ (z .- transpose(zz))
112+
r = ((C*(w.*f))./(C*w))
113+
114+
r = r[1] / (sqrt(z) * exp(-z))
115+
end
116+
end
117+
check && (r = conj(r))
118+
return r
119+
end
120+
121+
_besselj0(z::Complex{T}) where T <: Float64 = besseli0(z/im)
122+
123+
#=
124+
# works for x > 20.0
125+
function i0_large_argument(z::T) where T
126+
c = exp(z) / sqrt(2 * pi * z)
127+
p = evalpoly(1/z, (1.0, 1/8, 9/128, 75/1024, 3675/32768, 59535/262144, 2401245/4194304, 57972915/33554432, 13043905875/2147483648, 418854310875/17179869184, 30241281245175/274877906944, 1212400457192925/2199023255552, 213786613951685775/70368744177664, 10278202593831046875/562949953421312, 1070401384414690453125/9007199254740992, 60013837619516978071875/72057594037927936, 57673297952355815927071875/9223372036854775808, 3694483615889146090857721875/73786976294838206464, 502860269940467106811189921875/1180591620717411303424, 36232405765710498380237842265625/9444732965739290427392, 11021897833929133607268351617203125/302231454903657293676544, 882276678992136837800861860405640625/2417851639229258349412352, T(148302689041496455735799416353639046875)/T(38685626227668133590597632)))
128+
return c*p
129+
end
130+
function i1_large_argument(z::T) where T
131+
c = exp(z) / sqrt(2 * pi * z)
132+
p = evalpoly(1/z, (1.0, -3/8, -15/128, -105/1024, -4725/32768, -72765/262144, -2837835/4194304, -66891825/33554432, -14783093325/2147483648, -468131288625/17179869184, -33424574007825/274877906944, -1327867167401775/2199023255552, -232376754295310625/70368744177664, -11100458801337530625/562949953421312, -1149690375852815671875/9007199254740992, -64152722972587114490625/72057594037927936, -61394155884765868567528125/9223372036854775808, -3918391713821821611515765625/73786976294838206464, -531595142508493798628972203125/1180591620717411303424, -38190914185478633427818266171875/9444732965739290427392, -11587123363874217382000061956546875/302231454903657293676544, -925314565772241073791147804815671875/2417851639229258349412352, -155200488531798616467697063625901328125/38685626227668133590597632))
133+
return c*p
134+
end
135+
136+
function i0_small_argument(z::T) where T
137+
return evalpoly(z*z, (1.0, 1/4, 1/64, 1/2304, 1/147456, 1/14745600, 1/2123366400, 1/416179814400, 1/106542032486400, 1/34519618525593600, 1/13807847410237440000, 1/6682998146554920960000, 1/3849406932415634472960000, 1/2602199086312968903720960000, 1/2040124083669367620517232640000, 1/1836111675302430858465509376000000, 1/1880178355509689199068681601024000000, 1/2173486178969200714123395930783744000000, 1/2816838087944084125503921126295732224000000, 1/4067514198991257477227662106371037331456000000, 1/6508022718386011963564259370193659730329600000000, 1/11480152075232925103727353529021615764301414400000000, 1/22225574417650943000816156432185848119687538278400000000, 1/47029315467749395389726987010505254621258830997094400000000, 1/108355542837694606977930978072204106647380346617305497600000000, 1/270888857094236517444827445180510266618450866543263744000000000000, 1/732483469582815543170813411768099760936291143132985163776000000000000, 1/2135921797303490123886091908715778902890224973375784737570816000000000000, 1/6698250756343745028506784225732682639463745516506460937022078976000000000000, 1/22532915544340358275896822135364744399156039917527734592142273675264000000000000, 1/81118495959625289793228559687313079836961743703099844531712185230950400000000000000, 1/311819498468799613965170583438031478893280942794715802379901640027773337600000000000000, 1/1277212665728203218801338709762176937546878741687155926548077117553759590809600000000000000, 1/5563538371912053221098631419724042739954203798789251216043423924064176777566617600000000000000, 1/25725801431721334094360071684803973629548238365601497622984792224872753419468039782400000000000000, 1/126056427015434537062364351255539470784786367991447338352625481901876491755393394933760000000000000000))
138+
end
139+
140+
function k0_large_argument(z::T) where T
141+
c = exp(-z) * sqrt(pi / (2*z))
142+
p = evalpoly(1/z, (1.0, -1/8, 9/128, -75/1024, 3675/32768, -59535/262144, 2401245/4194304, -57972915/33554432, 13043905875/2147483648, -418854310875/17179869184, 30241281245175/274877906944, -1212400457192925/2199023255552, 213786613951685775/70368744177664, -10278202593831046875/562949953421312, 1070401384414690453125/9007199254740992, -60013837619516978071875/72057594037927936, 57673297952355815927071875/9223372036854775808, -3694483615889146090857721875/73786976294838206464, 502860269940467106811189921875/1180591620717411303424, -36232405765710498380237842265625/9444732965739290427392, 11021897833929133607268351617203125/302231454903657293676544, -882276678992136837800861860405640625/2417851639229258349412352, T(148302689041496455735799416353639046875)/T(38685626227668133590597632)))
143+
return c*p
144+
end
145+
=#

0 commit comments

Comments
 (0)