BesselJ 0.1.0.1 → 0.2.0.0

raw patch · 8 files changed

+301/−39 lines, 8 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

• BesselJ.cabal +6/−3
• CHANGELOG.md +6/−0
• README.md +5/−2
• src/Math/AngerJ.hs +47/−0
• src/Math/AngerWeber.hs +51/−0
• src/Math/BesselJ.hs +9/−5
• src/Math/WeberE.hs +47/−0
• tests/Main.hs +130/−29

@@ -1,8 +1,8 @@ cabal-version:       2.2
 name:                BesselJ
-version:             0.1.0.1
-synopsis:            Bessel J-function
-description:         Computation of the Bessel J-function of a complex variable. 
+version:             0.2.0.0
+synopsis:            Bessel J-function, Anger J-function, Weber E-function, and Anger-Weber function.
+description:         Computation of Bessel J-function, Anger J-function, Weber E-function, Anger-Weber function, of a complex variable. 
 homepage:            https://github.com/stla/BesselJ#readme
 license:             BSD-3-Clause
 license-file:        LICENSE
@@ -17,6 +17,9 @@ library
   hs-source-dirs:      src
   exposed-modules:     Math.BesselJ
+                     , Math.AngerJ
+                     , Math.WeberE
+                     , Math.AngerWeber
   build-depends:       base >= 4.7 && < 5
                      , gamma >= 0.10.0.0
                      , numerical-integration >= 0.1.2.3

@@ -9,3 +9,9 @@ ## 0.1.0.1 - 2023-09-22  More unit tests.+++## 0.2.0.0 - 2023-09-22++Added implementations of Anger J-function, Weber E-function, and Anger-Weber function.+

@@ -1,9 +1,12 @@ # BesselJ
 
-*Computation of the Bessel J-function of a complex variable.*
+*Computation of the Bessel J-function, Anger J-function, Weber E-function, and Anger-Weber 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.
+complex number with real part larger than -0.5, or any integer. 
+There is no restriction for the Anger J-function and the Weber E-function. 
+For the Anger-Weber function, the order must be a non-integer 
+complex number with real part larger than -0.5.
 
 ___
 

@@ -0,0 +1,47 @@+module Math.AngerJ
+  ( AngerResult(..), angerJ )
+  where
+import Data.Complex          ( imagPart, realPart, Complex(..) )
+import Numerical.Integration ( integration, IntegralResult(..) )
+import Foreign.C             ( CDouble )
+
+
+-- | Data type to store the result of a computation of the Anger 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 AngerResult = AngerResult {
+    _result :: Complex Double
+  , _errors :: (Double, Double)
+  , _codes  :: (Int, Int)
+} deriving Show
+
+
+cpxdbl2cpxcdbl :: Complex Double -> Complex CDouble
+cpxdbl2cpxcdbl z = realToFrac (realPart z) :+ realToFrac (imagPart z)
+
+
+-- | Anger-J function. It is computed with two integrals. The field @_errors@ 
+-- in the result provides the error estimates of the integrals. The field 
+-- @_codes@ provides the codes indicating success (0) or failure of each integral.
+angerJ :: Complex Double  -- ^ order, complex number 
+       -> 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 AngerResult  -- ^ result
+angerJ nu z err subdiv = do
+  let z' = cpxdbl2cpxcdbl z
+      x' = realPart z'
+      y' = imagPart z'
+      nu' = cpxdbl2cpxcdbl nu
+      a' = realPart nu'
+      b' = imagPart nu'
+      reintegrand t = cosh(b' * t - y' * sin t) * cos(a' * t - x' * sin t)
+      imintegrand t = -sinh(b' * t - y' * sin t) * sin(a' * t - x' * sin t)
+  re <- integration reintegrand 0 pi 0.0 err subdiv
+  im <- integration imintegrand 0 pi 0.0 err subdiv
+  return AngerResult {
+      _result = (_value re :+ _value im) / pi
+    , _errors = (_error re, _error im)
+    , _codes  = (_code re, _code im)
+  }

