6 STATIC double *
d3_np_fs(
long n,
double a[],
const double b[],
double x[] );
78 for( i = 0; i < n; i++ )
88 for( i = 1; i < n; i++ )
90 xmult = a[2+(i-1)*3] / a[1+(i-1)*3];
91 a[1+i*3] = a[1+i*3] - xmult * a[0+i*3];
92 x[i] = b[i] - xmult * x[i-1];
95 x[n-1] = x[n-1] / a[1+(n-1)*3];
96 for( i = n-2; 0 <= i; i-- )
98 x[i] = ( x[i] - a[0+(i+1)*3] * x[i+1] ) / a[1+i*3];
106 int ibcbeg,
double ybcbeg,
int ibcend,
double ybcend )
228 for( i = 0; i < n - 1; i++ )
232 fprintf(
ioQQQ,
"SPLINE_CUBIC_SET - Fatal error!\n" );
233 fprintf(
ioQQQ,
" The knots must be strictly increasing\n" );
239 double *a = (
double*)
MALLOC((
size_t)(3*n)*
sizeof(
double));
240 double *b = (
double*)
MALLOC((
size_t)n*
sizeof(double));
250 else if( ibcbeg == 1 )
252 b[0] = ( y[1] - y[0] ) / ( t[1] - t[0] ) - ybcbeg;
253 a[1+0*3] = ( t[1] - t[0] ) / 3.0;
254 a[0+1*3] = ( t[1] - t[0] ) / 6.0;
256 else if( ibcbeg == 2 )
264 fprintf(
ioQQQ,
"SPLINE_CUBIC_SET - Fatal error!\n" );
265 fprintf(
ioQQQ,
" IBCBEG must be 0, 1 or 2, but I found %d.\n", ibcbeg );
271 for( i = 1; i < n-1; i++ )
273 b[i] = ( y[i+1] - y[i] ) / ( t[i+1] - t[i] )
274 - ( y[i] - y[i-1] ) / ( t[i] - t[i-1] );
275 a[2+(i-1)*3] = ( t[i] - t[i-1] ) / 6.0;
276 a[1+ i *3] = ( t[i+1] - t[i-1] ) / 3.0;
277 a[0+(i+1)*3] = ( t[i+1] - t[i] ) / 6.0;
288 else if( ibcend == 1 )
290 b[n-1] = ybcend - ( y[n-1] - y[n-2] ) / ( t[n-1] - t[n-2] );
291 a[2+(n-2)*3] = ( t[n-1] - t[n-2] ) / 6.0;
292 a[1+(n-1)*3] = ( t[n-1] - t[n-2] ) / 3.0;
294 else if( ibcend == 2 )
302 fprintf(
ioQQQ,
"SPLINE_CUBIC_SET - Fatal error!\n" );
303 fprintf(
ioQQQ,
" IBCEND must be 0, 1 or 2, but I found %d.\n", ibcend );
309 if( n == 2 && ibcbeg == 0 && ibcend == 0 )
316 if(
d3_np_fs( n, a, b, ypp ) == NULL )
318 fprintf(
ioQQQ,
"SPLINE_CUBIC_SET - Fatal error!\n" );
319 fprintf(
ioQQQ,
" The linear system could not be solved.\n" );
331 void spline_cubic_val(
long n,
const double t[],
double tval,
const double y[],
const double ypp[],
332 double *yval,
double *ypval,
double *yppval )
401 double dt = tval - t[ival];
402 double h = t[ival+1] - t[ival];
407 + dt * ( ( y[ival+1] - y[ival] ) / h
408 - ( ypp[ival+1] / 6.0 + ypp[ival] / 3.0 ) * h
409 + dt * ( 0.5 * ypp[ival]
410 + dt * ( ( ypp[ival+1] - ypp[ival] ) / ( 6.0 * h ) ) ) );
414 *ypval = ( y[ival+1] - y[ival] ) / h
415 - ( ypp[ival+1] / 6.0 + ypp[ival] / 3.0 ) * h
417 + dt * ( 0.5 * ( ypp[ival+1] - ypp[ival] ) / h ) );
421 *yppval = ypp[ival] + dt * ( ypp[ival+1] - ypp[ival] ) / h;
442 for(
long i=0; i < n; i++ )
445 for(
long j=0; j < n; j++ )
448 l *= (xval-x[j])/(x[i]-x[j]);
469 else if( xval >= x[n-1] )
474 double deriv = (y[ilo+1]-y[ilo])/(x[ilo+1]-x[ilo]);
475 yval = y[ilo] + deriv*(xval-x[ilo]);