dtable.c in tioga-1.11

- old
+ new

@@ -1743,10 +1743,76 @@
       rb_raise(rb_eRuntimeError, "corrupted data given to Dtable._load");
     }
   return ret;
 }
 
+/* The following function has been written by Benjamin ter Kuile <bterkuile@gmail.com> */
+
+PRIVATE
+/*
+ *  call-seq:
+ *     dtable.interpolate(Xs, Ys, nx, ny, x_start, x_end, y_start, y_end)  -> a_dtable
+ *  
+ *  Returns a copy of _dtable_ with the values interpolated given the proper X and Y axis to create a uniform spaced result in the X- and Y 
+ *  direction consisting of nx- and ny values for each direction.
+ */ VALUE dtable_interpolate(VALUE ary, VALUE x_vec, VALUE y_vec, VALUE nx_val, VALUE ny_val, VALUE xstart_val, VALUE xend_val, VALUE ystart_val, VALUE yend_val)
+{
+   Dtable *d = Get_Dtable(ary);
+	int nx = NUM2DBL(rb_Integer(nx_val));
+	int ny = NUM2DBL(rb_Integer(ny_val));
+   int i, j, num_cols = d->num_cols, num_rows = d->num_rows, last_row = num_rows - 1;
+	
+   long xsrc_len, ysrc_len;
+   double *xsrc = Dvector_Data_for_Read(x_vec, &xsrc_len);
+   double *ysrc = Dvector_Data_for_Read(y_vec, &ysrc_len);
+	if(xsrc_len != num_cols) rb_raise(rb_eArgError, "Number of x values (%d) do not match the number of columns (%d)", xsrc_len, num_cols);
+	if(ysrc_len != num_rows) rb_raise(rb_eArgError, "Number of y values (%d) do not match the number of rows (%d)", ysrc_len, num_rows);
+   VALUE new = dtable_init(dtable_alloc(cDtable), nx, ny);
+   Dtable *d2 = Get_Dtable(new);
+   double **src, **dest;
+	xstart_val = rb_Float(xstart_val);
+	double xstart = NUM2DBL(rb_Float(xstart_val));
+	if(xstart < xsrc[0]) rb_raise(rb_eArgError, "The start x value %g is smaller than the bound (%g)", xstart, xsrc[0]);
+	double xend = NUM2DBL(rb_Float(xend_val));
+	if(xend > xsrc[xsrc_len-1]) rb_raise(rb_eArgError, "The end x value %g is bigger than the bound (%g)", xend, xsrc[xsrc_len-1]);
+	double ystart = NUM2DBL(rb_Float(ystart_val));
+	if(ystart < ysrc[0]) rb_raise(rb_eArgError, "The start y value %g is smaller than the bound (%g)", ystart, ysrc[0]);
+	double yend = NUM2DBL(rb_Float(yend_val));
+	if(yend > ysrc[ysrc_len-1]) rb_raise(rb_eArgError, "The end y value %g is bigger than the bound (%g)", yend, ysrc[ysrc_len-1]);
+	double dx = (xend-xstart)/(nx-1);
+	double dy = (yend-ystart)/(ny-1);
+	double xcurrent = xstart;
+	double ycurrent = ystart;
+	double intvalue;
+	int isrc = 1;
+	int jsrc = 1;
+   src = d->ptr; dest = d2->ptr;
+   for (i = 1; i < ny+1; i++) {
+		while(ysrc[isrc] < ycurrent && ycurrent < ysrc[ysrc_len-1]){
+			isrc++;
+		}
+      for (j = 1; j < nx+1; j++) {
+			while(xsrc[jsrc] < xcurrent && xcurrent < xsrc[xsrc_len-1]){
+				jsrc++;
+			}
+			intvalue = (
+				( src[isrc-1][jsrc-1]*(ysrc[isrc]-ycurrent)*(xsrc[jsrc]-xcurrent)) +
+				( src[isrc][jsrc - 1] * (ycurrent - ysrc[isrc - 1]) * (xsrc[jsrc] - xcurrent) ) +
+				( src[isrc - 1][jsrc] * (ysrc[isrc] - ycurrent) * (xcurrent - xsrc[jsrc - 1]) ) +
+				( src[isrc][jsrc] * (ycurrent - ysrc[isrc - 1]) * (xcurrent - xsrc[jsrc - 1]) )
+			);
+			intvalue = intvalue / ( (ysrc[isrc] - ysrc[isrc - 1]) * (xsrc[jsrc] - xsrc[jsrc - 1]) );
+         dest[i-1][j-1] = intvalue;
+			xcurrent += dx;
+      }
+		xcurrent = xstart;
+		jsrc = 1;
+		ycurrent += dy;
+   }
+   return new;
+}
+
 /* 
  * Document-class: Dobjects::Dtable
  *
  * Dtables are a specialized implementation of two-dimensional arrays of double precision floating point numbers. 
  * They are intended for use in applications needing efficient processing of large 2D tables of numeric data.
@@ -1898,9 +1964,11 @@
    rb_define_method(cDtable, "safe_log10!", dtable_safe_log10_bang, -1);
    rb_define_method(cDtable, "safe_inv!", dtable_safe_inv_bang, -1);
    rb_define_method(cDtable, "safe_sqrt!", dtable_safe_sqrt_bang, 0);
    rb_define_method(cDtable, "safe_asin!", dtable_safe_asin_bang, 0);
    rb_define_method(cDtable, "safe_acos!", dtable_safe_acos_bang, 0);
+
+   rb_define_method(cDtable, "interpolate", dtable_interpolate, 8);
 
    /* Marshal : */
    rb_define_method(cDtable, "_dump", dtable_dump, 1);
    rb_define_singleton_method(cDtable, "_load", dtable_load, 1);
    /* modified by Vincent Fourmond, for splitting out the libraries */