diff --git a/ext/f2c_recipes/.cvsignore b/ext/f2c_recipes/.cvsignore
new file mode 100644
index 000000000..756fc5489
--- /dev/null
+++ b/ext/f2c_recipes/.cvsignore
@@ -0,0 +1,3 @@
+Makefile
+*.d
+.depends
diff --git a/ext/f2c_recipes/Makefile.in b/ext/f2c_recipes/Makefile.in
new file mode 100644
index 000000000..2ec83d9c2
--- /dev/null
+++ b/ext/f2c_recipes/Makefile.in
@@ -0,0 +1,75 @@
+#/bin/sh
+#
+# $Source$
+# $Author$
+# $Revision$
+# $Date$
+#
+
+.SUFFIXES :
+.SUFFIXES : .c .d .o
+
+# the directory where the Cantera libraries are located
+CANTERA_LIBDIR=@buildlib@
+
+# the directory where Cantera include files may be found.
+CANTERA_INCDIR=@ctroot@/build/include/cantera
+
+# the C++ compiler
+CXX = @CXX@
+
+# the C compiler
+CC = @CC@
+
+# C++ compile flags
+CXX_FLAGS = @CXXFLAGS@ $(CXX_OPT)
+
+# Local include files
+CXX_INCLUDES=-I../f2c_libs
+
+# How to compile the dependency file
+.c.d:
+	g++ -MM  $(CXX_FLAGS) $(CXX_INCLUDES) $*.c > $*.d
+
+# How to compile a C file
+.c.o:
+	@CC@ -c $< @DEFS@ $(CXX_FLAGS) $(CXX_INCLUDES)
+
+# -----------------------------------------------
+
+LIB = $(CANTERA_LIBDIR)/librecipes.a
+
+all:  $(LIB)
+
+# list of object files
+OBJS = simp1.o simp2.o simp3.o simplx.o splint.o \
+       splie2.o spline.o splin2.o
+
+# list of source files
+SRCS = $(OBJS:.o=.c)
+
+# List of dependency files to be created
+DEPENDS=$(OBJS:.o=.d)
+
+# rule to make library
+$(LIB): $(OBJS)
+	@ARCHIVE@ $(LIB) $(OBJS) > /dev/null
+
+# ------------------------------------------------
+# Utility Targets
+
+clean:
+	$(RM) $(OBJS) $(LIB) *.d .depends
+
+# depends target
+depends:
+	$(RM) *.d .depends
+	@MAKE@ .depends
+
+.depends: $(DEPENDS)
+	cat *.d > .depends
+
+
+ifeq ($(wildcard .depends), .depends)
+include .depends
+endif
diff --git a/ext/f2c_recipes/simp1.c b/ext/f2c_recipes/simp1.c
new file mode 100644
index 000000000..1bd8421dc
--- /dev/null
+++ b/ext/f2c_recipes/simp1.c
@@ -0,0 +1,61 @@
+/* simp1.f -- translated by f2c (version 20031025).
+   You must link the resulting object file with libf2c:
+	on Microsoft Windows system, link with libf2c.lib;
+	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
+	or, if you install libf2c.a in a standard place, with -lf2c -lm
+	-- in that order, at the end of the command line, as in
+		cc *.o -lf2c -lm
+	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
+
+		http://www.netlib.org/f2c/libf2c.zip
+*/
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+#include "f2c.h"
+
+/* Subroutine */ int simp1_(doublereal *a, integer *mp, integer *np, integer *
+	mm, integer *ll, integer *nll, integer *iabf, integer *kp, doublereal 
+	*bmax)
+{
+    /* System generated locals */
+    integer a_dim1, a_offset, i__1;
+    doublereal d__1;
+
+    /* Local variables */
+    static integer k;
+    static doublereal test;
+
+    /* Parameter adjustments */
+    --ll;
+    a_dim1 = *mp;
+    a_offset = 1 + a_dim1;
+    a -= a_offset;
+
+    /* Function Body */
+    *kp = ll[1];
+    *bmax = a[*mm + 1 + (*kp + 1) * a_dim1];
+    if (*nll < 2) {
+	return 0;
+    }
+    i__1 = *nll;
+    for (k = 2; k <= i__1; ++k) {
+	if (*iabf == 0) {
+	    test = a[*mm + 1 + (ll[k] + 1) * a_dim1] - *bmax;
+	} else {
+	    test = (d__1 = a[*mm + 1 + (ll[k] + 1) * a_dim1], abs(d__1)) - 
+		    abs(*bmax);
+	}
+	if (test > (float)0.) {
+	    *bmax = a[*mm + 1 + (ll[k] + 1) * a_dim1];
+	    *kp = ll[k];
+	}
+/* L11: */
+    }
+    return 0;
+} /* simp1_ */
+
+#ifdef __cplusplus
+	}
+#endif
diff --git a/ext/f2c_recipes/simp2.c b/ext/f2c_recipes/simp2.c
new file mode 100644
index 000000000..f4130a29d
--- /dev/null
+++ b/ext/f2c_recipes/simp2.c
@@ -0,0 +1,89 @@
+/* simp2.f -- translated by f2c (version 20031025).
+   You must link the resulting object file with libf2c:
+	on Microsoft Windows system, link with libf2c.lib;
+	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
+	or, if you install libf2c.a in a standard place, with -lf2c -lm
+	-- in that order, at the end of the command line, as in
+		cc *.o -lf2c -lm
+	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
+
+		http://www.netlib.org/f2c/libf2c.zip
+*/
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+#include "f2c.h"
+
+/* Subroutine */ int simp2_(doublereal *a, integer *m, integer *n, integer *
+	mp, integer *np, integer *l2, integer *nl2, integer *ip, integer *kp, 
+	doublereal *q1)
+{
+    /* System generated locals */
+    integer a_dim1, a_offset, i__1, i__2;
+
+    /* Local variables */
+    static integer i__, k;
+    static doublereal q, q0;
+    static integer ii;
+    static doublereal qp;
+
+    /* Parameter adjustments */
+    --l2;
+    a_dim1 = *mp;
+    a_offset = 1 + a_dim1;
+    a -= a_offset;
+
+    /* Function Body */
+    *ip = 0;
+    if (*nl2 < 1) {
+	return 0;
+    }
+    i__1 = *nl2;
+    for (i__ = 1; i__ <= i__1; ++i__) {
+	if (a[l2[i__] + 1 + (*kp + 1) * a_dim1] < -1e-6) {
+	    goto L2;
+	}
+/* L11: */
+    }
+    return 0;
+L2:
+    *q1 = -a[l2[i__] + 1 + a_dim1] / a[l2[i__] + 1 + (*kp + 1) * a_dim1];
+    *ip = l2[i__];
+    if (i__ + 1 > *nl2) {
+	return 0;
+    }
+    i__1 = *nl2;
+    for (++i__; i__ <= i__1; ++i__) {
+	ii = l2[i__];
+	if (a[ii + 1 + (*kp + 1) * a_dim1] < -1e-6) {
+	    q = -a[ii + 1 + a_dim1] / a[ii + 1 + (*kp + 1) * a_dim1];
+	    if (q < *q1) {
+		*ip = ii;
+		*q1 = q;
+	    } else if (q == *q1) {
+		i__2 = *n;
+		for (k = 1; k <= i__2; ++k) {
+		    qp = -a[*ip + 1 + (k + 1) * a_dim1] / a[*ip + 1 + (*kp + 
+			    1) * a_dim1];
+		    q0 = -a[ii + 1 + (k + 1) * a_dim1] / a[ii + 1 + (*kp + 1) 
+			    * a_dim1];
+		    if (q0 != qp) {
+			goto L6;
+		    }
+/* L12: */
+		}
+L6:
+		if (q0 < qp) {
+		    *ip = ii;
+		}
+	    }
+	}
+/* L13: */
+    }
+    return 0;
+} /* simp2_ */
+
+#ifdef __cplusplus
+	}
+#endif
diff --git a/ext/f2c_recipes/simp3.c b/ext/f2c_recipes/simp3.c
new file mode 100644
index 000000000..366169f7f
--- /dev/null
+++ b/ext/f2c_recipes/simp3.c
@@ -0,0 +1,65 @@
+/* simp3.f -- translated by f2c (version 20031025).
+   You must link the resulting object file with libf2c:
+	on Microsoft Windows system, link with libf2c.lib;
+	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
+	or, if you install libf2c.a in a standard place, with -lf2c -lm
+	-- in that order, at the end of the command line, as in
+		cc *.o -lf2c -lm
+	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
+
+		http://www.netlib.org/f2c/libf2c.zip
+*/
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+#include "f2c.h"
+
+/* Subroutine */ int simp3_(doublereal *a, integer *mp, integer *np, integer *
+	i1, integer *k1, integer *ip, integer *kp)
+{
+    /* System generated locals */
+    integer a_dim1, a_offset, i__1, i__2;
+
+    /* Local variables */
+    static integer ii, kk;
+    static doublereal piv;
+
+    /* Parameter adjustments */
+    a_dim1 = *mp;
+    a_offset = 1 + a_dim1;
+    a -= a_offset;
+
+    /* Function Body */
+    piv = 1. / a[*ip + 1 + (*kp + 1) * a_dim1];
+    if (*i1 >= 0) {
+	i__1 = *i1 + 1;
+	for (ii = 1; ii <= i__1; ++ii) {
+	    if (ii - 1 != *ip) {
+		a[ii + (*kp + 1) * a_dim1] *= piv;
+		i__2 = *k1 + 1;
+		for (kk = 1; kk <= i__2; ++kk) {
+		    if (kk - 1 != *kp) {
+			a[ii + kk * a_dim1] -= a[*ip + 1 + kk * a_dim1] * a[
+				ii + (*kp + 1) * a_dim1];
+		    }
+/* L11: */
+		}
+	    }
+/* L12: */
+	}
+    }
+    i__1 = *k1 + 1;
+    for (kk = 1; kk <= i__1; ++kk) {
+	if (kk - 1 != *kp) {
+	    a[*ip + 1 + kk * a_dim1] = -a[*ip + 1 + kk * a_dim1] * piv;
+	}
+/* L13: */
+    }
+    a[*ip + 1 + (*kp + 1) * a_dim1] = piv;
+    return 0;
+} /* simp3_ */
+
+#ifdef __cplusplus
+	}
+#endif
diff --git a/ext/f2c_recipes/simplx.c b/ext/f2c_recipes/simplx.c
new file mode 100644
index 000000000..e885f42ee
--- /dev/null
+++ b/ext/f2c_recipes/simplx.c
@@ -0,0 +1,197 @@
+/* simplx.f -- translated by f2c (version 20031025).
+   You must link the resulting object file with libf2c:
+	on Microsoft Windows system, link with libf2c.lib;
+	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
+	or, if you install libf2c.a in a standard place, with -lf2c -lm
+	-- in that order, at the end of the command line, as in
+		cc *.o -lf2c -lm
+	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
+
+		http://www.netlib.org/f2c/libf2c.zip
+*/
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+#include "f2c.h"
+
+/* Table of constant values */
+
+static integer c__0 = 0;
+static integer c__1 = 1;
+
+/* Subroutine */ int simplx_(doublereal *a, integer *m, integer *n, integer *
+	mp, integer *np, integer *m1, integer *m2, integer *m3, integer *
+	icase, integer *izrov, integer *iposv)
+{
+    /* System generated locals */
+    integer a_dim1, a_offset, i__1, i__2;
+
+    /* Builtin functions */
+    /* Subroutine */ int s_paus(char *, ftnlen);
+
+    /* Local variables */
+    static integer i__, k, l1[1000], l2[1000], l3[1000];
+    static doublereal q1;
+    static integer m12, kh, ip, ir, kp, is, nl1, nl2;
+    static doublereal bmax;
+    extern /* Subroutine */ int simp1_(doublereal *, integer *, integer *, 
+	    integer *, integer *, integer *, integer *, integer *, doublereal 
+	    *), simp2_(doublereal *, integer *, integer *, integer *, integer 
+	    *, integer *, integer *, integer *, integer *, doublereal *), 
+	    simp3_(doublereal *, integer *, integer *, integer *, integer *, 
+	    integer *, integer *);
+
+    /* Parameter adjustments */
+    --iposv;
+    --izrov;
+    a_dim1 = *mp;
+    a_offset = 1 + a_dim1;
+    a -= a_offset;
+
+    /* Function Body */
+    if (*m != *m1 + *m2 + *m3) {
+	s_paus("Bad input constraint counts.", (ftnlen)28);
+    }
+    nl1 = *n;
+    i__1 = *n;
+    for (k = 1; k <= i__1; ++k) {
+	l1[k - 1] = k;
+	izrov[k] = k;
+/* L11: */
+    }
+    nl2 = *m;
+    i__1 = *m;
+    for (i__ = 1; i__ <= i__1; ++i__) {
+	if (a[i__ + 1 + a_dim1] < 0.) {
+/*          write(*,*) 'The A matrix input to SIMPLX is invalid.' */
+	    s_paus("Bad input tableau.", (ftnlen)18);
+	}
+	l2[i__ - 1] = i__;
+	iposv[i__] = *n + i__;
+/* L12: */
+    }
+    i__1 = *m2;
+    for (i__ = 1; i__ <= i__1; ++i__) {
+	l3[i__ - 1] = 1;
+/* L13: */
+    }
+    ir = 0;
+    if (*m2 + *m3 == 0) {
+	goto L30;
+    }
+    ir = 1;
+    i__1 = *n + 1;
+    for (k = 1; k <= i__1; ++k) {
+	q1 = (float)0.;
+	i__2 = *m;
+	for (i__ = *m1 + 1; i__ <= i__2; ++i__) {
+	    q1 += a[i__ + 1 + k * a_dim1];
+/* L14: */
+	}
+	a[*m + 2 + k * a_dim1] = -q1;
+/* L15: */
+    }
+L10:
+    i__1 = *m + 1;
+    simp1_(&a[a_offset], mp, np, &i__1, l1, &nl1, &c__0, &kp, &bmax);
+    if (bmax <= 1e-6 && a[*m + 2 + a_dim1] < -1e-6) {
+	*icase = -1;
+	return 0;
+    } else if (bmax <= 1e-6 && a[*m + 2 + a_dim1] <= 1e-6) {
+	m12 = *m1 + *m2 + 1;
+	if (m12 <= *m) {
+	    i__1 = *m;
+	    for (ip = m12; ip <= i__1; ++ip) {
+		if (iposv[ip] == ip + *n) {
+		    simp1_(&a[a_offset], mp, np, &ip, l1, &nl1, &c__1, &kp, &
+			    bmax);
+		    if (bmax > (float)0.) {
+			goto L1;
+		    }
+		}
+/* L16: */
+	    }
+	}
+	ir = 0;
+	--m12;
+	if (*m1 + 1 > m12) {
+	    goto L30;
+	}
+	i__1 = m12;
+	for (i__ = *m1 + 1; i__ <= i__1; ++i__) {
+	    if (l3[i__ - *m1 - 1] == 1) {
+		i__2 = *n + 1;
+		for (k = 1; k <= i__2; ++k) {
+		    a[i__ + 1 + k * a_dim1] = -a[i__ + 1 + k * a_dim1];
+/* L17: */
+		}
+	    }
+/* L18: */
+	}
+	goto L30;
+    }
+    simp2_(&a[a_offset], m, n, mp, np, l2, &nl2, &ip, &kp, &q1);
+    if (ip == 0) {
+	*icase = -1;
+	return 0;
+    }
+L1:
+    i__1 = *m + 1;
+    simp3_(&a[a_offset], mp, np, &i__1, n, &ip, &kp);
+    if (iposv[ip] >= *n + *m1 + *m2 + 1) {
+	i__1 = nl1;
+	for (k = 1; k <= i__1; ++k) {
+	    if (l1[k - 1] == kp) {
+		goto L2;
+	    }
+/* L19: */
+	}
+L2:
+	--nl1;
+	i__1 = nl1;
+	for (is = k; is <= i__1; ++is) {
+	    l1[is - 1] = l1[is];
+/* L21: */
+	}
+    } else {
+	if (iposv[ip] < *n + *m1 + 1) {
+	    goto L20;
+	}
+	kh = iposv[ip] - *m1 - *n;
+	if (l3[kh - 1] == 0) {
+	    goto L20;
+	}
+	l3[kh - 1] = 0;
+    }
+    a[*m + 2 + (kp + 1) * a_dim1] += (float)1.;
+    i__1 = *m + 2;
+    for (i__ = 1; i__ <= i__1; ++i__) {
+	a[i__ + (kp + 1) * a_dim1] = -a[i__ + (kp + 1) * a_dim1];
+/* L22: */
+    }
+L20:
+    is = izrov[kp];
+    izrov[kp] = iposv[ip];
+    iposv[ip] = is;
+    if (ir != 0) {
+	goto L10;
+    }
+L30:
+    simp1_(&a[a_offset], mp, np, &c__0, l1, &nl1, &c__0, &kp, &bmax);
+    if (bmax <= (float)0.) {
+	*icase = 0;
+	return 0;
+    }
+    simp2_(&a[a_offset], m, n, mp, np, l2, &nl2, &ip, &kp, &q1);
+    if (ip == 0) {
+	*icase = 1;
+	return 0;
+    }
+    simp3_(&a[a_offset], mp, np, m, n, &ip, &kp);
+    goto L20;
+} /* simplx_ */
+
+#ifdef __cplusplus
+	}
+#endif
diff --git a/ext/f2c_recipes/splie2.c b/ext/f2c_recipes/splie2.c
new file mode 100644
index 000000000..525d377d2
--- /dev/null
+++ b/ext/f2c_recipes/splie2.c
@@ -0,0 +1,65 @@
+/* splie2.f -- translated by f2c (version 20031025).
+   You must link the resulting object file with libf2c:
+	on Microsoft Windows system, link with libf2c.lib;
+	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
+	or, if you install libf2c.a in a standard place, with -lf2c -lm
+	-- in that order, at the end of the command line, as in
+		cc *.o -lf2c -lm
+	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
+
+		http://www.netlib.org/f2c/libf2c.zip
+*/
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+#include "f2c.h"
+
+/* Table of constant values */
+
+static doublereal c_b4 = 1e30;
+
+/* Subroutine */ int splie2_(doublereal *x1a, doublereal *x2a, doublereal *ya,
+	 integer *m, integer *n, doublereal *y2a)
+{
+    /* System generated locals */
+    integer ya_dim1, ya_offset, y2a_dim1, y2a_offset, i__1, i__2;
+
+    /* Local variables */
+    static integer j, k;
+    static doublereal ytmp[100], y2tmp[100];
+    extern /* Subroutine */ int spline_(doublereal *, doublereal *, integer *,
+	     doublereal *, doublereal *, doublereal *);
+
+    /* Parameter adjustments */
+    --x1a;
+    y2a_dim1 = *m;
+    y2a_offset = 1 + y2a_dim1;
+    y2a -= y2a_offset;
+    ya_dim1 = *m;
+    ya_offset = 1 + ya_dim1;
+    ya -= ya_offset;
+    --x2a;
+
+    /* Function Body */
+    i__1 = *m;
+    for (j = 1; j <= i__1; ++j) {
+	i__2 = *n;
+	for (k = 1; k <= i__2; ++k) {
+	    ytmp[k - 1] = ya[j + k * ya_dim1];
+/* L11: */
+	}
+	spline_(&x2a[1], ytmp, n, &c_b4, &c_b4, y2tmp);
+	i__2 = *n;
+	for (k = 1; k <= i__2; ++k) {
+	    y2a[j + k * y2a_dim1] = y2tmp[k - 1];
+/* L12: */
+	}
+/* L13: */
+    }
+    return 0;
+} /* splie2_ */
+
+#ifdef __cplusplus
+	}
+#endif
diff --git a/ext/f2c_recipes/splin2.c b/ext/f2c_recipes/splin2.c
new file mode 100644
index 000000000..b69ff0dde
--- /dev/null
+++ b/ext/f2c_recipes/splin2.c
@@ -0,0 +1,66 @@
+/* splin2.f -- translated by f2c (version 20031025).
+   You must link the resulting object file with libf2c:
+	on Microsoft Windows system, link with libf2c.lib;
+	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
+	or, if you install libf2c.a in a standard place, with -lf2c -lm
+	-- in that order, at the end of the command line, as in
+		cc *.o -lf2c -lm
+	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
+
+		http://www.netlib.org/f2c/libf2c.zip
+*/
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+#include "f2c.h"
+
+/* Table of constant values */
+
+static doublereal c_b4 = 1e30;
+
+/* Subroutine */ int splin2_(doublereal *x1a, doublereal *x2a, doublereal *ya,
+	 doublereal *y2a, integer *m, integer *n, doublereal *x1, doublereal *
+	x2, doublereal *y)
+{
+    /* System generated locals */
+    integer ya_dim1, ya_offset, y2a_dim1, y2a_offset, i__1, i__2;
+
+    /* Local variables */
+    static integer j, k;
+    static doublereal ytmp[100], y2tmp[100], yytmp[100];
+    extern /* Subroutine */ int spline_(doublereal *, doublereal *, integer *,
+	     doublereal *, doublereal *, doublereal *), splint_(doublereal *, 
+	    doublereal *, doublereal *, integer *, doublereal *, doublereal *)
+	    ;
+
+    /* Parameter adjustments */
+    --x1a;
+    y2a_dim1 = *m;
+    y2a_offset = 1 + y2a_dim1;
+    y2a -= y2a_offset;
+    ya_dim1 = *m;
+    ya_offset = 1 + ya_dim1;
+    ya -= ya_offset;
+    --x2a;
+
+    /* Function Body */
+    i__1 = *m;
+    for (j = 1; j <= i__1; ++j) {
+	i__2 = *n;
+	for (k = 1; k <= i__2; ++k) {
+	    ytmp[k - 1] = ya[j + k * ya_dim1];
+	    y2tmp[k - 1] = y2a[j + k * y2a_dim1];
+/* L11: */
+	}
+	splint_(&x2a[1], ytmp, y2tmp, n, x2, &yytmp[j - 1]);
+/* L12: */
+    }
+    spline_(&x1a[1], yytmp, m, &c_b4, &c_b4, y2tmp);
+    splint_(&x1a[1], yytmp, y2tmp, m, x1, y);
+    return 0;
+} /* splin2_ */
+
+#ifdef __cplusplus
+	}
+#endif
diff --git a/ext/f2c_recipes/spline.c b/ext/f2c_recipes/spline.c
new file mode 100644
index 000000000..a05f81379
--- /dev/null
+++ b/ext/f2c_recipes/spline.c
@@ -0,0 +1,69 @@
+/* spline.f -- translated by f2c (version 20031025).
+   You must link the resulting object file with libf2c:
+	on Microsoft Windows system, link with libf2c.lib;
+	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
+	or, if you install libf2c.a in a standard place, with -lf2c -lm
+	-- in that order, at the end of the command line, as in
+		cc *.o -lf2c -lm
+	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
+
+		http://www.netlib.org/f2c/libf2c.zip
+*/
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+#include "f2c.h"
+
+/* Subroutine */ int spline_(doublereal *x, doublereal *y, integer *n, 
+	doublereal *yp1, doublereal *ypn, doublereal *y2)
+{
+    /* System generated locals */
+    integer i__1;
+
+    /* Local variables */
+    static integer i__, k;
+    static doublereal p, u[100], qn, un, sig;
+
+    /* Parameter adjustments */
+    --y2;
+    --y;
+    --x;
+
+    /* Function Body */
+    if (*yp1 > 9.9e29) {
+	y2[1] = (float)0.;
+	u[0] = (float)0.;
+    } else {
+	y2[1] = -.5;
+	u[0] = 3. / (x[2] - x[1]) * ((y[2] - y[1]) / (x[2] - x[1]) - *yp1);
+    }
+    i__1 = *n - 1;
+    for (i__ = 2; i__ <= i__1; ++i__) {
+	sig = (x[i__] - x[i__ - 1]) / (x[i__ + 1] - x[i__ - 1]);
+	p = sig * y2[i__ - 1] + 2.;
+	y2[i__] = (sig - (float)1.) / p;
+	u[i__ - 1] = (((y[i__ + 1] - y[i__]) / (x[i__ + 1] - x[i__]) - (y[i__]
+		 - y[i__ - 1]) / (x[i__] - x[i__ - 1])) * 6. / (x[i__ + 1] - 
+		x[i__ - 1]) - sig * u[i__ - 2]) / p;
+/* L11: */
+    }
+    if (*ypn > 9.9e29) {
+	qn = 0.;
+	un = 0.;
+    } else {
+	qn = .5;
+	un = 3. / (x[*n] - x[*n - 1]) * (*ypn - (y[*n] - y[*n - 1]) / (x[*n] 
+		- x[*n - 1]));
+    }
+    y2[*n] = (un - qn * u[*n - 2]) / (qn * y2[*n - 1] + 1.);
+    for (k = *n - 1; k >= 1; --k) {
+	y2[k] = y2[k] * y2[k + 1] + u[k - 1];
+/* L12: */
+    }
+    return 0;
+} /* spline_ */
+
+#ifdef __cplusplus
+	}
+#endif
diff --git a/ext/f2c_recipes/splint.c b/ext/f2c_recipes/splint.c
new file mode 100644
index 000000000..c408fed20
--- /dev/null
+++ b/ext/f2c_recipes/splint.c
@@ -0,0 +1,68 @@
+/* splint.f -- translated by f2c (version 20031025).
+   You must link the resulting object file with libf2c:
+	on Microsoft Windows system, link with libf2c.lib;
+	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
+	or, if you install libf2c.a in a standard place, with -lf2c -lm
+	-- in that order, at the end of the command line, as in
+		cc *.o -lf2c -lm
+	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
+
+		http://www.netlib.org/f2c/libf2c.zip
+*/
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+#include "f2c.h"
+
+/* Subroutine */ int splint_(doublereal *xa, doublereal *ya, doublereal *y2a, 
+	integer *n, doublereal *x, doublereal *y)
+{
+    /* System generated locals */
+    doublereal d__1, d__2, d__3;
+
+    /* Builtin functions */
+    /* Subroutine */ int s_paus(char *, ftnlen);
+
+    /* Local variables */
+    static doublereal a, b, h__;
+    static integer k, khi, klo;
+
+    /* Parameter adjustments */
+    --y2a;
+    --ya;
+    --xa;
+
+    /* Function Body */
+    klo = 1;
+    khi = *n;
+L1:
+    if (khi - klo > 1) {
+	k = (khi + klo) / 2;
+	if (xa[k] > *x) {
+	    khi = k;
+	} else {
+	    klo = k;
+	}
+	goto L1;
+    }
+    h__ = xa[khi] - xa[klo];
+    if (h__ == (float)0.) {
+	s_paus("Bad XA input.", (ftnlen)13);
+    }
+    a = (xa[khi] - *x) / h__;
+    b = (*x - xa[klo]) / h__;
+/* Computing 3rd power */
+    d__1 = a;
+/* Computing 3rd power */
+    d__2 = b;
+/* Computing 2nd power */
+    d__3 = h__;
+    *y = a * ya[klo] + b * ya[khi] + ((d__1 * (d__1 * d__1) - a) * y2a[klo] + 
+	    (d__2 * (d__2 * d__2) - b) * y2a[khi]) * (d__3 * d__3) / 6.;
+    return 0;
+} /* splint_ */
+
+#ifdef __cplusplus
+	}
+#endif