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 {
    BesselResult -> Complex Double
_result :: Complex Double
  , BesselResult -> (Double, Double)
_errors :: (Double, Double)
  , BesselResult -> (Int, Int)
_codes  :: (Int, Int)
} deriving Int -> BesselResult -> ShowS
[BesselResult] -> ShowS
BesselResult -> String
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [BesselResult] -> ShowS
$cshowList :: [BesselResult] -> ShowS
show :: BesselResult -> String
$cshow :: BesselResult -> String
showsPrec :: Int -> BesselResult -> ShowS
$cshowsPrec :: Int -> BesselResult -> ShowS
Show


cpxdbl2cpxcdbl :: Complex Double -> Complex CDouble
cpxdbl2cpxcdbl :: Complex Double -> Complex CDouble
cpxdbl2cpxcdbl Complex Double
z = forall a b. (Real a, Fractional b) => a -> b
realToFrac (forall a. Complex a -> a
realPart Complex Double
z) forall a. a -> a -> Complex a
:+ forall a b. (Real a, Fractional b) => a -> b
realToFrac (forall a. Complex a -> a
imagPart Complex Double
z)


dbl2cdbl :: Double -> CDouble
dbl2cdbl :: Double -> CDouble
dbl2cdbl = forall a b. (Real a, Fractional b) => a -> b
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 :: Int -> Complex Double -> Double -> Int -> IO BesselResult
besselJn Int
n Complex Double
z Double
err Int
subdiv = do
  let z' :: Complex CDouble
z' = Complex Double -> Complex CDouble
cpxdbl2cpxcdbl Complex Double
z
      n' :: CDouble
n' = Double -> CDouble
dbl2cdbl forall a b. (a -> b) -> a -> b
$ forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
      a :: CDouble
a = forall a. Complex a -> a
realPart Complex CDouble
z'
      b :: CDouble
b = forall a. Complex a -> a
imagPart Complex CDouble
z'
  IntegralResult
re <- (CDouble -> CDouble)
-> Double -> Double -> Double -> Double -> Int -> IO IntegralResult
integration 
    (\CDouble
t -> (forall a. Floating a => a -> a
cos (CDouble
a forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sin CDouble
t forall a. Num a => a -> a -> a
+ CDouble
n' forall a. Num a => a -> a -> a
* CDouble
t) forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cosh (CDouble
b forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sin CDouble
t)) forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a
pi) 
      Double
0 forall a. Floating a => a
pi Double
0.0 Double
err Int
subdiv
  IntegralResult
im <- (CDouble -> CDouble)
-> Double -> Double -> Double -> Double -> Int -> IO IntegralResult
integration 
    (\CDouble
t -> -(forall a. Floating a => a -> a
sin (CDouble
a forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sin CDouble
t forall a. Num a => a -> a -> a
+ CDouble
n' forall a. Num a => a -> a -> a
* CDouble
t) forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sinh (CDouble
b forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sin CDouble
t)) forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a
pi) 
      Double
0 forall a. Floating a => a
pi Double
0.0 Double
err Int
subdiv
  forall (m :: * -> *) a. Monad m => a -> m a
return BesselResult {
      _result :: Complex Double
_result = IntegralResult -> Double
_value IntegralResult
re forall a. a -> a -> Complex a
:+ IntegralResult -> Double
_value IntegralResult
im
    , _errors :: (Double, Double)
_errors = (IntegralResult -> Double
_error IntegralResult
re, IntegralResult -> Double
_error IntegralResult
im)
    , _codes :: (Int, Int)
_codes  = (IntegralResult -> Int
_code IntegralResult
re, IntegralResult -> Int
_code IntegralResult
im)
  }


-- Re(t^z)

realPartTpowz :: CDouble -> Complex CDouble -> CDouble
realPartTpowz :: CDouble -> Complex CDouble -> CDouble
realPartTpowz CDouble
t Complex CDouble
z =
  let x :: CDouble
x = forall a. Complex a -> a
realPart Complex CDouble
z
      y :: CDouble
y = forall a. Complex a -> a
imagPart Complex CDouble
z
  in
    CDouble
tforall a. Floating a => a -> a -> a
**CDouble
x forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos (CDouble
y forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
log CDouble
t)

-- Im(t^z)

imagPartTpowz :: CDouble -> Complex CDouble -> CDouble
imagPartTpowz :: CDouble -> Complex CDouble -> CDouble
imagPartTpowz CDouble
t Complex CDouble
z =
  let x :: CDouble
x = forall a. Complex a -> a
realPart Complex CDouble
z
      y :: CDouble
y = forall a. Complex a -> a
imagPart Complex CDouble
z
  in
    CDouble
tforall a. Floating a => a -> a -> a
**CDouble
x forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sin (CDouble
y forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
log CDouble
t)

-- Re(cos(z * cos(t)))

reCosZcosT :: CDouble -> Complex CDouble -> CDouble
reCosZcosT :: CDouble -> Complex CDouble -> CDouble
reCosZcosT CDouble
t Complex CDouble
z =
  let x :: CDouble
x = forall a. Complex a -> a
realPart Complex CDouble
z
      y :: CDouble
y = forall a. Complex a -> a
imagPart Complex CDouble
z
  in
    forall a. Floating a => a -> a
cos (CDouble
x forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos CDouble
t) forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cosh (CDouble
y forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos CDouble
t)

-- Im(cos(z * cos(t)))

imCosZcosT :: CDouble -> Complex CDouble -> CDouble
imCosZcosT :: CDouble -> Complex CDouble -> CDouble
imCosZcosT CDouble
t Complex CDouble
z =
  let x :: CDouble
x = forall a. Complex a -> a
realPart Complex CDouble
z
      y :: CDouble
