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
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
|
// interpolate.cc -- classes for 1D interpolation
// $Id$
//
// interpolation::Lagrange_uniform_1d::Lagrange_uniform
// interpolation::Lagrange_uniform_1d::setup_data_x
// interpolation::Lagrange_uniform_1d::setup_data_y
// interpolation::Lagrange_uniform_1d::bdod
// interpolation::Lagrange_uniform_1d::interpolate
// interpolation::Lagrange_uniform_1d::Jacobian
//
// ***** template instantiations *****
//
#include <assert.h>
#include <jt/stdc.h>
#include <jt/util.hh>
#include <jt/linear_map.hh>
#include <jt/interpolate.hh>
using jtutil::error_exit;
//******************************************************************************
//******************************************************************************
//******************************************************************************
//
// This function constructs an interpolation::Lagrange_uniform_1d object.
//
namespace interpolation
{
template <typename fpt>
Lagrange_uniform_1d<fpt>::Lagrange_uniform_1d(int order_in)
: order(order_in), N_interp(order_in+1),
x_map_ptr(NULL), y_array(NULL)
{
}
}
//*****************************************************************************
//
// This function is used by client code to set up the x coordinates of
// the (x,y) data points to be interpolated.
//
namespace interpolation
{
template <typename fpt>
Lagrange_uniform_1d<fpt>::setup_data_x(const linear_map<fpt>& x_map_in)
{
x_map_ptr = &x_map_in;
if (x_map_ptr->N_points() < N_interp)
then error_exit(ERROR_EXIT,
"***** Lagrange_uniform_1d<fpt>::setup_data_x():\n"
" not enough data points for this order of interpolation!\n"
" order=%d needs minimum of N_points=%d data points,\n"
" but x_map only provides %d!\n"
" x_map i_[min,max]=[%d,%d]\n"
" x=[%g(%g)%g]\n"
,
order, N_points,
x_map_ptr->N_points(),
x_map_ptr->min_int(), x_map_ptr->max_int(),
double(x_map_ptr->min_fp()),
double(x_map_ptr->delta_fp()),
double(x_map_ptr->max_fp())); /*NOTREACHED*/
}
}
//*****************************************************************************
//
// This function is used by client code to set up the y coordinates of
// the (x,y) data points to be interpolated.
//
namespace interpolation
{
template <typename fpt>
void Lagrange_uniform_1d<fpt>::setup_data_y(const fpt y_array_in[])
{
y_array = y_array_in;
}
}
//*****************************************************************************
//
// This function does the actual interpolation.
//
namespace interpolation
{
template <typename fpt>
fpt Lagrange_uniform_1d<fpt>::interpolate(fpt x, char end_action = 'o')
const
{
assert(x_map_ptr != NULL);
assert(y_array != NULL);
// FIXME FIXME FIXME
}
}
//*****************************************************************************
//
// This function computes the backwards domain of dependence of
// interpolate(x,end_action), i.e. it computes the interval of [i]
// such that y[i] is needed to compute interpolate(x,end_action) .
// Equivalently, this function defines the interpolator's centering
// policy.
//
namespace interpolation
{
template <typename fpt>
void Lagrange_uniform_1d<fpt>::bdod(fpt x, char end_action = 'o',
int& i_lo, int& i_hi)
const
{
//
// Suppose we have a data array indexed by [i], from which we wish to
// select an N-point subarray [i_lo, i_hi] for interpolation (and thus
// centered as close as possible) to the point x . The problem is
// just how to choose the centering.
//
// The following diagram illustrates the issues:
//
// N i_lo i_hi [i-3] [i-2] [i-1] [i] [i+1] [i+2] [i+3] [i+4]
// - ---- ---- ----- ----- ----- ----- ----- ----- ----- -----
// 2 i i+1 *xxxxxx*
// 3 i-1 i+1 * xxx*xxx *
// 4 i-1 i+2 * *xxxxxx* *
// 5 i-2 i+2 * * xxx*xxx * *
// 6 i-2 i+3 * * *xxxxxx* * *
// 7 i-3 i+3 * * * xxx*xxx * * *
// 8 i-3 i+4 * * * *xxxxxx* * * *
//
// The diagram shows the range of x values relative to the grid, for
// which we'll choose each centering. If N is even, the centering is
// based on which grid zone contains x (i.e. [i] is chosen as floor(x)),
// while if N is odd, the centering is based on which grid point is closest
// to x (i.e. [i] is chosen as round(x)).
//
assert(x_map_ptr != NULL);
const enum linear_map<fpt>::noninteger_action floor_or_round
= jtutil::is_even(N_interp) ? linear_map<fpt>::floor
: linear_map<fpt>::round);
const int i_center = x_map_ptr->int_of_fp(x, floor_or_round);
i_hi = i_center + (N_interp >> 1);
i_lo = i_hi - N_interp + 1;
//
// now adjust these endpoints as needed if we're off the end of the
// data array
//
switch (end_action)
{
case 'e':
if ( (i_lo < x_map_ptr->min_int()) || (i_hi > x_map_ptr->max_int()) )
then error_exit(ERROR_EXIT,
"***** Lagrange_uniform<fpt>::bdod():\n"
" x=%g is out of range!\n"
" x_map i_[min,max]=[%d,%d]\n"
" x=[%g(%g)%g]\n"
" i_center=%d i_[lo,hi]=[%d,%d]\n"
,
double(x),
x_map_ptr->min_int(), x_map_ptr->max_int(),
double(x_map_ptr->min_fp()),
double(x_map_ptr->delta_fp()),
double(x_map_ptr->max_fp()),
i_center, i_lo, i_hi); /*NOTREACHED*/
case 'o':
if (i_lo < x_map_ptr->min_int())
then {
const int shift = x_map_ptr->min_int() - i_lo;
i_center += shift;
i_lo += shift;
i_hi += shift;
break;
}
if (i_hi > x_map_ptr->max_int())
then {
const int shift = x_map_ptr->max_int() - i_hi;
i_center += shift;
i_lo += shift;
i_hi += shift;
break;
}
break;
case '*':
default:
error_exit(ERROR_EXIT,
"***** Lagrange_uniform<fpt>::bdod(): invalid end_action='%c'='0x%02x'!\n",
int(end_action), unsigned int(end_action)); /*NOTREACHED*/
}
// final sanity check (bug firewall)
if (! ( (i_lo >= x_map_ptr->min_int())
&& (i_hi <= x_map_ptr->max_int())
&& (jtutil::how_many_in_range(i_lo,i_hi) == N_points) ) )
then error_exit(PANIC_EXIT,
"***** Lagrange_uniform<fpt>::bdod():\n"
" inconsistent i_lo, i_hi, and/or N_points!\n"
" (this should never happen!)\n"
" x=%g\n"
" x_map i_[min,max]=[%d,%d]\n"
" x=[%g(%g)%g]\n"
" order=%d N_points=%d\n"
" i_center=%d i_[lo,hi]=[%d,%d]\n"
,
double(x),
x_map_ptr->min_int(), x_map_ptr->max_int(),
double(x_map_ptr->min_fp()),
double(x_map_ptr->delta_fp()),
double(x_map_ptr->max_fp()),
order, N_points,
i_center, i_lo, i_hi); /*NOTREACHED*/
return i_center;
}
}
//*****************************************************************************
//
// This function computes the Jacobian
// d interpolate(x,end_action)
// ---------------------------
// d y[i]
//
namespace interpolation
{
template <typename fpt>
fpt Lagrange_uniform_1d<fpt>::Jacobian(fpt x, char end_action = 'o', int i)
const
{
assert(x_map_ptr != NULL);
// FIXME FIXME FIXME
}
}
|
