diff --git a/.astylerc b/.astylerc
new file mode 100644
index 0000000..c2859e0
--- /dev/null
+++ b/.astylerc
@@ -0,0 +1,4 @@
+--style=kr
+--indent=spaces=4
+--convert-tabs
+
diff --git a/README.md b/README.md
index 0a14dd0..120be30 100644
--- a/README.md
+++ b/README.md
@@ -1,35 +1,69 @@
-## spline
-c++ cubic spline library
+## C++ cubic spline interpolation
 
-http://kluge.in-chemnitz.de/opensource/spline/
-
-It generates a piecewise polynomial function of degree 3 and is
-twice continuously differentiable everywhere. Boundary conditions
-default to zero-curvature at the end points. It extrapolates linearly,
-if default boundary conditions are used, or otherwise extrapolation
-is a quadratic function.
+This is a lightweight implementation of **cubic splines**
+with the following features:
 
+* available spline types:
+  * cubic (classical) splines: global, twice continuously differentiable (C<sup>2</sup>)
+  * cubic Hermite splines: local, continuously differentiable (C<sup>1</sup>)
+* boundary conditions: first and second order derivatives can be specified
+* extrapolation
+  * is at most a quadratic function
+  * linear when second order derivatives are set to zero at boundary points or when first derivatives are specified
+* monotonicity can be enforced (when input is monotonic as well)
 
 ### Usage
-The library is a header-only file and can be used like this:
+The library is a header-only file with no external dependencies and can
+be used like this:
 
 ```C++
+#include <vector>
 #include "spline.h"
 ...
-int main(int argc, char** argv)
-{
     std::vector<double> X, Y;
     ...
-    tk::spline s;
-    s.set_points(X,Y);    // X needs to be sorted, strictly increasing
-    double value=s(1.5);  // interpolated value at 1.5
+    // default cubic spline (C^2) with natural boundary conditions (f''=0)
+    tk::spline s(X,Y);			// X needs to be strictly increasing
+    double value=s(1.3);		// interpolated value at 1.3
+    double deriv=s.deriv(1,1.3);	// 1st order derivative at 1.3
     ...
 }
 ```
 
-Optionally, boundary condition can be modified via `set_boundary()`,
-which needs to be called before `set_points()`.
-
+The constructor can take more arguments to define the spline, e.g.:
 ```C++
-s.set_boundary(tk::spline::second_deriv,0.0,tk::spline::first_deriv,-2.0,false);
+    // cubic Hermite splines (C^1) with enforced monotonicity and
+    // left curvature equal to 0.0 and right slope equal 1.0
+    tk::spline s(X,Y,tk::spline::cspline_hermite, true,
+                 tk::spline::second_deriv, 0.0,
+                 tk::spline::first_deriv, 1.0);
+}
 ```
+This is identical to (must be called in that order):
+```C++
+    tk::spline s;
+    s.set_boundary(tk::spline::second_deriv, 0.0,
+                   tk::spline::first_deriv, 1.0);
+    s.set_points(X,Y);
+    s.make_monotonic();
+}
+```
+
+### Spline types
+* `tk::spline::cspline`: classical cubic spline
+  * is uniquely determined by the requirement to be twice continuously differentiable (C<sup>2</sup>)
+  * requires solving a sparse equation system
+  * is a global spline in the sense that changing an input point will impact the spline everywhere
+  * setting first order derivatives at the boundary will break C<sup>2</sup> at the boundary
+* `tk::spline::cspline_hermite`: cubic hermite spline
+  * is continuously differentiable (C<sup>1</sup>) and derivatives are specified on every grid point (equal to the finite differences of the 3 adjacent grid points)
+  * is a local spline in the sense that changing grid points will only affect the spline around that grid point in a few adjacent segments
+
+A function to enforce monotonicity is available as well:
+* `tk::spline::make_monotonic()`: will make the spline monotonic if input grid points are monotonic
+  * this function can only be called after `set_points(...)` has been called
+  * it will break C<sup>2</sup> if the original spline was C<sup>2</sup> and not already monotonic
+  * it will break boundary conditions if it was not monotonic in the first or last segment
+
+### References
+http://kluge.in-chemnitz.de/opensource/spline/
diff --git a/examples/Makefile b/examples/Makefile
index 6af1389..59ed62b 100644
--- a/examples/Makefile
+++ b/examples/Makefile
@@ -6,37 +6,54 @@ IALGLIB=-I$(ALGLIB)/src
 
 # compiler setting
 CC=g++
-CFLAGS = ${IALGLIB} -I../src
-CFLAGS += -O2 -DNDEBUG
-CFLAGS += -Wall
+CFLAGS = ${IALGLIB} -I../src -std=c++11 -Wall -Wextra
+CFLAGS += -O2
+#CFLAGS += -DNDEBUG			# disables all assert() checks
+#CFLAGS +=-D_GLIBCXX_USE_CXX11_ABI=0	# for new gcc & boost compiled with old
 LFLAGS=
 
 
-
-PROGS=demo interpol check_alglib interpol_alglib bench
+PROGS=simple_demo plot plot_avg_preserv plot_parametric
+PROGS_ALGLIB=plot_alglib bench
+PROGS_TESTS=tests/unit_tests tests/compile_units
 
 all:		${PROGS}
+more:		${PROGS_ALGLIB}
+tests:		${PROGS_TESTS}
+		./tests/unit_tests
 
-demo:		demo.o
-		${CC} $< -o $@
 
-interpol:	interpol.o
+simple_demo:	simple_demo.o
+		${CC} $< -o $@
+plot:		plot.o
+		${CC} $< -o $@
+plot_avg_preserv:	plot_avg_preserv.o
 		${CC} $< -o $@
+plot_parametric:	plot_parametric.o
+		${CC} $< -o $@
+
+tests/compile_units:	tests/compile_units.o tests/myfunc1.o tests/myfunc2.o
+		${CC} $? -o $@ 
+tests/unit_tests:	tests/unit_tests.o
+		${CC} $< -o $@ -lboost_unit_test_framework
+
+%.o: 		%.cpp ../src/spline.h
+		${CC} ${CFLAGS} -c $< -o $@
+
+tests/%.o:	tests/%.cpp ../src/spline.h
+		${CC} ${CFLAGS} -Werror -c $< -o $@
+
 
 check_alglib:
 ifeq ($(wildcard $(ALGLIB)),)
-  $(error Please install ALGLIB and adjust the corresponding variables within this Makefile)
+	$(error Please install ALGLIB and adjust the corresponding variables within this Makefile)
 endif
 
-interpol_alglib:	interpol_alglib.o
-		${CC} $< -o $@ ${LALGLIB}
+plot_alglib:	check_alglib plot_alglib.o
+		${CC} plot_alglib.o -o $@ ${LALGLIB}
+bench:		check_alglib bench.o
+		${CC} bench.o -o $@ ${LALGLIB}
 
-bench:		bench.o
-		${CC} $< -o $@ ${LALGLIB}
-
-
-%.o: %.cpp
-		${CC} ${CFLAGS} -c $< -o $@
 
 clean:
