1 files changed, 116 insertions, 0 deletions
diff --git a/gsl-1.9/doc/examples/bspline.c b/gsl-1.9/doc/examples/bspline.c
new file mode 100644
index 0000000..0d6a584
--- /dev/null
+++ b/gsl-1.9/doc/examples/bspline.c
@@ -0,0 +1,116 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <gsl/gsl_bspline.h>
+#include <gsl/gsl_multifit.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
+
+/* number of data points to fit */
+#define N        200
+
+/* number of fit coefficients */
+#define NCOEFFS  8
+
+/* nbreak = ncoeffs + 2 - k = ncoeffs - 2 since k = 4 */
+#define NBREAK   (NCOEFFS - 2)
+
+int
+main (void)
+{
+  const size_t n = N;
+  const size_t ncoeffs = NCOEFFS;
+  const size_t nbreak = NBREAK;
+  size_t i, j;
+  gsl_bspline_workspace *bw;
+  gsl_vector *B;
+  double dy;
+  gsl_rng *r;
+  gsl_vector *c, *w;
+  gsl_vector *x, *y;
+  gsl_matrix *X, *cov;
+  gsl_multifit_linear_workspace *mw;
+  double chisq;
+
+  gsl_rng_env_setup();
+  r = gsl_rng_alloc(gsl_rng_default);
+
+  /* allocate a cubic bspline workspace (k = 4) */
+  bw = gsl_bspline_alloc(4, nbreak);
+  B = gsl_vector_alloc(ncoeffs);
+
+  x = gsl_vector_alloc(n);
+  y = gsl_vector_alloc(n);
+  X = gsl_matrix_alloc(n, ncoeffs);
+  c = gsl_vector_alloc(ncoeffs);
+  w = gsl_vector_alloc(n);
+  cov = gsl_matrix_alloc(ncoeffs, ncoeffs);
+  mw = gsl_multifit_linear_alloc(n, ncoeffs);
+
+  printf("#m=0,S=0\n");
+  /* this is the data to be fitted */
+  for (i = 0; i < n; ++i)
+    {
+      double sigma;
+      double xi = (15.0/(N-1)) * i;
+      double yi = cos(xi) * exp(-0.1 * xi);
+
+      sigma = 0.1;
+      dy = gsl_ran_gaussian(r, sigma);
+      yi += dy;
+
+      gsl_vector_set(x, i, xi);
+      gsl_vector_set(y, i, yi);
+      gsl_vector_set(w, i, 1.0 / (sigma*sigma));
+
+      printf("%f %f\n", xi, yi);
+    }
+
+  /* use uniform breakpoints on [0, 15] */
+  gsl_bspline_knots_uniform(0.0, 15.0, bw);
+
+  /* construct the fit matrix X */
+  for (i = 0; i < n; ++i)
+    {
+      double xi = gsl_vector_get(x, i);
+
+      /* compute B_j(xi) for all j */
+      gsl_bspline_eval(xi, B, bw);
+
+      /* fill in row i of X */
+      for (j = 0; j < ncoeffs; ++j)
+        {
+          double Bj = gsl_vector_get(B, j);
+          gsl_matrix_set(X, i, j, Bj);
+        }
+    }
+
+  /* do the fit */
+  gsl_multifit_wlinear(X, w, y, c, cov, &chisq, mw);
+
+  /* output the smoothed curve */
+  {
+    double xi, yi, yerr;
+
+    printf("#m=1,S=0\n");
+    for (xi = 0.0; xi < 15.0; xi += 0.1)
+      {
+        gsl_bspline_eval(xi, B, bw);
+        gsl_multifit_linear_est(B, c, cov, &yi, &yerr);
+        printf("%f %f\n", xi, yi);
+      }
+  }
+
+  gsl_rng_free(r);
+  gsl_bspline_free(bw);
+  gsl_vector_free(B);
+  gsl_vector_free(x);
+  gsl_vector_free(y);
+  gsl_matrix_free(X);
+  gsl_vector_free(c);
+  gsl_vector_free(w);
+  gsl_matrix_free(cov);
+  gsl_multifit_linear_free(mw);
+
+  return 0;
+} /* main() */