Index: branches/eam_branches/psLib.20230123/src/math/psSpline.c
===================================================================
--- branches/eam_branches/psLib.20230123/src/math/psSpline.c	(revision 42325)
+++ branches/eam_branches/psLib.20230123/src/math/psSpline.c	(revision 42326)
@@ -72,5 +72,5 @@
  *****************************************************************************/
 
-static psF32 *calculateSecondDerivs(
+void calculateSecondDerivs(
     const psSpline1D *mySpline,
     psF32 dyLower, // if not NAN, lower-bound 1st derivative is defined
@@ -82,7 +82,8 @@
     psS32 n = mySpline->n; // n is the number of knots
     psF32 *u   = (psF32 *) psAlloc(n * sizeof(psF32));
-    psF32 *d2y = (psF32 *) psAlloc(n * sizeof(psF32));
-    psF32 *X = mySpline->xKnots;
-    psF32 *Y = mySpline->yKnots;
+
+    psF32 *X   = mySpline->xKnots;
+    psF32 *Y   = mySpline->yKnots;
+    psF32 *d2y = mySpline->d2yKnots;
 
     if (isfinite(dyLower)) {
@@ -120,6 +121,5 @@
     psFree(u);
     psTrace("psLib.math", 4, "---- %s() end ----\n", __func__);
-
-    return(d2y);
+    return;
 }
 
@@ -151,4 +151,24 @@
 
     return(tmpSpline);
+}
+
+/** Create an empty 1D spline **/
+psSpline1D *psSpline1DCreate(
+    int nKnots)			///< number of knots
+{
+    psSpline1D *spline = psSpline1DAlloc();
+    spline->n = nKnots; // number of knots
+
+    spline->xKnots   = (psF32 *) psAlloc( spline->n * sizeof(psF32));
+    spline->yKnots   = (psF32 *) psAlloc( spline->n * sizeof(psF32));
+    spline->d2yKnots = (psF32 *) psAlloc( spline->n * sizeof(psF32));
+
+    // x & y can both be F32 or F64. should knots be F64?
+    for (psS32 i = 0 ; i < spline->n ; i++) {
+	spline->xKnots[i]   = NAN;
+	spline->yKnots[i]   = NAN;
+	spline->d2yKnots[i] = NAN;
+    }
+    return(spline);
 }
 
@@ -187,10 +207,5 @@
     PS_ASSERT_LONG_LARGER_THAN_OR_EQUAL(y->n, (long)2, NULL);
 
-    psSpline1D *spline = psSpline1DAlloc();
-    spline->n = y->n; // number of knots
-
-    spline->xKnots   = (psF32 *) psAlloc( spline->n * sizeof(psF32));
-    spline->yKnots   = (psF32 *) psAlloc( spline->n * sizeof(psF32));
-    spline->d2yKnots = (psF32 *) psAlloc( spline->n * sizeof(psF32));
+    psSpline1D *spline = psSpline1DCreate(y->n);
 
     // x & y can both be F32 or F64. should knots be F64?
@@ -201,5 +216,5 @@
 
     // Generate the second derivatives at each data point.
-    spline->d2yKnots = calculateSecondDerivs(spline, dyLower, dyUpper);
+    calculateSecondDerivs(spline, dyLower, dyUpper);
 
     psTrace("psLib.math", 3, "---- %s() end ----\n", __func__);
Index: branches/eam_branches/psLib.20230123/src/math/psSpline.h
===================================================================
--- branches/eam_branches/psLib.20230123/src/math/psSpline.h	(revision 42325)
+++ branches/eam_branches/psLib.20230123/src/math/psSpline.h	(revision 42326)
@@ -51,4 +51,8 @@
  */
 psSpline1D *psSpline1DAlloc(void) PS_ATTR_MALLOC;
+
+/** Create an empty 1D spline **/
+psSpline1D *psSpline1DCreate(
+    int nKnots);			///< number of knots
 
 /** Evaluates 1-D spline polynomials at a specific coordinate.