@@ -0,0 +1,51 @@+module Math.AngerWeber
+  ( AngerWeberResult(..), angerWeber )
+  where
+import Data.Complex          ( realPart, imagPart, Complex(..) )
+import Math.BesselJ          ( besselJ, BesselResult(..) )
+import Math.AngerJ           ( angerJ, AngerResult(..))
+
+
+-- | Data type to store the result of a computation of the Anger-Weber function.
+-- It is based on a computation of the Bessel J-function and a computation of 
+-- the Anger J-function. 
+-- The fields are @_result@ for the value, @_besselResult@ for the result of 
+-- the computation of the Bessel J-function, and @_angerResult@ for the result 
+-- of the computation of the Anger J-function.
+data AngerWeberResult = AngerWeberResult {
+    _result       :: Complex Double
+  , _besselResult :: BesselResult
+  , _angerResult  :: AngerResult
+} deriving Show
+
+
+isInteger :: Complex Double -> Bool
+isInteger z = y == 0 && x == fromIntegral (floor x :: Int)
+  where
+    x = realPart z
+    y = imagPart z
+
+
+aResult :: AngerResult -> Complex Double
+aResult (AngerResult r  _  _) = r
+
+bResult :: BesselResult -> Complex Double
+bResult (BesselResult r _ _) = r
+
+
+-- | Anger-Weber function. 
+angerWeber :: Complex Double       -- ^ order, non-integer complex number with real part larger than -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 AngerWeberResult  -- ^ result
+angerWeber nu z err subdiv 
+  | isInteger nu = error "The order `nu` cannot be an integer."
+  | realPart nu <= -0.5 = error "The real part of the order `nu` must be larger than -0.5."
+  | otherwise = do
+      bessel <- besselJ nu z err subdiv
+      anger <- angerJ nu z err subdiv
+      let bresult = bResult bessel
+          aresult = aResult anger 
+          result = (aresult - bresult) / sin(pi*nu)
+      return (AngerWeberResult result bessel anger)

@@ -6,6 +6,7 @@ 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 
@@ -16,12 +17,15 @@   , _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.
@@ -81,8 +85,8 @@ 
 
 -- | 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.
+-- in the result provides the error estimates of the integrals. The field 
+-- @_codes@ provides the codes 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
@@ -115,8 +119,8 @@ 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.
+-- in the result provides the error estimates of the integrals. The field 
+-- @_codes@ provides the codes 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
@@ -125,4 +129,4 @@ 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."+  | otherwise = error "Invalid value of the order."

@@ -0,0 +1,47 @@+module Math.WeberE
+  ( WeberResult(..), weberE )
+  where
+import Data.Complex          ( imagPart, realPart, Complex(..) )
+import Numerical.Integration ( integration, IntegralResult(..) )
+import Foreign.C             ( CDouble )
+
+
+-- | Data type to store the result of a computation of the Weber E-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 WeberResult = WeberResult {
+    _result :: Complex Double
+  , _errors :: (Double, Double)
+  , _codes  :: (Int, Int)
+} deriving Show
+
+
+cpxdbl2cpxcdbl :: Complex Double -> Complex CDouble
+cpxdbl2cpxcdbl z = realToFrac (realPart z) :+ realToFrac (imagPart z)
+
+
+-- | Weber-E function. It is computed with two integrals. The field @_errors@ 
+-- in the result provides the error estimates of the integrals. The field 
+-- @_codes@ provides the codes indicating success (0) or failure of each integral.
+weberE :: Complex Double  -- ^ order, complex number 
+       -> 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 WeberResult  -- ^ result
+weberE nu z err subdiv = do
+  let z' = cpxdbl2cpxcdbl z
+      x' = realPart z'
+      y' = imagPart z'
+      nu' = cpxdbl2cpxcdbl nu
+      a' = realPart nu'
+      b' = imagPart nu'
+      reintegrand t = cosh(b' * t - y' * sin t) * sin(a' * t - x' * sin t)
+      imintegrand t = sinh(b' * t - y' * sin t) * cos(a' * t - x' * sin t)
+  re <- integration reintegrand 0 pi 0.0 err subdiv
+  im <- integration imintegrand 0 pi 0.0 err subdiv
+  return WeberResult {
+      _result = (_value re :+ _value im) / pi
+    , _errors = (_error re, _error im)
+    , _codes  = (_code re, _code im)
+  }

@@ -3,18 +3,30 @@ import           Data.Complex     ( Complex(..), conjugate )
 import           Test.Tasty       ( defaultMain, testGroup )
 import           Test.Tasty.HUnit ( testCase )
-import           Math.BesselJ     ( BesselResult(_result), besselJ )
+import           Math.BesselJ     ( BesselResult(..), besselJ )
+import           Math.AngerJ      ( AngerResult(..), angerJ )
+import           Math.WeberE      ( WeberResult(..), weberE )
+import           Math.AngerWeber  ( AngerWeberResult(..), angerWeber )
 
