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


-- | 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 :: Complex Double -> Complex Double -> Double -> Int -> IO AngerResult
angerJ Complex Double
nu Complex Double
z Double
err Int
subdiv = do
  let z' :: Complex CDouble
z' = Complex Double -> Complex CDouble
cpxdbl2cpxcdbl Complex Double
z
      x' :: CDouble
x' = forall a. Complex a -> a
realPart Complex CDouble
z'
      y' :: CDouble
y' = forall a. Complex a -> a
imagPart Complex CDouble
z'
      nu' :: Complex CDouble
nu' = Complex Double -> Complex CDouble
cpxdbl2cpxcdbl Complex Double
nu
      a' :: CDouble
a' = forall a. Complex a -> a
realPart Complex CDouble
nu'
      b' :: CDouble
b' = forall a. Complex a -> a
imagPart Complex CDouble
nu'
      reintegrand :: CDouble -> CDouble
reintegrand CDouble
t = forall a. Floating a => a -> a
cosh(CDouble
b' forall a. Num a => a -> a -> a
* CDouble
t forall a. Num a => a -> a -> a
- CDouble
y' forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sin CDouble
t) forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
cos(CDouble
a' forall a. Num a => a -> a -> a
* CDouble
t forall a. Num a => a -> a -> a
- CDouble
x' forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sin CDouble
t)
      imintegrand :: CDouble -> CDouble
imintegrand CDouble
t = -forall a. Floating a => a -> a
sinh(CDouble
b' forall a. Num a => a -> a -> a
* CDouble
t forall a. Num a => a -> a -> a
- CDouble
y' forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sin CDouble
t) forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sin(CDouble
a' forall a. Num a => a -> a -> a
* CDouble
t forall a. Num a => a -> a -> a
- CDouble
x' forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sin CDouble
t)
  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 AngerResult {
      _result :: Complex Double
_result = (IntegralResult -> Double
_value IntegralResult
re forall a. a -> a -> Complex a
:+ IntegralResult -> Double
_value IntegralResult
im) forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a
pi
    , _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)
  }