BesselJ (empty) → 0.1.0.0
raw patch · 8 files changed
+297/−0 lines, 8 filesdep +BesselJdep +basedep +gammasetup-changed
Dependencies added: BesselJ, base, gamma, numerical-integration, system-cxx-std-lib, tasty, tasty-hunit
Files
- BesselJ.cabal +57/−0
- CHANGELOG.md +6/−0
- LICENSE +30/−0
- README.md +10/−0
- Setup.hs +2/−0
- src/Math/BesselJ.hs +128/−0
- tests/Approx.hs +17/−0
- tests/Main.hs +47/−0
+ BesselJ.cabal view
@@ -0,0 +1,57 @@+cabal-version: 2.2 +name: BesselJ +version: 0.1.0.0 +synopsis: Bessel J-function +description: Computation of the Bessel J-function of a complex variable. +homepage: https://github.com/stla/BesselJ#readme +license: BSD-3-Clause +license-file: LICENSE +author: Stéphane Laurent +maintainer: laurent_step@outlook.fr +copyright: 2023 Stéphane Laurent +category: Math +build-type: Simple +extra-source-files: README.md + CHANGELOG.md + +library + hs-source-dirs: src + exposed-modules: Math.BesselJ + build-depends: base >= 4.7 && < 5 + , gamma >= 0.10.0.0 + , numerical-integration >= 0.1.2.3 + if impl(ghc >= 9.4) + build-depends: system-cxx-std-lib == 1.0 + elif os(darwin) || os(freebsd) + extra-libraries: c++11 + else + extra-libraries: stdc++ + default-language: Haskell2010 + ghc-options: -Wall + -Wcompat + -Widentities + -Wincomplete-record-updates + -Wincomplete-uni-patterns + -Wmissing-export-lists + -Wmissing-home-modules + -Wpartial-fields + -Wredundant-constraints + -optcxx-std=c++11 + if os(darwin) || os(freebsd) + ghc-options: -optcxx-stdlib=libc++ + +test-suite unit-tests + type: exitcode-stdio-1.0 + main-is: Main.hs + hs-source-dirs: tests/ + other-modules: Approx + Build-Depends: base >= 4.7 && < 5 + , tasty + , tasty-hunit + , BesselJ + , gamma >= 0.10.0.0 + Default-Language: Haskell2010 + +source-repository head + type: git + location: https://github.com/stla/BesselJ
+ CHANGELOG.md view
@@ -0,0 +1,6 @@+# Changelog for `BesselJ`+++## 0.1.0.0 - 2023-09-21++First release.
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright Stéphane Laurent (c) 2023++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * Redistributions in binary form must reproduce the above+ copyright notice, this list of conditions and the following+ disclaimer in the documentation and/or other materials provided+ with the distribution.++ * Neither the name of Stéphane Laurent nor the names of other+ contributors may be used to endorse or promote products derived+ from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ README.md view
@@ -0,0 +1,10 @@+# BesselJ + +*Computation of the Bessel J-function of a complex variable.* + +The order of the Bessel J-function implemented in this package can be a +complex number with real part larger than -0.5, or any integer. + +___ + +
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ src/Math/BesselJ.hs view
@@ -0,0 +1,128 @@+module Math.BesselJ + ( BesselResult(..), besselJ ) + where +import Data.Complex ( imagPart, realPart, Complex(..) ) +import Numerical.Integration ( integration, IntegralResult(..) ) +import Math.Gamma ( Gamma(gamma) ) +import Foreign.C ( CDouble ) + +-- | Data type to store the result of a computation of the Bessel J-function. +-- The fields are @_result@ for the value, @_errors@ for the error estimates +-- of the integrals used for the computation, and @_codes@ for the convergence +-- codes of these integrals (0 for success). +data BesselResult = BesselResult { + _result :: Complex Double + , _errors :: (Double, Double) + , _codes :: (Int, Int) +} deriving Show + +cpxdbl2cpxcdbl :: Complex Double -> Complex CDouble +cpxdbl2cpxcdbl z = realToFrac (realPart z) :+ realToFrac (imagPart z) + +dbl2cdbl :: Double -> CDouble +dbl2cdbl = realToFrac + +-- | Bessel-J function. It is computed with two integrals. The field @_errors@ +-- in the result are the error estimates of the integrals. The field @_codes@ +-- is the code indicating success (0) or failure. +besselJn :: Int -- ^ order + -> Complex Double -- ^ variable + -> Double -- ^ target relative error accuracy for the integrals + -> Int -- ^ number of subdivisions for the integrals + -> IO BesselResult -- ^ result +besselJn n z err subdiv = do + let z' = cpxdbl2cpxcdbl z + n' = dbl2cdbl $ fromIntegral n + a = realPart z' + b = imagPart z' + re <- integration + (\t -> (cos (a * sin t + n' * t) * cosh (b * sin t)) / pi) + 0 pi 0.0 err subdiv + im <- integration + (\t -> -(sin (a * sin t + n' * t) * sinh (b * sin t)) / pi) + 0 pi 0.0 err subdiv + return BesselResult { + _result = _value re :+ _value im + , _errors = (_error re, _error im) + , _codes = (_code re, _code im) + } + + +-- Re(t^z) +realPartTpowz :: CDouble -> Complex CDouble -> CDouble +realPartTpowz t z = + let x = realPart z + y = imagPart z + in + t**x * cos (y * log t) + +-- Im(t^z) +imagPartTpowz :: CDouble -> Complex CDouble -> CDouble +imagPartTpowz t z = + let x = realPart z + y = imagPart z + in + t**x * sin (y * log t) + +-- Re(cos(z * cos(t))) +reCosZcosT :: CDouble -> Complex CDouble -> CDouble +reCosZcosT t z = + let x = realPart z + y = imagPart z + in + cos (x * cos t) * cosh (y * cos t) + +-- Im(cos(z * cos(t))) +imCosZcosT :: CDouble -> Complex CDouble -> CDouble +imCosZcosT t z = + let x = realPart z + y = imagPart z + in -sin (x * cos t) * sinh (y * cos t) + + +-- | Bessel-J function. It is computed with two integrals. The field @_errors@ +-- in the result are the error estimates of the integrals. The field @_codes@ +-- provides the code indicating success (0) or failure of each integral. +besselJnu :: Complex Double -- ^ order, complex number with real part > -0.5 + -> Complex Double -- ^ the variable, a complex number + -> Double -- ^ target relative accuracy for the integrals, e.g. 1e-5 + -> Int -- ^ number of subdivisions for the integrals, e.g. 5000 + -> IO BesselResult -- ^ result +besselJnu nu z err subdiv = do + let z' = cpxdbl2cpxcdbl z + nu' = cpxdbl2cpxcdbl nu + reintegrand t = reCosZcosT t z' * realPartTpowz (sin t) (2*nu') + - imCosZcosT t z' * imagPartTpowz (sin t) (2*nu') + imintegrand t = reCosZcosT t z' * imagPartTpowz (sin t) (2*nu') + + imCosZcosT t z' * realPartTpowz (sin t) (2*nu') + cst = (z/2)**nu / (sqrt pi * gamma (nu + 0.5)) + re <- integration reintegrand 0 pi 0.0 err subdiv + im <- integration imintegrand 0 pi 0.0 err subdiv + return BesselResult { + _result = cst * (_value re :+ _value im) + , _errors = (_error re, _error im) + , _codes = (_code re, _code im) + } + + +isInteger :: Complex Double -> Bool +isInteger z = y == 0 && x == fromIntegral (floor x :: Int) + where + x = realPart z + y = imagPart z + +asInteger :: Complex Double -> Int +asInteger z = floor (realPart z) :: Int + +-- | Bessel-J function. It is computed with two integrals. The field @_errors@ +-- in the result are the error estimates of the integrals. The field @_codes@ +-- provides the code indicating success (0) or failure of each integral. +besselJ :: Complex Double -- ^ order, integer or complex number with real part > -0.5 + -> Complex Double -- ^ the variable, a complex number + -> Double -- ^ target relative accuracy for the integrals, e.g. 1e-5 + -> Int -- ^ number of subdivisions for the integrals, e.g. 5000 + -> IO BesselResult -- ^ result +besselJ nu z err subdiv + | isInteger nu = besselJn (asInteger nu) z err subdiv + | realPart nu > -0.5 = besselJnu nu z err subdiv + | otherwise = error "Invalid value of the order."
+ tests/Approx.hs view
@@ -0,0 +1,17 @@+module Approx (assertAreClose) where +import Data.Complex ( magnitude, Complex(..) ) +import Test.Tasty.HUnit ( Assertion, assertBool ) + +areClose :: Double -> Complex Double -> Complex Double -> Bool +areClose epsilon z1 z2 = + let maxmod = if magnitude z2 < epsilon + then 1.0 + else max (magnitude z1) (magnitude z2) + in + magnitude (z1 - z2) < 2.0 * epsilon * maxmod + +-- assert approximate equality +assertAreClose :: + String -> Double -> Complex Double -> Complex Double -> Assertion +assertAreClose prefix epsilon z1 z2 = + assertBool prefix (areClose epsilon z1 z2)
+ tests/Main.hs view
@@ -0,0 +1,47 @@+module Main where +import Approx ( assertAreClose ) +import Data.Complex ( Complex(..) ) +import Math.BesselJ +-- import Math.Gamma ( gamma ) +import Test.Tasty ( defaultMain, testGroup ) +import Test.Tasty.HUnit ( testCase ) + +i_ :: Complex Double +i_ = 0.0 :+ 1.0 + + +main :: IO () +main = defaultMain $ + testGroup "Tests" + [ + testCase "nu = 1+2i -- z = 3+4i" $ do + my <- _result <$> besselJ (1 :+ 2) (3 :+ 4) 1e-5 5000 + let wolfram = 0.31925 :+ (-0.66956) + assertAreClose "" 1e-5 my wolfram, + + -- testCase "Relation Bessel-J" $ do + -- let nu = 0.5 :+ 2 + -- z = 3 :+ 4 + -- y = 2 * sin (nu * pi) / (pi * z) + -- x1 <- _result <$> jnu (nu-1) z + -- x2 <- _result <$> jnu (-nu) z + -- x3 <- _result <$> jnu (1-nu) z + -- x4 <- _result <$> jnu nu z + -- assertAreClose "" 1e-3 (x1*x2 + x3*x4) y + + testCase "recurrence relation" $ do + let nu = 0.5 :+ 2 + z = 3 :+ 4 + x1 <- _result <$> besselJ nu z 1e-5 5000 + x2 <- _result <$> besselJ (nu+1) z 1e-5 5000 + x3 <- _result <$> besselJ (nu+2) z 1e-5 5000 + let y = 2*(nu+1)/z * x2 - x3 + assertAreClose "" 1e-7 x1 y, + + testCase "elementary equality" $ do + let z = 3 :+ 4 + s = sqrt(2 / pi / z) * sin z + x <- _result <$> besselJ 0.5 z 1e-5 5000 + assertAreClose "" 1e-9 x s + + ]