-i_ :: Complex Double
-i_ = 0.0 :+ 1.0
 
+aResult :: AngerResult -> Complex Double
+aResult (AngerResult r  _  _) = r
 
+bResult :: BesselResult -> Complex Double
+bResult (BesselResult r _ _) = r
+
+wResult :: WeberResult -> Complex Double
+wResult (WeberResult r  _  _) = r
+
+awResult :: AngerWeberResult -> Complex Double
+awResult (AngerWeberResult r  _  _) = r
+
 main :: IO ()
 main = defaultMain $
   testGroup "Tests"
   [ 
     testCase "nu = 1+2i -- z = 3+4i" $ do
-      my <- _result <$> besselJ (1 :+ 2) (3 :+ 4) 1e-5 5000
+      my <- bResult <$> besselJ (1 :+ 2) (3 :+ 4) 1e-5 5000
       let wolfram = 0.31925 :+ (-0.66956) 
       assertAreClose "" 1e-5 my wolfram,
 
@@ -22,60 +34,149 @@     --   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 
+    --   x1 <- bResult <$> jnu (nu-1) z 
+    --   x2 <- bResult <$> jnu (-nu) z
+    --   x3 <- bResult <$> jnu (1-nu) z
+    --   x4 <- bResult <$> 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
+      x1 <- bResult <$> besselJ nu z 1e-8 10000
+      x2 <- bResult <$> besselJ (nu+1) z 1e-8 10000
+      x3 <- bResult <$> besselJ (nu+2) z 1e-8 10000
       let y = 2*(nu+1)/z * x2 - x3
-      assertAreClose "" 1e-7 x1 y,
+      assertAreClose "" 1e-8 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,
+      x <- bResult <$> besselJ 0.5 z 1e-5 5000
+      assertAreClose "" 1e-10 x s,
 
     testCase "remove square root" $ do
       let z = 2 :+ 1
           nu = (-0.3) :+ 1
-      x <- _result <$> besselJ nu (sqrt (z*z)) 1e-5 5000
-      y <- _result <$> besselJ nu z 1e-5 5000
-      assertAreClose "" 1e-6 x (z**(-nu) * (z*z)**(nu/2) * y),
+      x <- bResult <$> besselJ nu (sqrt (z*z)) 1e-5 5000
+      y <- bResult <$> besselJ nu z 1e-5 5000
+      assertAreClose "" 1e-7 x (z**(-nu) * (z*z)**(nu/2) * y),
 
     testCase "remove square root --- integer nu" $ do
       let z = 2 :+ 1
           nu = -4
-      x <- _result <$> besselJ nu (sqrt (z*z)) 1e-5 5000
-      y <- _result <$> besselJ nu z 1e-5 5000
-      assertAreClose "" 1e-6 x (z**(-nu) * (z*z)**(nu/2) * y),
+      x <- bResult <$> besselJ nu (sqrt (z*z)) 1e-5 5000
+      y <- bResult <$> besselJ nu z 1e-5 5000
+      assertAreClose "" 1e-7 x (z**(-nu) * (z*z)**(nu/2) * y),
 
     testCase "remove minus sign" $ do
       let z = 2 :+ 1
           nu = (-0.3) :+ 1
-      x <- _result <$> besselJ nu (-z) 1e-5 5000
-      y <- _result <$> besselJ nu z 1e-5 5000
-      assertAreClose "" 1e-6 x ((-z)**nu * z**(-nu) * y), 
+      x <- bResult <$> besselJ nu (-z) 1e-5 5000
+      y <- bResult <$> besselJ nu z 1e-5 5000
+      assertAreClose "" 1e-7 x ((-z)**nu * z**(-nu) * y), 
 
     testCase "conjugate" $ do
       let z = 2 :+ 5
           nu = 0.3 :+ (-1)
-      x <- _result <$> besselJ (conjugate nu) (conjugate z) 1e-5 5000
-      y <- _result <$> besselJ nu z 1e-5 5000
-      assertAreClose "" 1e-6 x (conjugate y), 
+      x <- bResult <$> besselJ (conjugate nu) (conjugate z) 1e-5 5000
+      y <- bResult <$> besselJ nu z 1e-5 5000
+      assertAreClose "" 1e-7 x (conjugate y), 
 
     testCase "conjugate --- integer nu" $ do
       let z = 2 :+ 5
           nu = 7
-      x <- _result <$> besselJ nu (conjugate z) 1e-5 5000
-      y <- _result <$> besselJ nu z 1e-5 5000
-      assertAreClose "" 1e-6 x (conjugate y)
+      x <- bResult <$> besselJ nu (conjugate z) 1e-5 5000
+      y <- bResult <$> besselJ nu z 1e-5 5000
+      assertAreClose "" 1e-7 x (conjugate y),
+
+
+    -- Anger --------
+
+    testCase "Anger at z = 0" $ do
+      let nu = (-2.3) :+ 1
+          y = sin(pi*nu) / (pi*nu)
+      x <- aResult <$> angerJ nu 0 1e-5 5000
+      assertAreClose "" 1e-8 x y,
+
+    testCase "Anger is Bessel when nu is integer" $ do
+      let z = 2 :+ 6
+          nu = 10
+      x <- aResult <$> angerJ nu z 1e-5 5000
+      y <- bResult <$> besselJ nu z 1e-5 5000
+      assertAreClose "" 1e-8 x y,
+
+    testCase "Anger - remove minus sign" $ do
+      let z = 2 :+ 1
+          nu = (-0.3) :+ 1
+      x <- aResult <$> angerJ nu (-z) 1e-5 5000
+      y <- aResult <$> angerJ (-nu) z 1e-5 5000
+      assertAreClose "" 1e-8 x y,
+
+    testCase "Anger - recurrence relation" $ do
+      let z = (-2) :+ 1
+          nu = (-0.3) :+ 1
+      x1 <- aResult <$> angerJ (nu-1) z 1e-5 5000
+      x2 <- aResult <$> angerJ (nu+1) z 1e-5 5000
+      y <- aResult <$> angerJ nu z 1e-5 5000
+      let x = x1 + x2 
+          y' = 2*nu/z * y - 2*sin(pi*nu)/(pi*z)
+      assertAreClose "" 1e-8 x y',
+
+
+    -- Weber --------
+
+    testCase "Weber at z = 0" $ do
+      let nu = (-2.3) :+ 1
+          y = (1 - cos(pi*nu)) / (pi*nu)
+      x <- wResult <$> weberE nu 0 1e-5 5000
+      assertAreClose "" 1e-7 x y,
+
+    testCase "Weber - remove minus sign" $ do
+      let z = 2 :+ 1
+          nu = (-0.3) :+ 1
+      x <- wResult <$> weberE nu (-z) 1e-5 5000
+      y <- wResult <$> weberE (-nu) z 1e-5 5000
+      assertAreClose "" 1e-7 x (-y),
+
+    testCase "Weber - recurrence relation" $ do
+      let z = (-2) :+ 1
+          nu = (-0.3) :+ 1
+      x1 <- wResult <$> weberE (nu-1) z 1e-5 5000
+      x2 <- wResult <$> weberE (nu+1) z 1e-5 5000
+      y <- wResult <$> weberE nu z 1e-5 5000
+      let x = x1 + x2 
+          y' = 2*nu/z * y - 2*(1-cos(pi*nu))/(pi*z)
+      assertAreClose "" 1e-7 x y',
+    
+
+    -- Relations between Anger and Weber --------
+
+    testCase "Relation 1 between Anger and Weber" $ do
+      let z = (-7) :+ 6
+          nu = (-3.3) :+ 9
+      w1 <- wResult <$> weberE nu z 1e-5 5000
+      w2 <- wResult <$> weberE (-nu) z 1e-5 5000
+      a <- aResult <$> angerJ nu z 1e-5 5000
+      assertAreClose "" 1e-7 (sin(pi*nu)*a) (cos(pi*nu)*w1 - w2),
+
+    testCase "Relation 2 between Anger and Weber" $ do
+      let z = (-7) :+ 6
+          nu = (-3.3) :+ 9
+      a1 <- aResult <$> angerJ nu z 1e-5 5000
+      a2 <- aResult <$> angerJ (-nu) z 1e-5 5000
+      w <- wResult <$> weberE nu z 1e-5 5000
+      assertAreClose "" 1e-7 (sin(pi*nu)*w) (a2 - cos(pi*nu)*a1),
+    
+
+    -- Anger-Weber --------
+
+    testCase "A value of Anger-Weber" $ do
+      let z = 2.5 :+ 0.5
+          nu = 1.0 / 3.0
+          wolfram = 0.102015 :+ (-0.0162118)
+      aw <- awResult <$> angerWeber nu z 1e-6 5000
+      assertAreClose "" 1e-5 aw wolfram
+
 
   ]