/* * 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; }