1 files changed, 173 insertions, 0 deletions
diff --git a/gsl-1.9/interpolation/poly.c b/gsl-1.9/interpolation/poly.c
new file mode 100644
index 0000000..a0cd68a
--- /dev/null
+++ b/gsl-1.9/interpolation/poly.c
@@ -0,0 +1,173 @@
+/* interpolation/interp_poly.c
+ * 
+ * Copyright (C) 2001 DAN, HO-JIN
+ * 
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or (at
+ * your option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ */
+
+#include <config.h>
+#include <math.h>
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_poly.h>
+#include <gsl/gsl_interp.h>
+
+typedef struct
+{
+  double *d;
+  double *coeff;
+  double *work;
+}
+polynomial_state_t;
+
+static void *
+polynomial_alloc (size_t size)
+{
+  polynomial_state_t *state =
+    (polynomial_state_t *) malloc (sizeof (polynomial_state_t));
+
+  if (state == 0)
+    {
+      GSL_ERROR_NULL ("failed to allocate space for polynomial state",
+                      GSL_ENOMEM);
+    }
+
+  state->d = (double *) malloc (sizeof (double) * size);
+
+  if (state->d == 0)
+    {
+      free (state);
+      GSL_ERROR_NULL ("failed to allocate space for d", GSL_ENOMEM);
+    }
+
+  state->coeff = (double *) malloc (sizeof (double) * size);
+
+  if (state->coeff == 0)
+    {
+      free (state->d);
+      free (state);
+      GSL_ERROR_NULL ("failed to allocate space for d", GSL_ENOMEM);
+    }
+
+  state->work = (double *) malloc (sizeof (double) * size);
+
+  if (state->work == 0)
+    {
+      free (state->coeff);
+      free (state->d);
+      free (state);
+      GSL_ERROR_NULL ("failed to allocate space for d", GSL_ENOMEM);
+    }
+
+  return state;
+}
+
+static int
+polynomial_init (void *vstate,
+                 const double xa[], const double ya[], size_t size)
+{
+  polynomial_state_t *state = (polynomial_state_t *) vstate;
+
+  int status = gsl_poly_dd_init (state->d, xa, ya, size);
+
+  return status;
+}
+
+static int
+polynomial_eval (const void *vstate,
+                 const double xa[], const double ya[], size_t size, double x,
+                 gsl_interp_accel * acc, double *y)
+{
+  const polynomial_state_t *state = (const polynomial_state_t *) vstate;
+
+  *y = gsl_poly_dd_eval (state->d, xa, size, x);
+
+  return GSL_SUCCESS;
+}
+
+
+static int
+polynomial_deriv (const void *vstate,
+                  const double xa[], const double ya[], size_t size, double x,
+                  gsl_interp_accel * acc, double *y)
+{
+  const polynomial_state_t *state = (const polynomial_state_t *) vstate;
+
+  gsl_poly_dd_taylor (state->coeff, x, state->d, xa, size, state->work);
+
+  *y = state->coeff[1];
+
+  return GSL_SUCCESS;
+}
+
+static int
+polynomial_deriv2 (const void *vstate,
+                   const double xa[], const double ya[], size_t size,
+                   double x, gsl_interp_accel * acc, double *y)
+{
+  const polynomial_state_t *state = (const polynomial_state_t *) vstate;
+
+  gsl_poly_dd_taylor (state->coeff, x, state->d, xa, size, state->work);
+
+  *y = 2.0 * state->coeff[2];
+
+  return GSL_SUCCESS;
+}
+
+static int
+polynomial_integ (const void *vstate, const double xa[], const double ya[],
+                  size_t size, gsl_interp_accel * acc, double a, double b,
+                  double *result)
+{
+  const polynomial_state_t *state = (const polynomial_state_t *) vstate;
+  size_t i;
+  double sum;
+
+  gsl_poly_dd_taylor (state->coeff, 0.0, state->d, xa, size, state->work);
+
+  sum = state->coeff[0] * (b - a);
+
+  for (i = 1; i < size; i++)
+    {
+      sum += state->coeff[i] * (pow (b, i + 1) - pow (a, i + 1)) / (i + 1.0);
+    }
+
+  *result = sum;
+
+  return GSL_SUCCESS;
+}
+
+static void
+polynomial_free (void *vstate)
+{
+  polynomial_state_t *state = (polynomial_state_t *) vstate;
+  free (state->d);
+  free (state->coeff);
+  free (state->work);
+  free (state);
+}
+
+static const gsl_interp_type polynomial_type = {
+  "polynomial",
+  3,
+  &polynomial_alloc,
+  &polynomial_init,
+  &polynomial_eval,
+  &polynomial_deriv,
+  &polynomial_deriv2,
+  &polynomial_integ,
+  &polynomial_free,
+};
+
+const gsl_interp_type *gsl_interp_polynomial = &polynomial_type;