Merge pull request #9 from cdhf/master

Support for integration over infinite intervals

2 files changed, 145 insertions, 1 deletions

diff --git a/lib/Numeric/GSL/Integration.hs b/lib/Numeric/GSL/Integration.hs
index a0e922b..2330fc6 100644
--- a/lib/Numeric/GSL/Integration.hs
+++ b/lib/Numeric/GSL/Integration.hs
@@ -17,7 +17,10 @@ Numerical integration routines.
 
 module Numeric.GSL.Integration (
     integrateQNG,
-    integrateQAGS
+    integrateQAGS,
+    integrateQAGI,
+    integrateQAGIU,
+    integrateQAGIL
 ) where
 
 import Foreign.C.Types
@@ -94,3 +97,102 @@ integrateQNG prec f a b = unsafePerformIO $ do
 foreign import ccall "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 "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 "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 "gsl-aux.h integrate_qagil" 
+ c_integrate_qagil :: FunPtr (Double-> Ptr() -> Double) -> Double -> Double -> CInt
+                     -> Ptr Double -> Ptr Double -> IO CInt
+
diff --git a/lib/Numeric/GSL/gsl-aux.c b/lib/Numeric/GSL/gsl-aux.c
index 33d7dab..5c17836 100644
--- a/lib/Numeric/GSL/gsl-aux.c
+++ b/lib/Numeric/GSL/gsl-aux.c
@@ -761,6 +761,48 @@ int integrate_qags(double f(double,void*), double a, double b, double prec, int
     OK
 }
 
+int integrate_qagi(double f(double,void*), double prec, int w, 
+               double *result, double* error) {
+    DEBUGMSG("integrate_qagi");
+    gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w);
+    gsl_function F;
+    F.function = f;
+    F.params = NULL;
+    int res = gsl_integration_qagi (&F, 0, prec, w,wk, result, error); 
+    CHECK(res,res);
+    gsl_integration_workspace_free (wk); 
+    OK
+}
+
+
+int integrate_qagiu(double f(double,void*), double a, double prec, int w, 
+               double *result, double* error) {
+    DEBUGMSG("integrate_qagiu");
+    gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w);
+    gsl_function F;
+    F.function = f;
+    F.params = NULL;
+    int res = gsl_integration_qagiu (&F, a, 0, prec, w,wk, result, error); 
+    CHECK(res,res);
+    gsl_integration_workspace_free (wk); 
+    OK
+}
+
+
+int integrate_qagil(double f(double,void*), double b, double prec, int w, 
+               double *result, double* error) {
+    DEBUGMSG("integrate_qagil");
+    gsl_integration_workspace * wk = gsl_integration_workspace_alloc (w);
+    gsl_function F;
+    F.function = f;
+    F.params = NULL;
+    int res = gsl_integration_qagil (&F, b, 0, prec, w,wk, result, error); 
+    CHECK(res,res);
+    gsl_integration_workspace_free (wk); 
+    OK
+}
+
+
 int polySolve(KRVEC(a), CVEC(z)) {
     DEBUGMSG("polySolve");
     REQUIRES(an>1,BAD_SIZE);