Hi Brian
> What do you think about returning NAN (instead of 0 and y0,yn) for
> values outside the range? It's a more obvious failure.
You're right but I oriented on gsl_interp_eval_..._e functions and did
like that. Returning NaN is a good idea, though. It makes possible to
write a workaround to have any number as output instead of NaN. I
attached the patch to original interp.c.
Regards,
Evgeny Kurbatov
diff --git a/interpolation/interp.c b/interpolation/interp.c
index bb75e10..352861e 100644
--- a/interpolation/interp.c
+++ b/interpolation/interp.c
@@ -122,14 +122,9 @@ gsl_interp_eval_e (const gsl_interp * interp,
const double xa[], const double ya[], double x,
gsl_interp_accel * a, double *y)
{
- if (x < interp->xmin)
+ if (x < interp->xmin || x > interp->xmax)
{
- *y = ya[0];
- return GSL_EDOM;
- }
- else if (x > interp->xmax)
- {
- *y = ya[interp->size - 1];
+ *y = GSL_NAN;
return GSL_EDOM;
}
@@ -142,7 +137,14 @@ gsl_interp_eval (const gsl_interp * interp,
gsl_interp_accel * a)
{
double y;
- int status = interp->type->eval (interp->state, xa, ya, interp->size, x, a, &y);
+ int status;
+
+ if (x < interp->xmin || x > interp->xmax)
+ {
+ GSL_ERROR_VAL("interpolation error", GSL_EDOM, GSL_NAN);
+ }
+
+ status = interp->type->eval (interp->state, xa, ya, interp->size, x, a, &y);
DISCARD_STATUS(status);
@@ -156,14 +158,9 @@ gsl_interp_eval_deriv_e (const gsl_interp * interp,
gsl_interp_accel * a,
double *dydx)
{
- if (x < interp->xmin)
+ if (x < interp->xmin || x > interp->xmax)
{
- *dydx = 0.0;
- return GSL_EDOM;
- }
- else if (x > interp->xmax)
- {
- *dydx = 0.0;
+ *dydx = GSL_NAN;
return GSL_EDOM;
}
@@ -176,7 +173,14 @@ gsl_interp_eval_deriv (const gsl_interp * interp,
gsl_interp_accel * a)
{
double dydx;
- int status = interp->type->eval_deriv (interp->state, xa, ya, interp->size, x, a, &dydx);
+ int status;
+
+ if (x < interp->xmin || x > interp->xmax)
+ {
+ GSL_ERROR_VAL("interpolation error", GSL_EDOM, GSL_NAN);
+ }
+
+ status = interp->type->eval_deriv (interp->state, xa, ya, interp->size, x, a, &dydx);
DISCARD_STATUS(status);
@@ -190,14 +194,9 @@ gsl_interp_eval_deriv2_e (const gsl_interp * interp,
gsl_interp_accel * a,
double * d2)
{
- if (x < interp->xmin)
- {
- *d2 = 0.0;
- return GSL_EDOM;
- }
- else if (x > interp->xmax)
+ if (x < interp->xmin || x > interp->xmax)
{
- *d2 = 0.0;
+ *d2 = GSL_NAN;
return GSL_EDOM;
}
@@ -210,7 +209,14 @@ gsl_interp_eval_deriv2 (const gsl_interp * interp,
gsl_interp_accel * a)
{
double d2;
- int status = interp->type->eval_deriv2 (interp->state, xa, ya, interp->size, x, a, &d2);
+ int status;
+
+ if (x < interp->xmin || x > interp->xmax)
+ {
+ GSL_ERROR_VAL("interpolation error", GSL_EDOM, GSL_NAN);
+ }
+
+ status = interp->type->eval_deriv2 (interp->state, xa, ya, interp->size, x, a, &d2);
DISCARD_STATUS(status);
@@ -227,7 +233,7 @@ gsl_interp_eval_integ_e (const gsl_interp * interp,
{
if (a > b || a < interp->xmin || b > interp->xmax)
{
- *result = 0.0;
+ *result = GSL_NAN;
return GSL_EDOM;
}
else if(a == b)
@@ -247,7 +253,18 @@ gsl_interp_eval_integ (const gsl_interp * interp,
gsl_interp_accel * acc)
{
double result;
- int status = interp->type->eval_integ (interp->state, xa, ya, interp->size, acc, a, b, &result);
+ int status;
+
+ if (a > b || a < interp->xmin || b > interp->xmax)
+ {
+ GSL_ERROR_VAL("interpolation error", GSL_EDOM, GSL_NAN);
+ }
+ else if(a == b)
+ {
+ return 0.0;
+ }
+
+ status = interp->type->eval_integ (interp->state, xa, ya, interp->size, acc, a, b, &result);
DISCARD_STATUS(status);
_______________________________________________
Bug-gsl mailing list
[email protected]
http://lists.gnu.org/mailman/listinfo/bug-gsl
