Added wrapper for QAWC

Cargo.toml

@@ -19,6 +19,7 @@ travis-ci = { repository = "jhod0/integrators", branch = "master" }
 default = ["cuba", "gsl"]
 cuba = []
 gsl = []
+no_gslcblas=[]
 
 [build-dependencies]
 bindgen = "0.43.*"

build.rs

@@ -15,7 +15,9 @@ fn main() {
 
     if env::var("CARGO_FEATURE_GSL").is_ok() {
         println!("cargo:rustc-link-lib=gsl");
-        println!("cargo:rustc-link-lib=gslcblas");
+        if env::var("CARGO_FEATURE_NO_GSLCBLAS").is_err(){
+            println!("cargo:rustc-link-lib=gslcblas");
+        }
         bindings = bindings.header("gsl_wrapper.h");
         any_bindings = true;
     }

src/gsl/mod.rs

@@ -52,6 +52,9 @@ pub use self::qng::QNG;
 mod qag;
 pub use self::qag::{QAG, QAGRule};
 
+mod qawc;
+pub use self::qawc::QAWC;
+
 mod qags;
 pub use self::qags::QAGS;
 

src/gsl/qawc.rs

@@ -0,0 +1,84 @@
+use ::bindings;
+use ::{IntegrationResult, Integrator, Real};
+use ::ffi::LandingPad;
+use ::traits::{IntegrandInput, IntegrandOutput};
+
+use super::{make_gsl_function, GSLIntegrationError, GSLIntegrationWorkspace};
+
+/// Quadrature Adaptive Weighed integration for Cauchy principal values
+/// Calculates the principal value of a function f(x) with a singular weight
+/// 1 / (x - c) centered at a point c
+#[derive(Debug, Clone)]
+pub struct QAWC {
+    range_low: Real,
+    range_high: Real,
+    c: Real,
+    wkspc: GSLIntegrationWorkspace
+}
+
+impl QAWC {
+    /// Creates a new QAGS with enough memory for `nintervals` subintervals.
+    /// This will create a QAG to integrate the range [0, 1].
+    /// To change the integration bounds, see `with_range`.
+    /// The singularity is set to 0.5 by default and can be changed with `with_singularity`
+    pub fn new(nintervals: usize) -> Self{
+        QAWC {
+            range_low: 0.0,
+            range_high: 1.0,
+            c: 0.5,
+            wkspc: GSLIntegrationWorkspace::new(nintervals)
+        }
+    }
+
+    /// Discards the old workspace and allocates a new one with enough memory
+    /// for `nintervals` subintervals.
+    pub fn with_nintervals(self, nintervals: usize) -> Self{
+        QAWC {
+            wkspc: GSLIntegrationWorkspace::new(nintervals),
+            ..self
+        }
+    }
+
+    pub fn with_range(self, range_low: Real, range_high: Real) -> Self {
+        QAWC { range_low, range_high, ..self }
+    }
+
+    pub fn with_singularity(self, c: Real) -> Self{
+        QAWC {c, ..self}
+    }
+}
+
+impl Integrator for QAWC {
+    type Success = IntegrationResult;
+    type Failure = GSLIntegrationError;
+
+    fn integrate<A, B, F: FnMut(A) -> B>(&mut self, fun: F, epsrel: Real, epsabs: Real) -> Result<Self::Success, Self::Failure>
+        where A: IntegrandInput,
+              B: IntegrandOutput
+    {
+        let mut value: Real = 0.0;
+        let mut error: Real = 0.0;
+
+        let mut lp = LandingPad::new(fun);
+        let retcode = unsafe {
+            let mut gslfn = make_gsl_function(&mut lp, self.range_low, self.range_high)?;
+            bindings::gsl_integration_qawc(&mut gslfn.function,
+                                           self.range_low, self.range_high,
+                                           self.c,
+                                           epsabs, epsrel,
+                                           self.wkspc.nintervals,
+                                           self.wkspc.wkspc,
+                                           &mut value,
+                                           &mut error)
+        };
+        lp.maybe_resume_unwind();
+
+        if retcode != bindings::GSL_SUCCESS {
+            Err(GSLIntegrationError::GSLError(retcode.into()))
+        } else {
+            Ok(IntegrationResult {
+                value, error
+            })
+        }
+    }
+}

src/gsl/test.rs

@@ -1,7 +1,7 @@
 //use std::intrinsics::unchecked_div;
 use ::Real;
 use ::Integrator;
-use super::{GSLIntegrationError, QNG, QAG, QAGS, QAGP};
+use super::{GSLIntegrationError, QNG, QAG, QAGS, QAGP, QAWC};
 
 fn nan(_: Real) -> Real {
     ::std::f64::NAN
@@ -124,3 +124,22 @@ fn test_qagp_singularity() {
                    .expect("integration should succeed");
     assert!((res2.value - 2f64).abs() <= res2.error);
 }
+
+#[test]
+fn test_qawc(){
+    let mut qawc = QAWC::new(1000)
+        .with_range(-1., 1.)
+        .with_singularity(0.0);
+    // Trivial example
+    let res1 = qawc.integrate(|x: Real| x, 1e-6, 1e-10)
+        .expect("integration should succeed");
+    assert!((res1.value - 2.0f64).abs() <= res1.error);
+
+    let mut qawc = qawc.with_range(0.0, 2.0)
+        .with_singularity(1.0);
+    // Analytical value: 2 E SinhIntegral[1] ~ 5.74781168
+    let res2 = qawc.integrate(|x: Real| x.exp(), 1e-6, 1e-10)
+        .expect("integration should succeed");
+    println!("Result: {},   err:{}", res2.value, res2.error);
+    assert!((res2.value - 5.747811685312522940687587).abs() <= res2.error);
+}