-
Notifications
You must be signed in to change notification settings - Fork 17
/
spline.c
97 lines (71 loc) · 2.35 KB
/
spline.c
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
/* ------- file: -------------------------- spline.c ----------------
Version: rh2.0
Author: Han Uitenbroek ([email protected])
Last modified: Tue Feb 16 14:57:56 1999 --
-------------------------- ----------RH-- */
/* --- Cubic spline interpolation routines. -- -------------- */
#include <stdio.h>
#include <stdlib.h>
#include "rh.h"
/* --- Function prototypes -- -------------- */
/* --- Global variables -- -------------- */
static bool_t ascend;
static int Ntable;
static double *xtable, xmin, xmax;
/* ------- begin -------------------------- splineCoef.c ------------ */
static double *M = NULL, *ytable;
void splineCoef(int N, double *x, double *y)
{
register int j;
static double *u = NULL;
double p, *q, hj, hj1, D, D1, mu;
ascend = (x[1] > x[0]) ? TRUE : FALSE;
xmin = (ascend) ? x[0] : x[N-1];
xmax = (ascend) ? x[N-1] : x[0];
q = M = (double *) realloc(M, N * sizeof(double));
u = (double *) realloc(u, N * sizeof(double));
hj = x[1] - x[0];
D = (y[1] - y[0]) / hj;
q[0] = u[0] = 0.0;
for (j = 1; j < N-1; j++) {
hj1 = x[j+1] - x[j];
mu = hj / (hj + hj1);
D1 = (y[j+1] - y[j]) / hj1;
p = mu*q[j-1] + 2;
q[j] = (mu - 1) / p;
u[j] = ((D1 - D) * 6/(hj + hj1) - mu*u[j-1]) / p;
hj = hj1; D = D1;
}
M[N - 1] = 0.0;
for (j = N-2; j >= 0; j--) {
M[j] = q[j]*M[j+1] + u[j];
}
Ntable = N;
xtable = x; ytable = y;
}
/* ------- end ---------------------------- splineCoef.c ------------ */
/* ---------------------------------------- splineEval.c ------------ */
void splineEval(int N, double *x, double *y, bool_t hunt)
{
register int n;
int j = 0;
double hj, fx, fx1;
for (n = 0; n < N; n++) {
if (x[n] <= xmin)
y[n] = (ascend) ? ytable[0] : ytable[Ntable-1];
else if (x[n] >= xmax)
y[n] = (ascend) ? ytable[Ntable-1] : ytable[0];
else {
if (hunt)
Hunt(Ntable, xtable, x[n], &j);
else
Locate(Ntable, xtable, x[n], &j);
hj = xtable[j+1] - xtable[j];
fx = (x[n] - xtable[j]) / hj;
fx1 = 1 - fx;
y[n] = fx1*ytable[j] + fx*ytable[j+1] +
(fx1*(SQ(fx1) - 1) * M[j] + fx*(SQ(fx) - 1) * M[j+1]) * SQ(hj)/6.0;
}
}
}
/* ------- end ---------------------------- splineEval.c ------------ */