47 #define UNUSED __attribute__ ((unused))
52 #define RETURN_IF_NULL(x) if (!x) { return ; }
56 #include <gsl/gsl_math.h>
57 #include <gsl/gsl_errno.h>
59 #include <gsl/gsl_interp.h>
83 GSL_ERROR_NULL(
"failed to allocate space for state", GSL_ENOMEM);
86 state->
a = (
double *) malloc (size *
sizeof (
double));
91 GSL_ERROR_NULL(
"failed to allocate space for a", GSL_ENOMEM);
94 state->
b = (
double *) malloc (size *
sizeof (
double));
99 GSL_ERROR_NULL(
"failed to allocate space for b", GSL_ENOMEM);
102 state->
c = (
double *) malloc (size *
sizeof (
double));
104 if (state->
c == NULL)
107 GSL_ERROR_NULL(
"failed to allocate space for c", GSL_ENOMEM);
110 state->
d = (
double *) malloc (size *
sizeof (
double));
112 if (state->
d == NULL)
115 GSL_ERROR_NULL(
"failed to allocate space for d", GSL_ENOMEM);
118 state->
y_prime = (
double *) malloc (size *
sizeof (
double));
122 GSL_ERROR_NULL(
"failed to allocate space for y_prime", GSL_ENOMEM);
130 const double y_array[],
size_t size)
134 double *
a = state->
a;
135 double *b = state->
b;
136 double *
c = state->
c;
137 double *
d = state->
d;
138 double *y_prime = state->
y_prime;
145 double h0 = (x_array[1] - x_array[0]);
146 double s0 = (y_array[1] - y_array[0]) / h0;
152 for (
i = 1;
i < (size - 1);
i++)
157 double hi = (x_array[
i+1] - x_array[
i]);
158 double him1 = (x_array[
i] - x_array[
i - 1]);
161 double si = (y_array[
i+1] - y_array[
i]) / hi;
162 double sim1 = (y_array[
i] - y_array[
i - 1]) / him1;
165 pi = (sim1*hi + si*him1) / (him1 + hi);
170 GSL_MIN(fabs(si), 0.5*fabs(
pi)));
180 y_prime[size-1] = (y_array[size - 1] - y_array[size - 2]) /
181 (x_array[size - 1] - x_array[size - 2]);
184 for (
i = 0;
i < (size - 1);
i++)
186 double hi = (x_array[
i+1] - x_array[
i]);
187 double si = (y_array[
i+1] - y_array[
i]) / hi;
190 a[
i] = (y_prime[
i] + y_prime[
i+1] - 2*si) / hi / hi;
191 b[
i] = (3*si - 2*y_prime[
i] - y_prime[
i+1]) / hi;
226 const double x_array[],
const double UNUSED y_array[],
size_t size,
227 double x, gsl_interp_accel *
a,
double *y)
235 index = gsl_interp_accel_find (
a, x_array, size,
x);
239 index = gsl_interp_bsearch (x_array,
x, 0, size - 1);
244 const double x_lo = x_array[index];
245 const double delx =
x - x_lo;
246 const double a_ = state->
a[index];
247 const double b = state->
b[index];
248 const double c = state->
c[index];
249 const double d = state->
d[index];
252 *
y =
d + delx*(
c + delx*(b + delx*a_));
260 const double x_array[],
const double UNUSED y_array[],
size_t size,
261 double x, gsl_interp_accel *
a,
double *dydx)
271 index = gsl_interp_accel_find (
a, x_array, size,
x);
275 index = gsl_interp_bsearch (x_array,
x, 0, size - 1);
280 double x_lo = x_array[index];
281 double delx =
x - x_lo;
282 double a_ = state->
a[index];
283 double b = state->
b[index];
284 double c = state->
c[index];
287 *dydx =
c + delx*(2*b + delx*3*a_);
294 const double x_array[],
const double UNUSED y_array[],
size_t size,
295 double x, gsl_interp_accel *
a,
double *y_pp)
305 index = gsl_interp_accel_find (
a, x_array, size,
x);
309 index = gsl_interp_bsearch (x_array,
x, 0, size - 1);
314 const double x_lo = x_array[index];
315 const double delx =
x - x_lo;
316 const double a_ = state->
a[index];
317 const double b = state->
b[index];
318 *y_pp = 6*a_*delx + 2*b;
325 const double x_array[],
const double UNUSED y_array[],
size_t size,
326 gsl_interp_accel * acc,
double a,
double b,
333 size_t i, index_a, index_b;
340 index_a = gsl_interp_accel_find (acc, x_array, size,
a);
341 index_b = gsl_interp_accel_find (acc, x_array, size, b);
345 index_a = gsl_interp_bsearch (x_array,
a, 0, size - 1);
346 index_b = gsl_interp_bsearch (x_array, b, 0, size - 1);
353 for(
i=index_a;
i<=index_b;
i++)
355 const double x_hi = x_array[
i + 1];
356 const double x_lo = x_array[
i];
357 const double dx = x_hi - x_lo;
364 double x1 = (
i == index_a) ?
a-x_lo : 0.0;
365 double x2 = (
i == index_b) ? b-x_lo : x_hi-x_lo;
367 *result += (1.0/4.0)*state->
a[
i]*(x2*x2*x2*x2 - x1*x1*x1*x1)
368 +(1.0/3.0)*state->
b[
i]*(x2*x2*x2 - x1*x1*x1)
369 +(1.0/2.0)*state->
c[
i]*(x2*x2 - x1*x1)
370 +state->
d[
i]*(x2-x1);
385 if ((x < 0 && y > 0) || (
x > 0 &&
y < 0))
static vector d(const double L_norm, const double J_norm, const vector roots)
Internal function that returns the coefficients "d_0", "d_2" and "d_4" from 1703.03967 corresponding ...
static int steffen_eval_integ(const void *vstate, const double x_array[], const double UNUSED y_array[], size_t size, gsl_interp_accel *acc, double a, double b, double *result)
static int steffen_eval_deriv2(const void *vstate, const double x_array[], const double UNUSED y_array[], size_t size, double x, gsl_interp_accel *a, double *y_pp)
static int steffen_eval(const void *vstate, const double x_array[], const double UNUSED y_array[], size_t size, double x, gsl_interp_accel *a, double *y)
static void * steffen_alloc(size_t size)
static int steffen_init(void *vstate, const double x_array[], const double y_array[], size_t size)
#define RETURN_IF_NULL(x)
const gsl_interp_type * lal_gsl_interp_steffen
static const gsl_interp_type steffen_type
static void steffen_free(void *vstate)
static double steffen_copysign(const double x, const double y)
static int steffen_eval_deriv(const void *vstate, const double x_array[], const double UNUSED y_array[], size_t size, double x, gsl_interp_accel *a, double *dydx)