/* Returns an interpolant that can be fed to Special_Paths#append_interpolant_to_path to make nice splines. Can be used this way: f = Function.new(x,y) t.append_interpolant_to_path(f.make_interpolant) t.stroke */ static VALUE function_make_interpolant(VALUE self) { VALUE x_vec = get_x_vector(self); VALUE y_vec = get_y_vector(self); VALUE cache; VALUE a_vec,b_vec,c_vec; VALUE ret_val; double *x, *y, *a, *b, *c, *y2; double delta_x; long size = function_sanity_check(self); long i; function_ensure_spline_data_present(self); cache = get_spline_vector(self); x = Dvector_Data_for_Read(x_vec,NULL); y = Dvector_Data_for_Read(y_vec,NULL); y2 = Dvector_Data_for_Read(cache,NULL); a_vec = rb_funcall(cDvector, idNew, 1, LONG2NUM(size)); a = Dvector_Data_for_Write(a_vec, NULL); b_vec = rb_funcall(cDvector, idNew, 1, LONG2NUM(size)); b = Dvector_Data_for_Write(b_vec, NULL); c_vec = rb_funcall(cDvector, idNew, 1, LONG2NUM(size)); c = Dvector_Data_for_Write(c_vec, NULL); /* from my computations, the formula is the following: A = (y_2n+1 - y_2n)/(6 * delta_x) B = 0.5 * y_2n C = (y_n+1 - y_n)/delta_x - (2 * y_2n + y_2n+1) * delta_x/6 */ for(i = 0; i < size - 1; i++) { delta_x = x[i+1] - x[i]; a[i] = (y2[i+1] - y2[i]) / (6.0 * delta_x); b[i] = 0.5 * y2[i]; c[i] = (y[i+1] - y[i])/delta_x - (2 * y2[i] + y2[i+1]) * (delta_x / 6.0); } a[i] = b[i] = c[i] = 0.0; ret_val = rb_ary_new(); rb_ary_push(ret_val, x_vec); rb_ary_push(ret_val, y_vec); rb_ary_push(ret_val, a_vec); rb_ary_push(ret_val, b_vec); rb_ary_push(ret_val, c_vec); return ret_val; }