y = forall a. Complex a -> a
imagPart Complex CDouble
z
  in -forall a. Floating a => a -> a
sin (CDouble
x forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos CDouble
t) forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sinh (CDouble
y forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos CDouble
t)


-- | Bessel-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.

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 :: Complex Double
-> Complex Double -> Double -> Int -> IO BesselResult
besselJnu Complex Double
nu Complex Double
z Double
err Int
subdiv = do
  let z' :: Complex CDouble
z' = Complex Double -> Complex CDouble
cpxdbl2cpxcdbl Complex Double
z
      nu' :: Complex CDouble
nu' = Complex Double -> Complex CDouble
cpxdbl2cpxcdbl Complex Double
nu
      reintegrand :: CDouble -> CDouble
reintegrand CDouble
t = CDouble -> Complex CDouble -> CDouble
reCosZcosT CDouble
t Complex CDouble
z' forall a. Num a => a -> a -> a
* CDouble -> Complex CDouble -> CDouble
realPartTpowz (forall a. Floating a => a -> a
sin CDouble
t) (Complex CDouble
2forall a. Num a => a -> a -> a
*Complex CDouble
nu')
        forall a. Num a => a -> a -> a
- CDouble -> Complex CDouble -> CDouble
imCosZcosT CDouble
t Complex CDouble
z' forall a. Num a => a -> a -> a
* CDouble -> Complex CDouble -> CDouble
imagPartTpowz (forall a. Floating a => a -> a
sin CDouble
t) (Complex CDouble
2forall a. Num a => a -> a -> a
*Complex CDouble
nu')
      imintegrand :: CDouble -> CDouble
imintegrand CDouble
t = CDouble -> Complex CDouble -> CDouble
reCosZcosT CDouble
t Complex CDouble
z' forall a. Num a => a -> a -> a
* CDouble -> Complex CDouble -> CDouble
imagPartTpowz (forall a. Floating a => a -> a
sin CDouble
t) (Complex CDouble
2forall a. Num a => a -> a -> a
*Complex CDouble
nu')
        forall a. Num a => a -> a -> a
+ CDouble -> Complex CDouble -> CDouble
imCosZcosT CDouble
t Complex CDouble
z' forall a. Num a => a -> a -> a
* CDouble -> Complex CDouble -> CDouble
realPartTpowz (forall a. Floating a => a -> a
sin CDouble
t) (Complex CDouble
2forall a. Num a => a -> a -> a
*Complex CDouble
nu')
      cst :: Complex Double
cst = (Complex Double
zforall a. Fractional a => a -> a -> a
/Complex Double
2)forall a. Floating a => a -> a -> a
**Complex Double
nu forall a. Fractional a => a -> a -> a
/ (forall a. Floating a => a -> a
sqrt forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* forall a. Gamma a => a -> a
gamma (Complex Double
nu forall a. Num a => a -> a -> a
+ Complex Double
0.5))
  IntegralResult
re <- (CDouble -> CDouble)
-> Double -> Double -> Double -> Double -> Int -> IO IntegralResult
integration CDouble -> CDouble
reintegrand Double
0 forall a. Floating a => a
pi Double
0.0 Double
err Int
subdiv
  IntegralResult
im <- (CDouble -> CDouble)
-> Double -> Double -> Double -> Double -> Int -> IO IntegralResult
integration CDouble -> CDouble
imintegrand Double
0 forall a. Floating a => a
pi Double
0.0 Double
err Int
subdiv
  forall (m :: * -> *) a. Monad m => a -> m a
return BesselResult {
      _result :: Complex Double
_result = Complex Double
cst forall a. Num a => a -> a -> a
* (IntegralResult -> Double
_value IntegralResult
re forall a. a -> a -> Complex a
:+ IntegralResult -> Double
_value IntegralResult
im)
    , _errors :: (Double, Double)
_errors = (IntegralResult -> Double
_error IntegralResult
re, IntegralResult -> Double
_error IntegralResult
im)
    , _codes :: (Int, Int)
_codes  = (IntegralResult -> Int
_code IntegralResult
re, IntegralResult -> Int
_code IntegralResult
im)
  }


isInteger :: Complex Double -> Bool
isInteger :: Complex Double -> Bool
isInteger Complex Double
z = Double
y forall a. Eq a => a -> a -> Bool
== Double
0 Bool -> Bool -> Bool
&& Double
x forall a. Eq a => a -> a -> Bool
== forall a b. (Integral a, Num b) => a -> b
fromIntegral (forall a b. (RealFrac a, Integral b) => a -> b
floor Double
x :: Int)
  where
    x :: Double
x = forall a. Complex a -> a
realPart Complex Double
z
    y :: Double
y = forall a. Complex a -> a
imagPart Complex Double
z

asInteger :: Complex Double -> Int
asInteger :: Complex Double -> Int
asInteger Complex Double
z = forall a b. (RealFrac a, Integral b) => a -> b
floor (forall a. Complex a -> a
realPart Complex Double
z) :: Int

-- | Bessel-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.

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 :: Complex Double
-> Complex Double -> Double -> Int -> IO BesselResult
besselJ Complex Double
nu Complex Double
z Double
err Int
subdiv
  | Complex Double -> Bool
isInteger Complex Double
nu = Int -> Complex Double -> Double -> Int -> IO BesselResult
besselJn (Complex Double -> Int
asInteger Complex Double
nu) Complex Double
z Double
err Int
subdiv
  | forall a. Complex a -> a
realPart Complex Double
nu forall a. Ord a => a -> a -> Bool
> -Double
0.5 = Complex Double
-> Complex Double -> Double -> Int -> IO BesselResult
besselJnu Complex Double
nu Complex Double
z Double
err Int
subdiv
  | Bool
otherwise = forall a. HasCallStack => String -> a
error String
"Invalid value of the order."