version 1.1, 2000/01/09 17:00:52 |
version 1.1.1.3, 2003/09/15 07:09:25 |
Line 188 static int solve_five_diag __PROTO((five_diag m[], dou |
|
Line 188 static int solve_five_diag __PROTO((five_diag m[], dou |
|
static spline_coeff *cp_approx_spline __PROTO((struct curve_points * plot, int first_point, int num_points)); |
static spline_coeff *cp_approx_spline __PROTO((struct curve_points * plot, int first_point, int num_points)); |
static spline_coeff *cp_tridiag __PROTO((struct curve_points * plot, int first_point, int num_points)); |
static spline_coeff *cp_tridiag __PROTO((struct curve_points * plot, int first_point, int num_points)); |
static void do_cubic __PROTO((struct curve_points * plot, spline_coeff * sc, int first_point, int num_points, struct coordinate * dest)); |
static void do_cubic __PROTO((struct curve_points * plot, spline_coeff * sc, int first_point, int num_points, struct coordinate * dest)); |
static int compare_points __PROTO((struct coordinate * p1, struct coordinate * p2)); |
int compare_points __PROTO((SORTFUNC_ARGS p1, SORTFUNC_ARGS p2)); |
|
|
|
|
/* |
/* |
Line 251 struct curve_points *plot; |
|
Line 251 struct curve_points *plot; |
|
*/ |
*/ |
|
|
|
|
|
/* HBB 990205: rewrote the 'bezier' interpolation routine, |
|
* to prevent numerical overflow and other undesirable things happening |
|
* for large data files (num_data about 1000 or so), where binomial |
|
* coefficients would explode, and powers of 'sr' (0 < sr < 1) become |
|
* extremely small. Method used: compute logarithms of these |
|
* extremely large and small numbers, and only go back to the |
|
* real numbers once they've cancelled out each other, leaving |
|
* a reasonable-sized one. */ |
|
|
/* |
/* |
* cp_binomial() computes the binomial coefficients needed for BEZIER stuff |
* cp_binomial() computes the binomial coefficients needed for BEZIER stuff |
* and stores them into an array which is hooked to sdat. |
* and stores them into an array which is hooked to sdat. |
|
|
|
|
n = points - 1; |
n = points - 1; |
e = n / 2; |
e = n / 2; |
coeff[0] = 1.0; |
/* HBB 990205: calculate these in 'logarithmic space', |
|
* as they become _very_ large, with growing n (4^n) */ |
|
coeff[0] = 0.0; |
|
|
for (k = 0; k < e; k++) { |
for (k = 0; k < e; k++) { |
coeff[k + 1] = coeff[k] * ((double) (n - k)) / ((double) (k + 1)); |
coeff[k + 1] = coeff[k] + log(((double) (n-k)) / ((double) (k+1))); |
} |
} |
|
|
for (k = n; k >= e; k--) |
for (k = n; k >= e; k--) |
Line 322 unsigned int exponent; |
|
Line 333 unsigned int exponent; |
|
return (y); |
return (y); |
} |
} |
|
|
|
|
static void eval_bezier(cp, first_point, num_points, sr, px, py, c) |
static void eval_bezier(cp, first_point, num_points, sr, px, py, c) |
struct curve_points *cp; |
struct curve_points *cp; |
int first_point; /* where to start in plot->points (to find x-range) */ |
int first_point; /* where to start in plot->points (to find x-range) */ |
|
|
*px = this_points[n].x; |
*px = this_points[n].x; |
*py = this_points[n].y; |
*py = this_points[n].y; |
} else { |
} else { |
|
/* HBB 990205: do calculation in 'logarithmic space', |
|
* to avoid over/underflow errors, which would exactly cancel |
|
* out each other, anyway, in an exact calculation |
|
*/ |
unsigned int i; |
unsigned int i; |
double lx = 0.0, ly = 0.0; |
double lx = 0.0, ly = 0.0; |
double sr_to_the_i = 1.0; |
double log_dsr_to_the_n = n * log(1 -sr); |
double dsr_to_the_di = s_pow(1 - sr, n); |
double log_sr_over_dsr = log(sr) - log(1 - sr); |
double reciprocal_dsr = 1.0 / (1 - sr); |
|
|
|
for (i = 0; i <= n; i++) { |
for (i = 0; i <= n; i++) { |
double u = c[i] * sr_to_the_i * s_pow(1 - sr, n - i); |
double u = exp(c[i] + log_dsr_to_the_n + i * log_sr_over_dsr); |
|
|
lx += this_points[i].x * u; |
lx += this_points[i].x * u; |
ly += this_points[i].y * u; |
ly += this_points[i].y * u; |
sr_to_the_i *= sr; |
|
dsr_to_the_di *= reciprocal_dsr; |
|
} |
} |
|
|
*px = lx; |
*px = lx; |
Line 909 struct curve_points *plot; |
|
Line 923 struct curve_points *plot; |
|
* (MGR 1992) |
* (MGR 1992) |
*/ |
*/ |
|
|
static int compare_points(p1, p2) |
int compare_points(argp1, argp2) |
struct coordinate *p1; |
SORTFUNC_ARGS argp1; |
struct coordinate *p2; |
SORTFUNC_ARGS argp2; |
{ |
{ |
|
const struct coordinate *p1 = argp1; |
|
const struct coordinate *p2 = argp2; |
|
|
if (p1->x > p2->x) |
if (p1->x > p2->x) |
return (1); |
return (1); |
if (p1->x < p2->x) |
if (p1->x < p2->x) |
Line 921 struct coordinate *p2; |
|
Line 938 struct coordinate *p2; |
|
} |
} |
|
|
void sort_points(plot) |
void sort_points(plot) |
struct curve_points *plot; |
struct curve_points *plot; |
{ |
{ |
int first_point, num_points; |
int first_point, num_points; |
|
|
Line 929 struct curve_points *plot; |
|
Line 946 struct curve_points *plot; |
|
while ((num_points = next_curve(plot, &first_point)) > 0) { |
while ((num_points = next_curve(plot, &first_point)) > 0) { |
/* Sort this set of points, does qsort handle 1 point correctly? */ |
/* Sort this set of points, does qsort handle 1 point correctly? */ |
qsort((char *) (plot->points + first_point), num_points, |
qsort((char *) (plot->points + first_point), num_points, |
sizeof(struct coordinate), (sortfunc) compare_points); |
sizeof(struct coordinate), compare_points); |
first_point += num_points; |
first_point += num_points; |
} |
} |
return; |
return; |