-		rm -f *.o ${PROGS} *.csv
+		rm -f *.o tests/*.o ${PROGS} ${PROGS_ALGLIB} ${PROGS_TESTS}  *.csv
diff --git a/examples/README.md b/examples/README.md
new file mode 100644
index 0000000..91043fe
--- /dev/null
+++ b/examples/README.md
@@ -0,0 +1,47 @@
+## Spline example use cases
+
+### Prerequisites
+* C++11 compatible compiler
+* For a few examples [ALGLIB](https://www.alglib.net/) is required
+
+### Building
+Adjust the variables `CC`, `CFLAGS` and `ALGLIB` in the `Makefile` and
+execute as follows:
+```
+$ make 		# builds all self-contained examples
+$ make more	# builds all examples which require ALGLIB
+$ make tests	# builds and runs several tests
+```
+
+### Example files (self contained)
+
+* `simple_demo.cpp`: sets up a spline and evaluates it at a single point
+* `plot*.cpp`: there are several examples to plot splines in different situations.
+	They all produce numerical output which can be used with [Gnuplot](http://www.gnuplot.info/) for visualisation. Shell scripts are provided to do this in one step: `./plot*.sh` (or manually: `./plot > plot.csv`, `gnuplot plot.gp`).
+  * `plot.cpp`: plot splines, e.g.
+	```
+	$ ./plot.sh -t hermite -x 0.0 0.5 1.5 2.0 -y 1.1 0.9 0.5 0.7
+	```
+  * `plot_avg_preserv.cpp`: plots **average preserving** (area preserving/mass preserving) splines. These are splines which do not interpolate points but which have an average value in the interval [x<sub>i</sub>, x<sub>i+1</sub>] equal to y<sub>i</sub>.
+	```
+	$ ./plot_avg_preserv.sh -t cspline 
+	```
+  * `plot_parametric.cpp`: plots a **parametric spline**, i.e. connects points in a 2D-plane, where the x-coordinates do no longer need to be strictly increasing.
+	```
+	$ ./plot_parametric.sh -t hermite
+	```
+
+### Example files (which require ALGLIB)
+* `bench.cpp`: performs a benchmark of creating a spline and interpolating values. E.g. on a 2.4GHz CPU, spline size 50 and 1 million spline evaluations call:
+	```
+	$ ./bench 2400 50 1000
+	                                tk                      alglib
+	random access:   loops=1e+06,   0.064s ( 153 cycl)      0.078s ( 188 cycl)
+	spline creation: loops=2e+04,   0.200s (2.4e+04 cycl)   0.202s (2.4e+04 cycl)
+	grid transform:  loops=2e+04,   0.250s (3.0e+04 cycl)   0.193s (2.3e+04 cycl)
+	```
+* `plot_alglib.cpp`: plots ALGLIB spline against this spline implementation (does not have any command line options)
+	```
+	$ ./plot_alglib.sh
+	```
+
diff --git a/examples/demo.cpp b/examples/demo.cpp
deleted file mode 100644
index 59eb293..0000000
--- a/examples/demo.cpp
+++ /dev/null
@@ -1,29 +0,0 @@
-#include <cstdio>
-#include <cstdlib>
-#include <vector>
-#include "spline.h"
-
-int main(int argc, char** argv)
-{
-
-    std::vector<double> X(5), Y(5);
-    X[0]=0.1;
-    X[1]=0.4;
-    X[2]=1.2;
-    X[3]=1.8;
-    X[4]=2.0;
-    Y[0]=0.1;
-    Y[1]=0.7;
-    Y[2]=0.6;
-    Y[3]=1.1;
-    Y[4]=0.9;
-
-    tk::spline s;
-    s.set_points(X,Y);    // currently it is required that X is already sorted
-
-    double x=1.5;
-
-    printf("spline at %f is %f\n", x, s(x));
-
-    return EXIT_SUCCESS;
-}
diff --git a/examples/interpol.cpp b/examples/interpol.cpp
deleted file mode 100644
index ea531ad..0000000
--- a/examples/interpol.cpp
+++ /dev/null
@@ -1,73 +0,0 @@
-#include <cstdio>
-#include <cstdlib>
-#include <cmath>
-#include <vector>
-#include "spline.h"
-
-template<class Function>
-double deriv1(const Function& f, double x)
-{
-    double dx=1e-8*(1.0+fabs(x));
-    double x0=x-dx;
-    double x2=x+dx;
-    return (f(x2)-f(x0)) / (2.0*dx);
-}
-
-template<class Function>
-double deriv2(const Function& f, double x)
-{
-    double dx=1e-6*(1.0+fabs(x));
-    double x0=x-dx;
-    double x2=x+dx;
-    return (f(x0)-2.0*f(x)+f(x2)) / (dx*dx);
-}
-
-
-
-int main(int argc, char** argv)
-{
-    if(argc<1) {
-        printf("usage: %s <>\n", argv[0]);
-        exit(EXIT_FAILURE);
-    }
-    std::vector<double> X(5), Y(5);
-    X[0]=0.1;
-    X[1]=0.4;
-    X[2]=1.2;
-    X[3]=1.8;
-    X[4]=2.0;
-    Y[0]=0.1;
-    Y[1]=0.7;
-    Y[2]=0.6;
-    Y[3]=1.1;
-    Y[4]=0.9;
-
-    tk::spline s;
-    // set_boundary() is optional and if omitted, natural boundary condition,
-    // f''(a)=f''(b)=0, will be used
-    // note, if natural boundary conditions are not used then extrapolation
-    // will be a quadratic function, unless the last argument is set to true,
-    // which forces linear extrapolation, but this will violate second order
-    // differentiability at the endpoint
-    s.set_boundary(tk::spline::second_deriv, 0.0,
-                   tk::spline::first_deriv, -2.0, false);
-    s.set_points(X,Y);
-
-    for(size_t i=0; i<X.size(); i++) {
-        printf("%f %f\n", X[i], Y[i]);
-    }
-    printf("\n");
-    for(int i=-50; i<250; i++) {
-        double x=0.01*i;
-        printf("%f %f %f %f %f\n", x, s(x),
-               s.deriv(1,x), s.deriv(2,x), s.deriv(3,x));
-        // checking analytic derivatives and finite differences are close
-        assert(fabs(s.deriv(1,x)-deriv1(s,x)) < 1e-8);
-        assert(fabs(s.deriv(2,x)-deriv2(s,x)) < 1e-8);
-    }
-    printf("\n");
-
-    return EXIT_SUCCESS;
-}
-
-
diff --git a/examples/interpol_alglib.cpp b/examples/interpol_alglib.cpp
deleted file mode 100644
index 402c1ac..0000000
--- a/examples/interpol_alglib.cpp
+++ /dev/null
@@ -1,47 +0,0 @@
-#include <cstdio>
-#include <cstdlib>
-#include <vector>
-#include "interpolation.h" // alglib
-
-int main(int argc, char** argv)
-{
-    if(argc<1) {
-        printf("usage: %s <>\n", argv[0]);
-        exit(EXIT_FAILURE);
-    }
-    std::vector<double> X(5), Y(5);
-    X[0]=0.1;
-    X[1]=0.4;
-    X[2]=1.2;
-    X[3]=1.8;
-    X[4]=2.0;
-    Y[0]=0.1;
-    Y[1]=0.7;
-    Y[2]=0.6;
-    Y[3]=1.1;
-    Y[4]=0.9;
-
-    alglib::real_1d_array AX, AY;
-    AX.setcontent(X.size(), &(X[0]));
-    AY.setcontent(Y.size(), &(Y[0]));
-
-    alglib::spline1dinterpolant spline;
-    alglib::spline1dbuildcubic(AX, AY, X.size(), 2,0.0,2,0.0, spline);
-    //alglib::spline1dbuildcubic(AX, AY, spline);
-
-
-    for(size_t i=0; i<X.size(); i++) {
-        printf("%f %f\n", X[i], Y[i]);
-    }
-
-    printf("\n");
-    for(int i=-50; i<250; i++) {
-        double x=0.01*i;
-        printf("%f %f\n", x, alglib::spline1dcalc(spline,x));
-    }
-    printf("\n");
-
-    return EXIT_SUCCESS;
-}
-
-
diff --git a/examples/plot.cpp b/examples/plot.cpp
new file mode 100644
index 0000000..8303a1d
--- /dev/null
+++ b/examples/plot.cpp
@@ -0,0 +1,171 @@
+// plot splines: generates output which can be plotted using gnuplot
+
+#include <cstdio>
+#include <cstdlib>
+#include <cmath>
+#include <vector>
+#include <getopt.h>
+#include "spline.h"
+
+void print_usage(const char* progname)
+{
+    printf( "usage: %s: \n"
+            "\t[-h] help\n"
+            "\t[-t <spline type: cspline, hermite, linear>]\n"
+            "\t[-m] force piece-wise monotonicity \n"
+            "\t[-x <x1> <x2> <x3> [x4] ... [xn]]\n"
+            "\t[-y <y1> <y2> <y3> [y4] ... [yn]] \n", progname);
+}
+
+void parse_args(int argc, char** argv,
+                std::vector<double>& X, std::vector<double>& Y,
+                tk::spline::spline_type& type, bool& make_monotonic)
+{
+    int opt;
+    bool x_provided=false;
+    bool y_provided=false;
+    while( (opt = getopt(argc, argv, "ht:mx:y:")) != -1) {
+        switch( opt ) {
+        case 'h':
+            print_usage(argv[0]);
+            exit(EXIT_SUCCESS);
+            break;
+        case 't':
+            if(std::string(optarg)=="cspline") {
+                type=tk::spline::cspline;
+            } else if (std::string(optarg)=="hermite") {
+                type=tk::spline::cspline_hermite;
+            } else if (std::string(optarg)=="linear") {
+                type=tk::spline::linear;
+            } else {
+                print_usage(argv[0]);
+                exit(EXIT_FAILURE);
+            }
+            break;
+        case 'm':
+            make_monotonic=true;
+            break;
+        case 'x':
+            x_provided=true;
+            X.resize(0);
+            optind--;
+            for( ; optind < argc && *argv[optind] != '-'; optind++) {
+                X.push_back(atof(argv[optind]));
+            }
+            break;
+        case 'y':
+            y_provided=true;
+            Y.resize(0);
+            optind--;
+            for( ; optind < argc && *argv[optind] != '-'; optind++) {
+                Y.push_back(atof(argv[optind]));
+            }
+            break;
+        default:
+            // we won't actually get here
+            print_usage(argv[0]);
+            exit(EXIT_FAILURE);
+            break;
+        }
+    }
+    // create uniform x-grid if only y-grid has been provided
+    if(x_provided==false && y_provided==true) {
+        X.resize(Y.size());
+        for(size_t i=0; i<X.size(); i++) {
+            X[i] = (double)i;
+        }
+    }
+    if(x_provided==true && y_provided==false) {
+        fprintf(stderr, "X vector provided but not Y\n");
+        exit(EXIT_FAILURE);
+    }
+    // check that X and Y have the correct size
+    if(X.size()!=Y.size() || X.size()<3) {
+        fprintf(stderr, "X and Y need to have the same number of elements and at least 3 points\n");
+        fprintf(stderr, "provided inputs:\n");
+        fprintf(stderr, "X = {");
+        for(size_t i=0; i<X.size(); i++) {
+            fprintf(stderr, "%.3f, ", X[i]);
+        }
+        fprintf(stderr, "}\n");
+        fprintf(stderr, "Y = {");
+        for(size_t i=0; i<Y.size(); i++) {
+            fprintf(stderr, "%.3f, ", Y[i]);
+        }
+        fprintf(stderr, "}\n");
+        exit(EXIT_FAILURE);
+    }
+    // check that X is monotonic
+    for(int i=1; i<(int)X.size(); i++) {
+        if(X[i-1]>=X[i]) {
+            fprintf(stderr, "input X needs to be strictly increasing, but");
+            fprintf(stderr, " X[%i]=%.3f, X[%i]=%.3f\n", i-1, X[i-1],i,X[i]);
+            exit(EXIT_FAILURE);
+        }
+    }
+}
+
+// finite differences for 1st and 2nd order derivatives
+template<class Function>
+double fderiv(const Function& f, int order, double x)
+{
+    if(order==1) {
+        double dx=2e-8*(1.0+fabs(x));
+        return (f(x+dx)-f(x-dx)) / (2.0*dx);
+    } else if(order==2) {
+        double dx=3e-6*(1.0+fabs(x));
+        return (f(x-dx)-2.0*f(x)+f(x+dx)) / (dx*dx);
+    } else {
+        assert(false);
+        return -1.0;
+    }
+}
+
+
+
+int main(int argc, char** argv)
+{
+    // default parameters
+    std::vector<double> X = {0.1, 0.4, 1.2, 1.8, 2.0};
+    std::vector<double> Y = {0.1, 0.7, 0.75, 1.1, 0.9};
+    tk::spline::spline_type type = tk::spline::cspline;
+    bool make_monotonic = false;
+
+    // override defaults with command line arguments if supplied
+    parse_args(argc, argv, X, Y, type, make_monotonic);
+
+    // setup spline
+    tk::spline s;
+    // boundary conditions: natural condition is 2nd deriv = 0
+    // extrapolation:
+    //  - at most quadratic and linear if 2nd deriv = 0
+    //  - linear if 1st derivative at boundary is specified (this breaks C^2)
+    s.set_boundary(tk::spline::second_deriv, 0.0,
+                   tk::spline::second_deriv, 0.0);
+    s.set_points(X,Y,type);     // this calculates all spline coefficients
+    if(make_monotonic) {
+        s.make_monotonic();     // adjusts spline coeffs to be monotonic
+    }
+
+    // evaluates spline and outputs data to be used with gnuplot
+    printf("# input x, input y\n");
+    for(size_t i=0; i<X.size(); i++) {
+        printf("%f %f\n", X[i], Y[i]);  // input grid points
+    }
+    printf("\n");
+    int n = 1000;    // number of grid points to plot the spline
+    double xmin = X[0] - 0.5;
+    double xmax = X.back() + 0.5;
+    printf("# x, s(x), s'(x), s''(x), s'''(x), fdiff(s,1,x), fdiff(s,2,x)\n");
+    for(int i=0; i<n; i++) {
+        double x = xmin + (double)i*(xmax-xmin)/(n-1);
+        printf("%f %f %f %f %f %f %f\n", x, s(x),
+               s.deriv(1,x), s.deriv(2,x), s.deriv(3,x),
+               fderiv(s,1,x), fderiv(s,2,x));
+    }
+    printf("\n");
+
+    return EXIT_SUCCESS;
+}
+
+
diff --git a/examples/plot.gp b/examples/plot.gp
new file mode 100644
index 0000000..e0a651a
--- /dev/null
+++ b/examples/plot.gp
@@ -0,0 +1,31 @@
+# gnuplot script file
+
+#set terminal x11 noreplotonresize
+#set terminal qt size 800,800
+
+set multiplot 
+set xzeroaxis
+set yzeroaxis
+set key top right
+set lmargin at screen 0.05	# align all plots to the left
+
+
+set size 1.0,0.4
+set origin 0.0,0.6
+plot	"plot.csv" every :::0::0 notitle with p ps 2, \
+	"plot.csv" every :::1::1 title "spline" with l lt 1
+
+set size 1.0,0.2
+set origin 0.0,0.4
+plot	"plot.csv" every :::1::1 using 1:3 title "1st derivative" with l lt 1,\
+	"plot.csv" every :::1::1 using 1:6 notitle with l lt 1 # should be same
+
+set origin 0.0,0.2
+plot	"plot.csv" every :::1::1 using 1:4 title "2nd derivative" with l lt 1,\
+	"plot.csv" every :::1::1 using 1:7 notitle with l lt 1 # should be same
+
+set origin 0.0,0.0
+plot	"plot.csv" every :::1::1 using 1:5 title "3rd derivative" with l lt 1
+
+unset multiplot
+pause -1
diff --git a/examples/plot.sh b/examples/plot.sh
index c668786..d962633 100755
--- a/examples/plot.sh
+++ b/examples/plot.sh
@@ -1,5 +1,6 @@
 #!/bin/bash
-
-./interpol > plot.csv
-./interpol_alglib >> plot.csv
-gnuplot plot.txt
+./plot $@ > plot.csv
+return_value=$?
+if [ $return_value -eq 0 ] ; then
+	gnuplot plot.gp
+fi
diff --git a/examples/plot.txt b/examples/plot.txt
deleted file mode 100644
index 7b01433..0000000
--- a/examples/plot.txt
+++ /dev/null
@@ -1,14 +0,0 @@
-set multiplot layout 3, 1
-
-set key top left
-plot	"plot.csv" every :::0::0 notitle with p ps 2, \
-	"plot.csv" every :::1::1 title "tk" with l lt 1, \
-	"plot.csv" every :::3::3 title "alglib" with l lt 2
-
-
-set key top right
-plot	"plot.csv" every :::1::1 using 1:3 title "tk first deriv" with l lt 1
-
-plot	"plot.csv" every :::1::1 using 1:4 title "tk second deriv" with l lt 1
-
-pause -1
diff --git a/examples/plot_alglib.cpp b/examples/plot_alglib.cpp
new file mode 100644
index 0000000..9449500
--- /dev/null
+++ b/examples/plot_alglib.cpp
@@ -0,0 +1,72 @@
+#include <cstdio>
+#include <cstdlib>
+#include <vector>
+#include "interpolation.h"  // alglib
+#include "spline.h"         // tk
+
+// wrap alglib spline into a class
+class alglib_spline
+{
+private:
+    alglib::spline1dinterpolant m_spline;
+public:
+    alglib_spline(const std::vector<double>& X, const std::vector<double>& Y)
+    {
+        alglib::real_1d_array AX, AY;
+        AX.setcontent(X.size(), &(X[0]));
+        AY.setcontent(Y.size(), &(Y[0]));
+        //alglib::spline1dbuildcubic(AX, AY, m_spline);
+        alglib::spline1dbuildcubic(AX, AY, X.size(), 2,0.0,2,0.0, m_spline);
+    }
+    double operator()(double x) const
+    {
+        return alglib::spline1dcalc(m_spline,x);
+    }
+};
+
+template<class Function>
+double fderiv(const Function& f, int order, double x)
+{
+    if(order==1) {
+        double dx=2e-8*(1.0+fabs(x));
+        return (f(x+dx)-f(x-dx)) / (2.0*dx);
+    } else if(order==2) {
+        double dx=3e-6*(1.0+fabs(x));
+        return (f(x-dx)-2.0*f(x)+f(x+dx)) / (dx*dx);
+    } else {
+        assert(false);
+        return -1.0;
+    }
+}
+
+
+int main(int argc, char** argv)
+{
+    if(argc<1) {
+        printf("usage: %s <>\n", argv[0]);
+        exit(EXIT_FAILURE);
+    }
+    std::vector<double> X = {0.1, 0.4, 1.2, 1.8, 2.0};
+    std::vector<double> Y = {0.1, 0.7, 0.6, 1.1, 0.9};
+
+    alglib_spline   aspline(X,Y);
+    tk::spline      tspline(X,Y);
+
+    printf("# input x, input y\n");
+    for(size_t i=0; i<X.size(); i++) {
+        printf("%f %f\n", X[i], Y[i]);
+    }
+    printf("\n");
+    printf("# x, alglib spline(x), tk spline(x), as'(x), ts'(x), as''(x), ts''(x)\n");
+    for(int i=-50; i<250; i++) {
+        double x=0.01*i;
+        printf("%f %f %f %f %f %f %f\n", x, aspline(x), tspline(x),
+               fderiv(aspline,1,x), fderiv(tspline,1,x),
+               fderiv(aspline,2,x), fderiv(tspline,2,x));
+    }
+    printf("\n");
+
+    return EXIT_SUCCESS;
+}
+
+
diff --git a/examples/plot_alglib.gp b/examples/plot_alglib.gp
new file mode 100644
index 0000000..1a4aa3d
--- /dev/null
+++ b/examples/plot_alglib.gp
@@ -0,0 +1,30 @@
+# gnuplot script file
+
+#set terminal x11 noreplotonresize
+#set terminal qt size 800,800
+
+set multiplot 
+set xzeroaxis
+set yzeroaxis
+set key top right
+set lmargin at screen 0.05	# align all plots to the left
+
+
+set size 1.0,0.5
+set origin 0.0,0.5
+plot	"plot.csv" every :::0::0 notitle with p ps 2, \
+        "plot.csv" every :::1::1 using 1:2 title "alglib spline" with l lt 1,\
+        "plot.csv" every :::1::1 using 1:3 title "tk spline" with l lt 2
+
+
+set size 1.0,0.25
+set origin 0.0,0.25
+plot	"plot.csv" every :::1::1 using 1:4 title "1st deriv alglib" with l lt 1,\
+	"plot.csv" every :::1::1 using 1:5 title "1st deriv tk"  with l lt 2
+
+set origin 0.0,0.0
+plot	"plot.csv" every :::1::1 using 1:6 title "2nd deriv alglib" with l lt 1,\
+	"plot.csv" every :::1::1 using 1:7 title "2nd deriv tk" with l lt 2
+
+unset multiplot
+pause -1
diff --git a/examples/plot_alglib.sh b/examples/plot_alglib.sh
new file mode 100755
index 0000000..9dd807a
--- /dev/null
+++ b/examples/plot_alglib.sh
@@ -0,0 +1,6 @@
+#!/bin/bash
+./plot_alglib $@ > plot.csv
+return_value=$?
+if [ $return_value -eq 0 ] ; then
+	gnuplot plot_alglib.gp
+fi
diff --git a/examples/plot_avg_preserv.cpp b/examples/plot_avg_preserv.cpp
new file mode 100644
index 0000000..e1fb214
--- /dev/null
+++ b/examples/plot_avg_preserv.cpp
@@ -0,0 +1,205 @@
+// plot average preserving splines:
+//      this only generates output which can then be plotted using gnuplot
+
+#include <cstdio>
+#include <cstdlib>
+#include <cmath>
+#include <vector>
+#include <getopt.h>
+#include "spline.h"
+
+void print_usage(const char* progname)
+{
+    printf( "usage: %s: \n"
+            "\t[-h] help\n"
+            "\t[-t <spline type: cspline, hermite, linear>]\n"
+            "\t[-m] force piece-wise positivity/negativity\n"
+            "\t[-x <x1> <x2> <x3> [x4] ... [xn]]\n"
+            "\t[-y <y1> <y2> [y3] ... [yn-1]] \n", progname);
+}
+
+void parse_args(int argc, char** argv,
+                std::vector<double>& X, std::vector<double>& Y,
+                tk::spline::spline_type& type, bool& make_monotonic)
+{
+    int opt;
+    bool x_provided=false;
+    bool y_provided=false;
+    while( (opt = getopt(argc, argv, "ht:mx:y:")) != -1) {
+        switch( opt ) {
+        case 'h':
+            print_usage(argv[0]);
+            exit(EXIT_SUCCESS);
+            break;
+        case 't':
+            if(std::string(optarg)=="cspline") {
+                type=tk::spline::cspline;
+            } else if (std::string(optarg)=="hermite") {
+                type=tk::spline::cspline_hermite;
+            } else if (std::string(optarg)=="linear") {
+                type=tk::spline::linear;
+            } else {
+                print_usage(argv[0]);
+                exit(EXIT_FAILURE);
+            }
+            break;
+        case 'm':
+            make_monotonic=true;
+            break;
+        case 'x':
+            x_provided=true;
+            X.resize(0);
+            optind--;
+            for( ; optind < argc && *argv[optind] != '-'; optind++) {
+                X.push_back(atof(argv[optind]));
+            }
+            break;
+        case 'y':
+            y_provided=true;
+            Y.resize(0);
+            optind--;
+            for( ; optind < argc && *argv[optind] != '-'; optind++) {
+                Y.push_back(atof(argv[optind]));
+            }
+            break;
+        default:
+            // we won't actually get here
+            print_usage(argv[0]);
+            exit(EXIT_FAILURE);
+            break;
+        }
+    }
+    // create uniform x-grid if only y-grid has been provided
+    if(x_provided==false && y_provided==true) {
+        X.resize(Y.size()+1);
+        for(size_t i=0; i<X.size(); i++) {
+            X[i] = (double)i;
+        }
+    }
+    if(x_provided==true && y_provided==false) {
+        fprintf(stderr, "X vector provided but not Y\n");
+        exit(EXIT_FAILURE);
+    }
+    // check that X and Y have the correct size
+    if(X.size()!=Y.size()+1 || X.size()<3) {
+        fprintf(stderr, "X needs to have one element more than Y and Y needs to have at least 2 elements\n");
+        fprintf(stderr, "provided inputs:\n");
+        fprintf(stderr, "X = {");
+        for(size_t i=0; i<X.size(); i++) {
+            fprintf(stderr, "%.3f, ", X[i]);
+        }
+        fprintf(stderr, "}\n");
+        fprintf(stderr, "Y = {");
+        for(size_t i=0; i<Y.size(); i++) {
+            fprintf(stderr, "%.3f, ", Y[i]);
+        }
+        fprintf(stderr, "}\n");
+        exit(EXIT_FAILURE);
+    }
+    // check that X is monotonic
+    for(int i=1; i<(int)X.size(); i++) {
+        if(X[i-1]>=X[i]) {
+            fprintf(stderr, "input X needs to be strictly increasing, but");
+            fprintf(stderr, " X[%i]=%.3f, X[%i]=%.3f\n", i-1, X[i-1],i,X[i]);
+            exit(EXIT_FAILURE);
+        }
+    }
+}
+
+// to generate some sample data points for a histogram
+double rand_unif()
+{
+    return (double)std::rand()/RAND_MAX;
+}
+double my_rand()
+{
+    double sum=0.0;
+    for(int i=0; i<10; i++) {
+        sum+=rand_unif()-0.5;
+    }
+    return sum;
+}
+int find_closest(const std::vector<double> X, double x)
+{
+    std::vector<double>::const_iterator it;
+    it=std::upper_bound(X.begin(),X.end(),x);       // *it > x
+    int idx = int(it-X.begin())-1;                    // m_x[idx] <= x
+    return idx;
+}
+void generate_histogram(std::vector<double>& X, std::vector<double>& avg)
+{
+    const int num_runs = 1000;
+    const int num_buckets = 20;
+    const double x_min = -3.0;
+    const double x_max = 3.0;
+    X.resize(num_buckets+1);
+    avg.resize(num_buckets);
+    for(int i=0; i<=num_buckets; i++) {
+        X[i] = x_min + (double)i * (x_max-x_min)/num_buckets;
+    }
+    for(int i=0; i<num_buckets; i++) {
+        avg[i] = 0.0;
+    }
+
+    for(int i=0; i<num_runs; i++) {
+        double x = my_rand();
+        int idx=find_closest(X, x);
+        if(0<=idx && idx<num_buckets) {
+            avg[idx]++;
+        }
+    }
+    for(int i=0; i<num_buckets; i++) {
+        avg[i]/=num_buckets;
+    }
+}
+
+
+int main(int argc, char** argv)
+{
+    // default parameters
+    std::vector<double> X, avg;
+    generate_histogram(X,avg);
+    tk::spline::spline_type type = tk::spline::cspline;
+    bool make_monotonic = false;
+
+    // override defaults with command line arguments if supplied
+    parse_args(argc, argv, X, avg, type, make_monotonic);
+
+    // build Y-vector as the integral over the averages avg
+    std::vector<double> Y(avg.size()+1);
+    Y[0]=0.0;
+    for(size_t i=1; i<Y.size(); i++) {
+        Y[i] = Y[i-1] + avg[i-1]*(X[i]-X[i-1]);
+    }
+
+    // apply spline to X-Y vectors
+    tk::spline s;
+    s.set_boundary(tk::spline::second_deriv, 0.0,
+                   tk::spline::second_deriv, 0.0);
+    s.set_points(X,Y,type);     // this calculates all spline coefficients
+    if(make_monotonic) {
+        s.make_monotonic();     // adjusts spline coeffs to be monotonic
+    }
+
+    // the derivative of our cubic spline s then preserves the averages avg
+    // this outputs the data to be used with gnuplot
+    printf("# input x, input avg\n");
+    for(size_t i=0; i<avg.size(); i++) {
+        printf("%f %f\n", X[i], avg[i]);  // input grid points
+    }
+    printf("%f %f\n", X.back(), avg.back());
+    printf("\n");
+    int n = 1000;    // number of grid points to plot the spline
+    double xmin = X[0] - 0.5;
+    double xmax = X.back() + 0.5;
+    printf("# x, s'(x), s''(x), s'''(x)\n");
+    for(int i=0; i<n; i++) {
+        double x = xmin + (double)i*(xmax-xmin)/(n-1);
+        printf("%f %f %f %f\n", x, s.deriv(1,x), s.deriv(2,x), s.deriv(3,x));
+    }
+    printf("\n");
+
+    return EXIT_SUCCESS;
+}
+
+
diff --git a/examples/plot_avg_preserv.gp b/examples/plot_avg_preserv.gp
new file mode 100644
index 0000000..1d26178
--- /dev/null
+++ b/examples/plot_avg_preserv.gp
@@ -0,0 +1,33 @@
+# gnuplot script file
+
+#set term x11 noreplotonresize
+#set terminal qt size 800,800
+
+
+set lmargin at screen 0.05      # align all plots to the left
+
+set multiplot
+
+set xzeroaxis
+set yzeroaxis
+set grid y
+set size 1.0,0.5
+set origin 0.0,0.5
+set key top right
+set title "average preserving spline interpolation"
+plot    "plot.csv" every :::0::0 using 1:2 notitle with steps lt 2, \
+        "plot.csv" every :::1::1 using 1:2 title "spline" with l lt 1
+
+
+set size 1.0,0.25
+set origin 0.0,0.25
+set notitle
+plot    "plot.csv" every :::1::1 using 1:3 title "1st derivative" with l lt 1
+        
+set size 1.0,0.25
+set origin 0.0,0.0
+plot    "plot.csv" every :::1::1 using 1:4 title "2nd derivative" with l lt 1
+
+unset multiplot
+
+pause -1
diff --git a/examples/plot_avg_preserv.sh b/examples/plot_avg_preserv.sh
new file mode 100755
index 0000000..c47f170
--- /dev/null
+++ b/examples/plot_avg_preserv.sh
@@ -0,0 +1,6 @@
+#!/bin/bash
+./plot_avg_preserv $@ > plot.csv
+return_value=$?
+if [ $return_value -eq 0 ] ; then
+	gnuplot plot_avg_preserv.gp
+fi
diff --git a/examples/plot_parametric.cpp b/examples/plot_parametric.cpp
new file mode 100644
index 0000000..0a4b3f5
--- /dev/null
+++ b/examples/plot_parametric.cpp
@@ -0,0 +1,156 @@
+// plot parametric splines
+//  this generates output which can be plotted using gnuplot
+
+#include <cstdio>
+#include <cstdlib>
+#include <cmath>
+#include <vector>
+#include <getopt.h>
+#include "spline.h"
+
+void print_usage(const char* progname)
+{
+    printf( "usage: %s: \n"
+            "\t[-h] help\n"
+            "\t[-t <spline type: cspline, hermite, linear>]\n"
+            "\t[-m] force piece-wise monotonicity \n"
+            "\t[-x <x1> <x2> <x3> [x4] ... [xn]]\n"
+            "\t[-y <y1> <y2> <y3> [y4] ... [yn]] \n", progname);
+}
+
+void parse_args(int argc, char** argv,
+                std::vector<double>& X, std::vector<double>& Y,
+                tk::spline::spline_type& type, bool& make_monotonic)
+{
+    int opt;
+    bool x_provided=false;
+    bool y_provided=false;
+    while( (opt = getopt(argc, argv, "ht:mx:y:")) != -1) {
+        switch( opt ) {
+        case 'h':
+            print_usage(argv[0]);
+            exit(EXIT_SUCCESS);
+            break;
+        case 't':
+            if(std::string(optarg)=="cspline") {
+                type=tk::spline::cspline;
+            } else if (std::string(optarg)=="hermite") {
+                type=tk::spline::cspline_hermite;
+            } else if (std::string(optarg)=="linear") {
+                type=tk::spline::linear;
+            } else {
+                print_usage(argv[0]);
+                exit(EXIT_FAILURE);
+            }
+            break;
+        case 'm':
+            make_monotonic=true;
+            break;
+        case 'x':
+            x_provided=true;
+            X.resize(0);
+            optind--;
+            for( ; optind < argc && *argv[optind] != '-'; optind++) {
+                X.push_back(atof(argv[optind]));
+            }
+            break;
+        case 'y':
+            y_provided=true;
+            Y.resize(0);
+            optind--;
+            for( ; optind < argc && *argv[optind] != '-'; optind++) {
+                Y.push_back(atof(argv[optind]));
+            }
+            break;
+        default:
+            // we won't actually get here
+            print_usage(argv[0]);
+            exit(EXIT_FAILURE);
+            break;
+        }
+    }
+    if(x_provided!=y_provided) {
+        fprintf(stderr, "need to provide both X and Y vectors or none\n");
+        exit(EXIT_FAILURE);
+    }
+
+    // check that X and Y have the correct size
+    if(X.size()!=Y.size() || X.size()<3) {
+        fprintf(stderr, "X and Y need to have the same number of elements and at least 3 points\n");
+        fprintf(stderr, "provided inputs:\n");
+        fprintf(stderr, "X = {");
+        for(size_t i=0; i<X.size(); i++) {
+            fprintf(stderr, "%.3f, ", X[i]);
+        }
+        fprintf(stderr, "}\n");
+        fprintf(stderr, "Y = {");
+        for(size_t i=0; i<Y.size(); i++) {
+            fprintf(stderr, "%.3f, ", Y[i]);
+        }
+        fprintf(stderr, "}\n");
+        exit(EXIT_FAILURE);
+    }
+    // check that two succesive points are not identical
+    for(int i=1; i<(int)X.size(); i++) {
+        if( X[i-1]==X[i] && Y[i-1]==Y[i]) {
+            fprintf(stderr, "two successive points must not be the same, but");
+            fprintf(stderr, " X[%i]=%.3f, X[%i]=%.3f, Y[%i]=%.3f, Y[%i]=%.3f\n",
+                    i-1,X[i-1], i,X[i], i-1,Y[i-1], i,Y[i]);
+            exit(EXIT_FAILURE);
+        }
+    }
+}
+
+double sqr(double x)
+{
+    return x*x;
+}
+
+int main(int argc, char** argv)
+{
+    // default parameters
+    std::vector<double> X = {1.0, 1.0, 0.0, 0.0, 1.0, 1.0, 2.0};
+    std::vector<double> Y = {1.0, 0.0, 0.0, 1.0, 1.0, 2.0, 2.0};
+    tk::spline::spline_type type = tk::spline::cspline;
+    bool make_monotonic = false;
+
+    // override defaults with command line arguments if supplied
+    parse_args(argc, argv, X, Y, type, make_monotonic);
+
+    // setup parametric spline
+    // introduce a "time variable"
+    std::vector<double> T(X.size());
+    T[0]=0.0;
+    for(size_t i=1; i<T.size(); i++) {
+        T[i] = T[i-1] + sqrt( sqr(X[i]-X[i-1]) + sqr(Y[i]-Y[i-1]) );
+    }
+    // define a spline for each coordinate x, y
+    tk::spline sx, sy;
+    sx.set_points(T,X,type);
+    sy.set_points(T,Y,type);
+    if(make_monotonic) {
+        // adjusts spline coeffs to be piecewise monotonic where possible
+        sx.make_monotonic();
+        sy.make_monotonic();
+    }
+
+    // evaluates spline and outputs data to be used with gnuplot
+    printf("# input x, input y\n");
+    for(size_t i=0; i<X.size(); i++) {
+        printf("%f %f\n", X[i], Y[i]);  // input grid points
+    }
+    printf("\n");
+    int n = 1000;    // number of grid points to plot the spline
+    double tmin = T[0] - 0.5;
+    double tmax = T.back() + 0.5;
+    printf("# t, sx(t), sy(t)\n");
+    for(int i=0; i<n; i++) {
+        double t = tmin + (double)i*(tmax-tmin)/(n-1);
+        printf("%f %f %f\n", t, sx(t), sy(t));
+    }
+    printf("\n");
+
+    return EXIT_SUCCESS;
+}
+
+
diff --git a/examples/plot_parametric.gp b/examples/plot_parametric.gp
new file mode 100644
index 0000000..cff4cdb
--- /dev/null
+++ b/examples/plot_parametric.gp
@@ -0,0 +1,10 @@
+# gnuplot script file
+
+set size ratio -1
+set xzeroaxis
+set yzeroaxis
+set key top right
+plot    "plot.csv" every :::0::0 using 1:2 notitle with p ps 2, \
+        "plot.csv" every :::1::1 using 2:3 title "parametric spline" with l lt 1
+
+pause -1
diff --git a/examples/plot_parametric.sh b/examples/plot_parametric.sh
new file mode 100755
index 0000000..e05253a
--- /dev/null
+++ b/examples/plot_parametric.sh
@@ -0,0 +1,6 @@
+#!/bin/bash
+./plot_parametric $@ > plot.csv
+return_value=$?
+if [ $return_value -eq 0 ] ; then
+	gnuplot plot_parametric.gp
+fi
diff --git a/examples/simple_demo.cpp b/examples/simple_demo.cpp
new file mode 100644
index 0000000..5a19789
--- /dev/null
+++ b/examples/simple_demo.cpp
@@ -0,0 +1,22 @@
+#include <cstdio>
+#include <cstdlib>
+#include <vector>
+#include "spline.h"
+
+int main(int argc, char** argv)
+{
+    double x=0.0;
+    if(argc>1)
+        x = atof(argv[1]);
+
+    std::vector<double> X = {-0.1, 0.4, 1.2, 1.8, 2.0}; // must be increasing
+    std::vector<double> Y = {0.1, 0.7, 0.6, 1.1, 0.9};
+
+    tk::spline s1(X,Y); // defaults to C^2 cubic spline (spline::cspline)
+    tk::spline s2(X,Y,tk::spline::cspline_hermite);
+    tk::spline s3(X,Y,tk::spline::linear);
+
+    printf("spline(%.3f): cubic (C^2) = %.3f, cubic hermite = %.3f, linear = %.3f\n", x, s1(x), s2(x), s3(x));
+
+    return EXIT_SUCCESS;
+}
diff --git a/examples/tests/compile_units.cpp b/examples/tests/compile_units.cpp
new file mode 100644
index 0000000..a8d0ada
--- /dev/null
+++ b/examples/tests/compile_units.cpp
@@ -0,0 +1,19 @@
+// test that we can include spline.h in more than one compilation unit
+#include <cstdio>
+#include <cstdlib>
+#include <vector>
+#include "spline.h"     // not needed here but we include it anyway
+#include "myfunc.h"
+
+
+int main(int argc, char** argv)
+{
+    double x=0.0;
+    if(argc>1)
+        x = atof(argv[1]);
+
+    printf("f1(%f)=%f\n", x, myfunc1(x));
+    printf("f2(%f)=%f\n", x, myfunc2(x));
+    return EXIT_SUCCESS;
+}
+
diff --git a/examples/tests/myfunc.h b/examples/tests/myfunc.h
new file mode 100644
index 0000000..75a138b
--- /dev/null
+++ b/examples/tests/myfunc.h
@@ -0,0 +1,10 @@
+#ifndef MYFUNC_H
+#define MYFUNC_H 
+
+#include "spline.h"     // not needed here but we include it anyway
+
+double myfunc1(double x);
+double myfunc2(double x);
+
+
+#endif /* MYFUNC_H */
diff --git a/examples/tests/myfunc1.cpp b/examples/tests/myfunc1.cpp
new file mode 100644
index 0000000..c2717b9
--- /dev/null
+++ b/examples/tests/myfunc1.cpp
@@ -0,0 +1,33 @@
+// test that we can include spline.h in more than one compilation unit
+#include <cstdio>
+#include <cstdlib>
+#include <vector>
+#include "spline.h"
+
+class myclass
+{
+protected:
+    tk::spline m_spline1;
+    tk::spline m_spline2;
+public:
+    myclass(const std::vector<double>& X, const std::vector<double>& Y):
+        m_spline1(X,Y)
+    {
+        std::vector<double> X2 = {0.0, 0.5, 1.0, 2.0};
+        std::vector<double> Y2 = {0.6, 0.5, 0.3, 1.1};
+        m_spline2.set_points(X2,Y2);
+    }
+    double getval(double x) const
+    {
+        return m_spline1(x)+m_spline2(x);
+    }
+};
+
+double myfunc1(double x)
+{
+
+    static std::vector<double> X = {0.1, 0.4, 1.2, 1.8, 2.0};
+    static std::vector<double> Y = {0.1, 0.7, 0.6, 1.1, 0.9};
+    static myclass s(X,Y);
+    return s.getval(x);
+}
diff --git a/examples/tests/myfunc2.cpp b/examples/tests/myfunc2.cpp
new file mode 100644
index 0000000..2f9129b
--- /dev/null
+++ b/examples/tests/myfunc2.cpp
@@ -0,0 +1,17 @@
+// test that we can include spline.h in more than one compilation unit
+#include <cstdio>
+#include <cstdlib>
+#include <vector>
+#include "spline.h"
+
+double myfunc2(double x)
+{
+
+    static std::vector<double> X = {0.1, 0.4, 1.2, 1.8, 2.0};
+    static std::vector<double> Y = {0.1, 0.7, 0.6, 1.1, 0.9};
+    tk::spline s;
+    s.set_boundary(tk::spline::second_deriv, 0.0,
+                   tk::spline::second_deriv, 0.1);
+    s.set_points(X,Y);
+    return s(x);
+}
diff --git a/examples/tests/unit_tests.cpp b/examples/tests/unit_tests.cpp
new file mode 100644
index 0000000..b70bda0
--- /dev/null
+++ b/examples/tests/unit_tests.cpp
@@ -0,0 +1,424 @@
+#include "spline.h"
+
+#include <cstdio>
+#include <cstdlib>
+#include <cassert>
+#include <cmath>
+#include <vector>
+#include <string>
+
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE SplineUnitTests
+#include <boost/test/unit_test.hpp>
+
+
+// setup different types of splines for testing
+// which are at least "cont_deriv" times continuously differentiable
+std::vector<tk::spline> setup_splines(const std::vector<double>& X,
+                                      const std::vector<double> Y,
+                                      int cont_deriv)
+{
+    assert(X.size()==Y.size());
+    std::vector<tk::spline> splines;
+
+    if(cont_deriv<=2) {
+        // C^2 splines
+        // default cubic spline
+        splines.push_back(tk::spline());
+        splines.back().set_points(X,Y);
+        // default cubic spline: non-zero curvature at boundary
+        splines.push_back(tk::spline());
+        splines.back().set_boundary(tk::spline::second_deriv, 0.2,
+                                    tk::spline::second_deriv, -0.3);
+        splines.back().set_points(X,Y);
+    }
+    if(cont_deriv<=1) {
+        // C^1 splines
+        // default cubic spline: left derivative given
+        splines.push_back(tk::spline());
+        splines.back().set_boundary(tk::spline::first_deriv, 0.5,
+                                    tk::spline::second_deriv, 0.0);
+        splines.back().set_points(X,Y);
+        // default cubic spline: right derivative given
+        splines.push_back(tk::spline());
+        splines.back().set_boundary(tk::spline::second_deriv, 0.0,
+                                    tk::spline::first_deriv, 1.2);
+        splines.back().set_points(X,Y);
+        // default cubic spline: left and right derivative given
+        splines.push_back(tk::spline());
+        splines.back().set_boundary(tk::spline::first_deriv, 0.0,
+                                    tk::spline::first_deriv, -1.2);
+        splines.back().set_points(X,Y);
+
+        // cubic hermite spline
+        splines.push_back(tk::spline());
+        splines.back().set_points(X,Y,tk::spline::cspline_hermite);
+        // cubic hermite spline: left derivative given
+        splines.push_back(tk::spline());
+        splines.back().set_boundary(tk::spline::first_deriv, 0.5,
+                                    tk::spline::second_deriv, 0.0);
+        splines.back().set_points(X,Y,tk::spline::cspline_hermite);
+        // cubic hermite spline: right derivative given
+        splines.push_back(tk::spline());
+        splines.back().set_boundary(tk::spline::second_deriv, 0.0,
+                                    tk::spline::first_deriv, 1.2);
+        splines.back().set_points(X,Y,tk::spline::cspline_hermite);
+        // cubic hermite spline: left and right derivative given
+        splines.push_back(tk::spline());
+        splines.back().set_boundary(tk::spline::first_deriv, 0.0,
+                                    tk::spline::first_deriv, -1.2);
+        splines.back().set_points(X,Y,tk::spline::cspline_hermite);
+    }
+    if(cont_deriv<=0) {
+        // C^0 splines
+        // linear interpolation
+        splines.push_back(tk::spline());
+        splines.back().set_points(X,Y,tk::spline::linear);
+    }
+
+    // add the same splines but with monotinicity requested
+    if(cont_deriv<=1) {
+        size_t n=splines.size();
+        for(size_t i=0; i<n; i++) {
+            splines.push_back(splines[i]);
+            splines.back().make_monotonic();
+        }
+    }
+
+    return splines;
+}
+
+
+std::vector<double> evaluation_grid(const std::vector<double>& X,
+                                    double margin, size_t more_points)
+{
+    std::vector<double> eval_grid=X;
+    for(size_t i=0; i<more_points; i++) {
+        double p = X[0]-margin + (X.back()-X[0]+2.0*margin)*i/more_points;
+        eval_grid.push_back(p);
+    }
+    return eval_grid;
+}
+
+
+// checks spline is C^0, i.e. continuous and goes through all grid points
+BOOST_AUTO_TEST_CASE( Continuity )
+{
+    const double dx = 1e-15;        // step size to check continuity
+    const double max_deriv = 10.0;  // max 1st deriv of below spline
+    // using well posed grid data for spline interpolation
+    std::vector<double> X = {-10.2,  0.1,  0.8, 2.4,  3.1,  3.2, 4.7, 19.1};
+    std::vector<double> Y = {  2.7, -1.2, -0.5, 1.5, -1.0, -0.7,  1.2, -1.3};
+    // points to check continuity of spline
+    std::vector<double> eval_grid = evaluation_grid(X,1.0,1000);
+
+    // setup all possible types of splines which are at least C^0
+    std::vector<tk::spline> splines=setup_splines(X,Y,0);
+
+    for(size_t k=0; k<splines.size(); k++) {
+        const tk::spline& s=splines[k];
+        BOOST_TEST_CONTEXT("spline " << k << ": " << s.info())
+        {
+        for(size_t i=0; i<eval_grid.size(); i++) {
+            double x = eval_grid[i];
+            if(i<Y.size()) {
+                double y = Y[i];
+                BOOST_CHECK_CLOSE( s(x), y, 0.0 );          // hit grid points
+            }
+            BOOST_CHECK_SMALL(s(x)-s(x+dx), max_deriv*dx);  // left continuous
+            BOOST_CHECK_SMALL(s(x)-s(x-dx), max_deriv*dx);  // right continous
+        }
+        } // BOOST_TEST_CONTEXT
+    }
+}
+
+// checks spline is C^1, i.e. differentiable
+BOOST_AUTO_TEST_CASE( Differentiability )
+{
+    const double dx = 1e-15;        // step size to check continuity
+    const double h  = 2e-8;         // step size for finite differences
+    const double max_deriv = 100.0; // max 2nd deriv of below spline
+    std::vector<double> X = {-10.2,  0.1,  0.8, 2.4,  3.1,  3.2, 4.7, 19.1};
+    std::vector<double> Y = {  2.7, -1.2, -0.5, 1.5, -1.0, -0.7,  1.2, -1.3};
+    // points to check continuity/differentiability of spline
+    std::vector<double> eval_grid = evaluation_grid(X,1.0,1000);
+
+    // setup all possible types of splines which are at least C^1
+    std::vector<tk::spline> splines=setup_splines(X,Y,1);
+
+    for(size_t k=0; k<splines.size(); k++) {
+        const tk::spline& s=splines[k];
+        BOOST_TEST_CONTEXT("spline " << k << ": " << s.info())
+        {
+        for(size_t i=0; i<eval_grid.size(); i++) {
+            double x = eval_grid[i];
+            double right_diff =   (s(x+h)-s(x)) / h;
+            double left_diff =    (s(x)-s(x-h)) / h;
+            double central_diff = (s(x+h)-s(x-h) ) / (2.0*h);
+    
+            BOOST_CHECK_SMALL( left_diff-right_diff, max_deriv*h );
+            BOOST_CHECK_SMALL( central_diff-right_diff, max_deriv*h);
+            BOOST_CHECK_SMALL( s.deriv(1,x)-central_diff, max_deriv*h );
+    
+            BOOST_CHECK_SMALL( s.deriv(1,x-dx)-s.deriv(1,x), max_deriv*dx );
+            BOOST_CHECK_SMALL( s.deriv(1,x+dx)-s.deriv(1,x), max_deriv*dx );
+        }
+        } // BOOST_TEST_CONTEXT
+    }
+}
+
+
+// checks spline is C^2
+BOOST_AUTO_TEST_CASE( Smoothness )
+{
+    const double dx = 1e-15;        // step size to check continuity
+    const double h  = 3e-6;         // step size for finite differences
+    const double max_deriv = 1000.0;    // max 3nd deriv of below spline
+    std::vector<double> X = {-10.2,  0.1,  0.8, 2.4,  3.1,  3.2, 4.7, 19.1};
+    std::vector<double> Y = {  2.7, -1.2, -0.5, 1.5, -1.0, -0.7,  1.2, -1.3};
+    // points to check continuity/differentiability of spline
+    std::vector<double> eval_grid = evaluation_grid(X,1.0,1000);
+
+    // setup all possible types of splines which are at least C^2
+    std::vector<tk::spline> splines=setup_splines(X,Y,2);
+
+    for(size_t k=0; k<splines.size(); k++) {
+        const tk::spline& s=splines[k];
+        BOOST_TEST_CONTEXT("spline " << k << ": " << s.info())
+        {
+        for(size_t i=0; i<eval_grid.size(); i++) {
+            double x = eval_grid[i];
+            double right_diff =   (s(x)-2.0*s(x+h)+s(x+2.0*h)) / (h*h);
+            double left_diff =    (s(x-2.0*h)-2.0*s(x-h)+s(x)) / (h*h);
+            double central_diff = (s(x-h)-2.0*s(x)+s(x+h)) / (h*h);
+    
+            BOOST_CHECK_SMALL( left_diff-right_diff, max_deriv*h );
+            BOOST_CHECK_SMALL( central_diff-right_diff, max_deriv*h );
+            BOOST_CHECK_SMALL( s.deriv(2,x)-central_diff, max_deriv*h );
+    
+            BOOST_CHECK_SMALL( s.deriv(2,x-dx)-s.deriv(2,x), max_deriv*dx );
+            BOOST_CHECK_SMALL( s.deriv(2,x+dx)-s.deriv(2,x), max_deriv*dx );
+        }
+        } // BOOST_TEST_CONTEXT
+    }
+}
+
+// checks the calculation of the third derivative
+BOOST_AUTO_TEST_CASE( ThirdDeriv )
+{
+    // 2nd derivative is piecewise linear, so we can use a larger
+    // step size for finite differences for 3rd derivative
+    const double h    = 1e-5;     // step size for finite differences
+    const double tol  = 1e-10;    // error tolerance in 3rd derivative
+    std::vector<double> X = {-10.2,  0.1,  0.8, 2.4,  3.1,  3.2, 4.7, 19.1};
+    std::vector<double> Y = {  2.7, -1.2, -0.5, 1.5, -1.0, -0.7,  1.2, -1.3};
+
+    // evaluation grid must not include any of the grid points X
+    std::vector<double> eval_grid = { -12.3, -4.3, 0.5, 3.5, 7.3, 22.1, 95.0};
+
+    // setup all possible types of splines which are at least C^0
+    std::vector<tk::spline> splines=setup_splines(X,Y,0);
+
+    for(size_t k=0; k<splines.size(); k++) {
+        const tk::spline& s=splines[k];
+        BOOST_TEST_CONTEXT("spline " << k << ": " << s.info())
+        {
+        for(size_t i=0; i<eval_grid.size(); i++) {
+            double x = eval_grid[i];
+            double right_diff =   (s.deriv(2,x+h)-s.deriv(2,x)) / h;
+            double left_diff =    (s.deriv(2,x)-s.deriv(2,x-h)) / h;
+            double central_diff = (s.deriv(2,x+h)-s.deriv(2,x-h) ) / (2.0*h);
+    
+            BOOST_CHECK_SMALL( left_diff-right_diff, tol );
+            BOOST_CHECK_SMALL( central_diff-right_diff, tol );
+            BOOST_CHECK_SMALL( s.deriv(3,x)-central_diff, tol );
+    
+            BOOST_CHECK_SMALL( s.deriv(3,x-h)-s.deriv(3,x), 0.0 );  // piecew. const
+            BOOST_CHECK_SMALL( s.deriv(3,x+h)-s.deriv(3,x), 0.0 );  // piecew. const
+        }
+        } // BOOST_TEST_CONTEXT
+    }
+}
+
+
+BOOST_AUTO_TEST_CASE( BoundaryTypes )
+{
+    const double dx = 1e-5;        // step size to check continuity
+    std::vector<double> X = {-10.2,  0.1,  0.8, 2.4,  3.1,  3.2, 4.7, 19.1};
+    std::vector<double> Y = {  2.7, -1.2, -0.5, 1.5, -1.0, -0.7,  1.2, -1.3};
+    {
+        // natural boundary conditions, i.e. 2nd derivative is zero
+        tk::spline s;
+        s.set_points(X,Y);
+        BOOST_CHECK_CLOSE( s.deriv(2,X[0]-dx), 0.0, 0.0 );      // 2nd deriv = 0
+        BOOST_CHECK_CLOSE( s.deriv(2,X.back()+dx), 0.0, 0.0 );  // 2nd deriv = 0
+    }
+
+    {
+        // left and right with given 2nd and 1st derivative, respectively
+        tk::spline s;
+        double deriv1=1.3;
+        double deriv2=-0.7;
+        s.set_boundary(tk::spline::second_deriv, deriv2,
+                       tk::spline::first_deriv, deriv1);
+        s.set_points(X,Y);
+        BOOST_CHECK_CLOSE( s.deriv(2,X[0]), deriv2, 0.0 );
+        BOOST_CHECK_CLOSE( s.deriv(2,X[0]-1.0), deriv2, 0.0 );
+        BOOST_CHECK_CLOSE( s.deriv(1,X.back()), deriv1, 1e-12);
+        BOOST_CHECK_CLOSE( s.deriv(1,X.back()+1.0), deriv1, 1e-12);
+        BOOST_CHECK_CLOSE( s.deriv(2,X.back()+1.0), 0.0, 0.0);
+    }
+    {
+        // left and right with given 1st and 2nd derivative, respectively
+        tk::spline s;
+        double deriv1=-4.7;
+        double deriv2=-1.2;
+        s.set_boundary(tk::spline::first_deriv, deriv1,
+                       tk::spline::second_deriv, deriv2);
+        s.set_points(X,Y);
+        BOOST_CHECK_CLOSE( s.deriv(2,X.back()), deriv2, 0.0 );
+        BOOST_CHECK_CLOSE( s.deriv(2,X.back()+1.0), deriv2, 0.0 );
+        BOOST_CHECK_CLOSE( s.deriv(1,X[0]), deriv1, 1e-12);
+        BOOST_CHECK_CLOSE( s.deriv(1,X[0]-1.0), deriv1, 1e-12);
+        BOOST_CHECK_CLOSE( s.deriv(2,X[0]-1.0), 0.0, 0.0);
+    }
+}
+
+
+// function to be approximated by a spline
+double myfunc(double x)
+{
+    return sin(x);
+}
+
+// this test checks how well a smooth function is approximated by a spline
+BOOST_AUTO_TEST_CASE( FunctionApproximation )
+{
+    double lo=0.0;
+    double hi=2.0*M_PI;
+    // regression results
+    std::vector< std::vector<double> > max_errors = {
+        {
+            1.0, 1.087710e-01, 2.001701e-02, 6.866929e-04,
+            3.189426e-05, 1.764343e-06, 1.043467e-07, 6.352541e-09,
+            3.919548e-10, 2.432621e-11
+        },
+        {
+            1.0, 1.895384e-01, 9.289323e-02, 7.816602e-03,
+            1.082736e-03, 1.319855e-04, 1.607416e-05, 1.977196e-06,
+            2.450258e-07, 3.049079e-08
+        }
+    };
+    std::vector< std::vector<double> > avg_errors = {
+        {
+            6.365561e-01, 4.767218e-02, 1.161861e-02, 2.479389e-04,
+            1.100872e-05, 6.015310e-07, 3.546300e-08, 2.157237e-09,
+            1.330842e-10, 8.264881e-12
+        },
+        {
+            6.365561e-01, 9.277315e-02, 3.244986e-02, 3.617337e-03,
+            3.095662e-04, 3.116561e-05, 3.514112e-06, 4.189284e-07,
+            5.122710e-08, 6.337061e-09
+        }
+    };
+    double tol=1e-6;        // 6 correct digits
+
+    // evaluation grid
+    std::vector<double> eval_grid(10000);
+    for(size_t i=0; i<eval_grid.size(); i++)
+        eval_grid[i] = lo + (double)i/(eval_grid.size()-1) * (hi-lo);
+
+    // spline sizes
+    std::vector<size_t> N = {3, 4, 5, 10, 20, 40, 80, 160, 320, 640};
+
+    // spline types
+    std::vector<tk::spline::spline_type> types =
+    {tk::spline::cspline, tk::spline::cspline_hermite};
+    std::vector<std::string> type_names = {"cspline", "hermite"};
+    for(size_t l=0; l<2; l++) {
+        tk::spline::spline_type type=types[l];
+        const std::vector<double>& avg_error = avg_errors[l];
+        const std::vector<double>& max_error = max_errors[l];
+        const std::string type_name = type_names[l];
+
+        for(size_t k=0; k<N.size(); k++) {
+            size_t n=N[k];
+            std::vector<double> X(n), Y(n);
+            for(size_t i=0; i<n; i++) {
+                X[i] = lo + (double)i/(n-1) * (hi-lo);
+                Y[i] = myfunc(X[i]);
+            }
+            tk::spline s;
+            s.set_points(X,Y,type);
+            for(size_t i=0; i<n; i++) {
+                X[i] = lo + (double)i/(n-1) * (hi-lo);
+                Y[i] = myfunc(X[i]);
+            }
+            double avg = 0.0;
+            double max = 0.0;
+            for(size_t i=0; i<eval_grid.size(); i++) {
+                double x = eval_grid[i];
+                double err = std::fabs(myfunc(x)-s(x));
+                avg += err;
+                if(err>max)
+                    max=err;
+            }
+            avg /= eval_grid.size();
+            printf("approx %s: n = %3lu, avg error = %e, max error = %e\n",
+                   type_name.c_str(), n, avg, max );
+            BOOST_CHECK_CLOSE(avg, avg_error[k], tol*100.0);
+            BOOST_CHECK_CLOSE(max, max_error[k], tol*100.0);
+        }
+    }
+}
+
+double rand_unif()
+{
+    return (double)std::rand()/RAND_MAX;
+}
+
+// this test is intended to detect any numerical changes
+// and will require re-basing whenever a numerical change is expected
+BOOST_AUTO_TEST_CASE( RegressionTest )
+{
+    std::vector<size_t> N = {3, 4, 5, 9, 13, 22, 86, 137, 657, 12357};
+
+    double sum=0.0;
+    double sum_mod=0.0;
+    for(size_t k=0; k<N.size(); k++) {
+        size_t n=N[k];
+        std::vector<double> X(n), Y(n);
+        X[0]=rand_unif();
+        Y[0]=rand_unif();
+        for(size_t i=1; i<n; i++) {
+            X[i] = X[i-1] + rand_unif();            // monotonic increasing
+            Y[i] = Y[i-1] + rand_unif()-0.5;
+        }
+
+        // evaluation grid
+        std::vector<double> eval_grid(1000);
+        for(size_t i=0; i<eval_grid.size(); i++)
+            eval_grid[i] = X[0]-5.0 + rand_unif()*(X.back()-X[0]+10.0);
+        eval_grid.insert(eval_grid.end(), X.begin(), X.end());  // add X points
+
+        // setup all possible types of splines which are at least C^0
+        std::vector<tk::spline> splines=setup_splines(X,Y,0);
+
+        for(size_t l=0; l<splines.size(); l++) {
+            const tk::spline& s = splines[l];
+            for(size_t i=0; i<eval_grid.size(); i++) {
+                double x = eval_grid[i];
+                double val = s(x)+0.1*s.deriv(1,x)
+                             + 0.01*s.deriv(2,x)+0.001*s.deriv(3,x);
+                sum += val;
+                sum_mod += std::fmod(val, 1.0);
+                sum_mod = std::fmod(sum_mod, 1.0);
+            }
+            //printf("spline %lu size %5lu: sum=%.16e, mod=%.16f\n",l,n,sum,sum_mod);
+        }
+    }
+    //printf("sum=%.16e, mod=%.16f\n", sum, sum_mod);
+    BOOST_CHECK_CLOSE(sum, -302474251493.73535, 0.0);
+    BOOST_CHECK_CLOSE(sum_mod, 0.25791635217581654, 0.0);
+}