1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
|
/* spl-array.c: Spline interpolation support routines for matrices
//
// Written and Copyright (C) 1993-1997 by Michael J. Gourlay
//
// PROVIDED AS IS. NO WARRANTEES, EXPRESS OR IMPLIED.
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "my_malloc.h"
#include "spline.h"
#include "spl-array.h"
/* --------------------------------------------------------------- */
#define NEAR(x1, x2) (((x2)!=0.0) && (((x1)/(x2)) >= 0.999) && (((x1)/(x2))<1.001))
#ifndef FALSE
#define FALSE 0
#endif
/* --------------------------------------------------------------- */
/* derivative_hack: compute 1st derivative of x,y data (len entries)
//
// Written and Copyright (C) 1994-1997 by Michael J. Gourlay
//
// PROVIDED AS IS. NO WARRANTEES, EXPRESS OR IMPLIED.
//
// Mathematically, it's a hack to prevent overshooting knots, but
// maintain smoothness and it works more intuitively. Besides, we all
// know that mathematicians are too worried about rigor, and the
// physicists end up creating the right math. -- MJG
*/
static void derivative_hack(const double *x, const double *y, double *yd, int len)
{
int indx;
if(x[0]==x[1]) {
yd[0] = 0.0; /* avoid division by zero */
} else {
yd[0] = (y[1]-y[0])/(x[1]-x[0]);
}
if(x[len-2] == x[len-1]) {
yd[len-1] = 0.0; /* avoid division by zero */
} else {
yd[len-1] = (y[len-1]-y[len-2])/(x[len-1]-x[len-2]);
}
for(indx=1; indx<(len-1); indx++) {
if(x[indx-1]==x[indx] || x[indx]==x[indx+1]) {
yd[indx] = 0.0; /* avoid division by zero */
} else {
if( (y[indx-1] - y[indx]) * (y[indx+1] - y[indx]) >= 0.0) {
/* There was a change in the sign of yd so force it zero */
/* This will prevent the spline from overshooting this knot */
yd[indx] = 0.0;
} else {
/* Set slope at this knot to slope between two adjacent knots */
yd[indx] = (y[indx-1]-y[indx+1]) / (x[indx-1]-x[indx+1]);
}
}
}
}
/* --------------------------------------------------------------- */
/* hermite3_array : cubic hermite interpolation for an array of points
//
// Uses derivative_hack to find derivatives.
//
//
// kx, ky (in): arrays of knots
// nk (in): number of knots
// sx (in): evaluation points array
// sy (out): spline values at evaluation pts sx. Must already be allocated.
// ns (in): number of evaluation points
*/
int
hermite3_array(const double *kx, const double *ky, long nk, double *sx,
double *sy, long ns)
{
register long xi;
double *kyd;
if((kyd=MY_CALLOC(nk, double))==NULL)
return 1;
/* Test bounds. */
/* As of 18jul94, this test was triggering for cases
// where the bounds were nearly equal, but slightly out of range, in
// which case the spline should work anyway, which is why I let it run
// even if the spline abcissas are out of range.
*/
if((sx[0] < kx[0]) || (sx[ns-1] > kx[nk-1])) {
if(!NEAR(sx[ns-1], kx[nk-1])) {
fprintf(stderr, "hermite3_array: out of range:\n");
fprintf(stderr,
"hermite3_array: eval=%.20g < knot=%.20g | %.20g>%.20g\n",
sx[0], kx[0], sx[ns-1], kx[nk-1]);
}
}
/* Find array of derivatives */
derivative_hack(kx, ky, kyd, nk);
/* Evaluate the spline */
for(xi=0; xi<ns; xi++) {
sy[xi]=hermite3_interp(sx[xi], kx, ky, kyd, nk, NULL, FALSE, NULL, NULL);
}
/* Free the derivatives array */
FREE(kyd);
return 0;
}
|
