add hankel, i1, k1

Author: heltonmc

src/hankel.jl

@@ -0,0 +1,12 @@
+function besselh(nu::Float64, k::Integer, x::AbstractFloat)
+    if k == 1
+        return complex(besselj(nu, x), bessely(nu, x))
+    elseif k == 2
+        return complex(besselj(nu, x), -bessely(nu, x))
+    else
+        throw(ArgumentError("k must be 1 or 2"))
+    end
+end
+
+hankelh1(nu, z) = besselh(nu, 1, z)
+hankelh2(nu, z) = besselh(nu, 2, z)

src/i1.jl

@@ -0,0 +1,86 @@
+const A_i1 = (
+    2.77791411276104639959E-18,
+-2.11142121435816608115E-17,
+ 1.55363195773620046921E-16,
+-1.10559694773538630805E-15,
+ 7.60068429473540693410E-15,
+-5.04218550472791168711E-14,
+ 3.22379336594557470981E-13,
+-1.98397439776494371520E-12,
+ 1.17361862988909016308E-11,
+-6.66348972350202774223E-11,
+ 3.62559028155211703701E-10,
+-1.88724975172282928790E-9,
+ 9.38153738649577178388E-9,
+-4.44505912879632808065E-8,
+ 2.00329475355213526229E-7,
+-8.56872026469545474066E-7,
+ 3.47025130813767847674E-6,
+-1.32731636560394358279E-5,
+ 4.78156510755005422638E-5,
+-1.61760815825896745588E-4,
+ 5.12285956168575772895E-4,
+-1.51357245063125314899E-3,
+ 4.15642294431288815669E-3,
+-1.05640848946261981558E-2,
+ 2.47264490306265168283E-2,
+-5.29459812080949914269E-2,
+ 1.02643658689847095384E-1,
+-1.76416518357834055153E-1,
+ 2.52587186443633654823E-1
+)
+
+const B_i1 = (
+    7.51729631084210481353E-18,
+ 4.41434832307170791151E-18,
+-4.65030536848935832153E-17,
+-3.20952592199342395980E-17,
+ 2.96262899764595013876E-16,
+ 3.30820231092092828324E-16,
+-1.88035477551078244854E-15,
+-3.81440307243700780478E-15,
+ 1.04202769841288027642E-14,
+ 4.27244001671195135429E-14,
+-2.10154184277266431302E-14,
+-4.08355111109219731823E-13,
+-7.19855177624590851209E-13,
+ 2.03562854414708950722E-12,
+ 1.41258074366137813316E-11,
+ 3.25260358301548823856E-11,
+-1.89749581235054123450E-11,
+-5.58974346219658380687E-10,
+-3.83538038596423702205E-9,
+-2.63146884688951950684E-8,
+-2.51223623787020892529E-7,
+-3.88256480887769039346E-6,
+-1.10588938762623716291E-4,
+-9.76109749136146840777E-3,
+ 7.78576235018280120474E-1
+)
+
+function besseli1(x)
+    z = abs(x)
+    if z <= 8.0
+        y = z / 2.0 - 2.0
+        z = chbevl(y, A_i1) * z * exp(z)
+    else
+        z = exp(z) * chbevl(32.0 / z - 2.0, B_i1) / sqrt(z)
+    end
+    if x < zero(x)
+        z = -z
+    end
+    return z
+end
+function besseli1(x)
+    z = abs(x)
+    if z <= 8.0
+        y = z / 2.0 - 2.0
+        z = chbevl(y, A_i1) * z
+    else
+        z = chbevl(32.0 / z - 2.0, B_i1) / sqrt(z)
+    end
+    if x < zero(x)
+        z = -z
+    end
+    return z
+end

src/k1.jl

@@ -0,0 +1,69 @@
+function besselk1(x)
+    z = 0.5 * x
+    if x <= zero(x)
+        return throw(DomainError(x, "NaN result for non-NaN input."))
+    end
+    if x <= 2.0
+        y = x * x - 2.0
+        y = log(z) * besseli1(x) + chbevl(y, A_k1) / x
+        return y
+    else
+        return exp(-x) * chbevl(8.0 / x - 2.0, B_k1) / sqrt(x)
+    end
+end
+function besselk1e(x)
+    z = 0.5 * x
+    if x <= zero(x)
+        return throw(DomainError(x, "NaN result for non-NaN input."))
+    end
+    if x <= 2.0
+        y = x * x - 2.0
+        y = log(z) * besseli1(x) + chbevl(y, A_k1) / x
+        return y * exp(x)
+    else
+        return chbevl(8.0 / x - 2.0, B_k1) / sqrt(x)
+    end
+end
+
+
+const A_k1 = (
+    -7.02386347938628759343E-18,
+-2.42744985051936593393E-15,
+-6.66690169419932900609E-13,
+-1.41148839263352776110E-10,
+-2.21338763073472585583E-8,
+-2.43340614156596823496E-6,
+-1.73028895751305206302E-4,
+-6.97572385963986435018E-3,
+-1.22611180822657148235E-1,
+-3.53155960776544875667E-1,
+ 1.52530022733894777053E0
+)
+
+const B_k1 = (
+    -5.75674448366501715755E-18,
+ 1.79405087314755922667E-17,
+-5.68946255844285935196E-17,
+ 1.83809354436663880070E-16,
+-6.05704724837331885336E-16,
+ 2.03870316562433424052E-15,
+-7.01983709041831346144E-15,
+ 2.47715442448130437068E-14,
+-8.97670518232499435011E-14,
+ 3.34841966607842919884E-13,
+-1.28917396095102890680E-12,
+ 5.13963967348173025100E-12,
+-2.12996783842756842877E-11,
+ 9.21831518760500529508E-11,
+-4.19035475934189648750E-10,
+ 2.01504975519703286596E-9,
+-1.03457624656780970260E-8,
+ 5.74108412545004946722E-8,
+-3.50196060308781257119E-7,
+ 2.40648494783721712015E-6,
+-1.93619797416608296024E-5,
+ 1.95215518471351631108E-4,
+-2.85781685962277938680E-3,
+ 1.03923736576817238437E-1,
+ 2.72062619048444266945E0
+)