Changeset 34522

Timestamp:

Oct 7, 2012, 2:36:26 PM (14 years ago)

Author:

eugene

Message:

rename spline functions to specify type, add dbl version

File:

• branches/eam_branches/ipp-20120905/Ohana/src/opihi/lib.data/spline.c (modified) (3 diffs)

branches/eam_branches/ipp-20120905/Ohana/src/opihi/lib.data/spline.c

@@ -2 +2 @@
 
 /* construct the natural spline for x, y in y2 */
-void spline_construct (float *x, float *y, int N, float *y2) {
+void spline_construct_flt (float *x, float *y, int N, float *y2) {
 
   int i;
@@ -27 +27 @@
 
 /* evaluate spline for x, y, y2 at X */
-float spline_apply (float *x, float *y, float *y2, int N, float X) {
+float spline_apply_flt (float *x, float *y, float *y2, int N, float X) {
 
   int i, lo, hi;
@@ -58 +58 @@
 
 }
+
+/* construct the natural spline for x, y in y2 */
+void spline_construct_dbl (opihi_flt *x, opihi_flt *y, int N, opihi_flt *y2) {
+
+  int i;
+  opihi_flt dy, dx, *tmp;
+  
+  ALLOCATE (tmp, opihi_flt, N);
+
+  y2[0] = tmp[0] = 0.0;
+  
+  for (i = 1; i < N-1; i++) {
+    dx = (x[i+0] - x[i-1]) / (x[i+1] - x[i-1]);
+    dy = dx * y2[i-1] + 2.0;
+    y2[i] = (dx - 1.0) / dy;
+    tmp[i] = (y[i+1] - y[i+0]) / (x[i+1] - x[i+0]) - (y[i+0] - y[i-1]) / (x[i+0] - x[i-1]);
+    tmp[i] = (6.0 * tmp[i] / (x[i+1] - x[i-1]) - dx*tmp[i-1]) / dy;
+  }
+  
+  y2[N-1] = 0;
+  for (i = N-2; i >= 1; i--)
+    y2[i] = y2[i]*y2[i+1] + tmp[i];
+
+  free (tmp);
+}
+
+/* evaluate spline for x, y, y2 at X */
+opihi_flt spline_apply_dbl (opihi_flt *x, opihi_flt *y, opihi_flt *y2, int N, opihi_flt X) {
+
+  int i, lo, hi;
+  opihi_flt dx, a, b, value;
+  
+  /* find correct element in array (x must be sorted) */
+  lo = 0;
+  hi = N-1;
+  while (hi - lo > 1) {
+    i = 0.5*(hi+lo);
+    if (x[i] > X) {
+      hi = i;
+    } else {
+      lo = i;
+    }
+  }
+
+  /* error condition: duplicate abssisca */
+  dx = x[hi] - x[lo];
+  if (dx == 0.0) {
+    return (HUGE_VAL);
+  }
+
+  /* evaluate spline */
+  a = (x[hi] - X) / dx;
+  b = (X - x[lo]) / dx;
+
+  value = a*y[lo] + b*y[hi] + ((a*a*a - a)*y2[lo] + (b*b*b - b)*y2[hi])*(dx*dx) / 6.0;
+  return (value);
+
+}