{-# OPTIONS #-}
-----------------------------------------------------------------------------
{- |
Module      :  Numeric.GSL.Integration
Copyright   :  (c) Alberto Ruiz 2006
License     :  GPL-style

Maintainer  :  Alberto Ruiz (aruiz at um dot es)
Stability   :  provisional
Portability :  uses ffi

Numerical integration routines.

<http://www.gnu.org/software/gsl/manual/html_node/Numerical-Integration.html#Numerical-Integration>
-}
-----------------------------------------------------------------------------

module Numeric.GSL.Integration (
    integrateQNG,
    integrateQAGS,
    integrateQAGI,
    integrateQAGIU,
    integrateQAGIL
) where

import Foreign.C.Types
import Foreign.Marshal.Alloc(malloc, free)
import Foreign.Ptr(Ptr, FunPtr, freeHaskellFunPtr)
import Foreign.Storable(peek)
import Data.Packed.Internal(check,(//))
import System.IO.Unsafe(unsafePerformIO)

{- | conversion of Haskell functions into function pointers that can be used in the C side
-}
foreign import ccall safe "wrapper" mkfun:: (Double -> Ptr() -> Double) -> IO( FunPtr (Double -> Ptr() -> Double)) 

--------------------------------------------------------------------
{- | Numerical integration using /gsl_integration_qags/ (adaptive integration with singularities). For example:

@\> let quad = integrateQAGS 1E-9 1000 
\> let f a x = x**(-0.5) * log (a*x)
\> quad (f 1) 0 1
(-3.999999999999974,4.871658632055187e-13)@
 
-}

integrateQAGS :: Double               -- ^ precision (e.g. 1E-9)
                 -> Int               -- ^ size of auxiliary workspace (e.g. 1000)
                 -> (Double -> Double) -- ^ function to be integrated on the interval (a,b)
                 -> Double           -- ^ a
                 -> Double           -- ^ b
                 -> (Double, Double) -- ^ result of the integration and error
integrateQAGS prec n f a b = unsafePerformIO $ do
    r <- malloc
    e <- malloc
    fp <- mkfun (\x _ -> f x) 
    c_integrate_qags fp a b prec (fromIntegral n) r e // check "integrate_qags"
    vr <- peek r
    ve <- peek e
    let result = (vr,ve)
    free r
    free e
    freeHaskellFunPtr fp
    return result

foreign import ccall safe "gsl-aux.h integrate_qags" 
 c_integrate_qags :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> Double -> CInt
                     -> Ptr Double -> Ptr Double -> IO CInt

-----------------------------------------------------------------
{- | Numerical integration using /gsl_integration_qng/ (useful for fast integration of smooth functions). For example:

@\> let quad = integrateQNG 1E-6 
\> quad (\\x -> 4\/(1+x*x)) 0 1 
(3.141592653589793,3.487868498008632e-14)@
 
-}

integrateQNG :: Double               -- ^ precision (e.g. 1E-9)
                 -> (Double -> Double) -- ^ function to be integrated on the interval (a,b)
                 -> Double           -- ^ a
                 -> Double           -- ^ b
                 -> (Double, Double) -- ^ result of the integration and error
integrateQNG prec f a b = unsafePerformIO $ do
    r <- malloc
    e <- malloc
    fp <- mkfun (\x _ -> f x) 
    c_integrate_qng fp a b prec r e  // check "integrate_qng"
    vr <- peek r
    ve <- peek e
    let result = (vr,ve)
    free r
    free e
    freeHaskellFunPtr fp
    return result

foreign import ccall safe "gsl-aux.h integrate_qng" 
 c_integrate_qng :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> Double 
                    -> Ptr Double -> Ptr Double -> IO CInt

--------------------------------------------------------------------
{- | Numerical integration using /gsl_integration_qagi/ (integration over the infinite integral -Inf..Inf using QAGS). 
For example:

@\> let quad = integrateQAGI 1E-9 1000 
\> let f a x = exp(-a * x^2)
\> quad (f 0.5) 
(2.5066282746310002,6.229215880648858e-11)@
 
-}

integrateQAGI :: Double               -- ^ precision (e.g. 1E-9)
                 -> Int               -- ^ size of auxiliary workspace (e.g. 1000)
                 -> (Double -> Double) -- ^ function to be integrated on the interval (-Inf,Inf)
                 -> (Double, Double) -- ^ result of the integration and error
integrateQAGI prec n f = unsafePerformIO $ do
    r <- malloc
    e <- malloc
    fp <- mkfun (\x _ -> f x) 
    c_integrate_qagi fp prec (fromIntegral n) r e // check "integrate_qagi"
    vr <- peek r
    ve <- peek e
    let result = (vr,ve)
    free r
    free e
    freeHaskellFunPtr fp
    return result

foreign import ccall safe "gsl-aux.h integrate_qagi" 
 c_integrate_qagi :: FunPtr (Double-> Ptr() -> Double) -> Double -> CInt
                     -> Ptr Double -> Ptr Double -> IO CInt

--------------------------------------------------------------------
{- | Numerical integration using /gsl_integration_qagiu/ (integration over the semi-infinite integral a..Inf). 
For example:

@\> let quad = integrateQAGIU 1E-9 1000 
\> let f a x = exp(-a * x^2)
\> quad (f 0.5) 0
(1.2533141373155001,3.114607940324429e-11)@
 
-}

integrateQAGIU :: Double               -- ^ precision (e.g. 1E-9)
                 -> Int               -- ^ size of auxiliary workspace (e.g. 1000)
                 -> (Double -> Double) -- ^ function to be integrated on the interval (a,Inf)
                 -> Double             -- ^ a
                 -> (Double, Double) -- ^ result of the integration and error
integrateQAGIU prec n f a = unsafePerformIO $ do
    r <- malloc
    e <- malloc
    fp <- mkfun (\x _ -> f x) 
    c_integrate_qagiu fp a prec (fromIntegral n) r e // check "integrate_qagiu"
    vr <- peek r
    ve <- peek e
    let result = (vr,ve)
    free r
    free e
    freeHaskellFunPtr fp
    return result

foreign import ccall safe "gsl-aux.h integrate_qagiu" 
 c_integrate_qagiu :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> CInt
                     -> Ptr Double -> Ptr Double -> IO CInt

--------------------------------------------------------------------
{- | Numerical integration using /gsl_integration_qagil/ (integration over the semi-infinite integral -Inf..b). 
For example:

@\> let quad = integrateQAGIL 1E-9 1000 
\> let f a x = exp(-a * x^2)
\> quad (f 0.5) 0 
(1.2533141373155001,3.114607940324429e-11)@
 
-}

integrateQAGIL :: Double               -- ^ precision (e.g. 1E-9)
                 -> Int               -- ^ size of auxiliary workspace (e.g. 1000)
                 -> (Double -> Double) -- ^ function to be integrated on the interval (a,Inf)
                 -> Double             -- ^ b
                 -> (Double, Double) -- ^ result of the integration and error
integrateQAGIL prec n f b = unsafePerformIO $ do
    r <- malloc
    e <- malloc
    fp <- mkfun (\x _ -> f x) 
    c_integrate_qagil fp b prec (fromIntegral n) r e // check "integrate_qagil"
    vr <- peek r
    ve <- peek e
    let result = (vr,ve)
    free r
    free e
    freeHaskellFunPtr fp
    return result

foreign import ccall safe "gsl-aux.h integrate_qagil" 
 c_integrate_qagil :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> CInt
                     -> Ptr Double -> Ptr Double -> IO CInt