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
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
|
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html><head><meta http-equiv="Content-Type" content="text/html;charset=iso-8859-1">
<title>Qwt User's Guide: qwt_spline.cpp Source File</title>
<link href="doxygen.css" rel="stylesheet" type="text/css">
</head><body>
<!-- Generated by Doxygen 1.3.8 -->
<div class="qindex"><a class="qindex" href="index.html">Main Page</a> | <a class="qindex" href="hierarchy.html">Class Hierarchy</a> | <a class="qindex" href="classes.html">Alphabetical List</a> | <a class="qindex" href="annotated.html">Class List</a> | <a class="qindex" href="files.html">File List</a> | <a class="qindex" href="functions.html">Class Members</a> | <a class="qindex" href="globals.html">File Members</a> | <a class="qindex" href="pages.html">Related Pages</a></div>
<h1>qwt_spline.cpp</h1><pre class="fragment"><div>00001 <span class="comment">/* -*- mode: C++ ; c-file-style: "stroustrup" -*- *****************************</span>
00002 <span class="comment"> * Qwt Widget Library</span>
00003 <span class="comment"> * Copyright (C) 1997 Josef Wilgen</span>
00004 <span class="comment"> * Copyright (C) 2002 Uwe Rathmann</span>
00005 <span class="comment"> * </span>
00006 <span class="comment"> * This library is free software; you can redistribute it and/or</span>
00007 <span class="comment"> * modify it under the terms of the Qwt License, Version 1.0</span>
00008 <span class="comment"> *****************************************************************************/</span>
00009
00010 <span class="preprocessor">#include "qwt_spline.h"</span>
00011 <span class="preprocessor">#include "<a class="code" href="qwt__math_8h.html">qwt_math.h</a>"</span>
00012 <span class="preprocessor">#include "qwt.h"</span>
00013
<a name="l00018"></a><a class="code" href="class_qwt_spline.html#a2">00018</a> <span class="keywordtype">double</span> <a class="code" href="class_qwt_spline.html#a2">QwtSpline::value</a>(<span class="keywordtype">double</span> x)<span class="keyword"> const</span>
00019 <span class="keyword"></span>{
00020 <span class="keywordflow">if</span> (!d_a)
00021 <span class="keywordflow">return</span> 0.0;
00022
00023 <span class="keyword">const</span> <span class="keywordtype">int</span> i = lookup(x);
00024
00025 <span class="keyword">const</span> <span class="keywordtype">double</span> delta = x - d_x[i];
00026 <span class="keywordflow">return</span>( ( ( ( d_a[i] * delta) + d_b[i] )
00027 * delta + d_c[i] ) * delta + d_y[i] );
00028 }
00029
<a name="l00031"></a><a class="code" href="class_qwt_spline.html#a0">00031</a> <a class="code" href="class_qwt_spline.html#a0">QwtSpline::QwtSpline</a>()
00032 {
00033 d_a = d_b = d_c = 0;
00034 d_xbuffer = d_ybuffer = d_x = d_y = 0;
00035 d_size = 0;
00036 d_buffered = 0;
00037 }
00038
<a name="l00061"></a><a class="code" href="class_qwt_spline.html#a5">00061</a> <span class="keywordtype">void</span> <a class="code" href="class_qwt_spline.html#a5">QwtSpline::copyValues</a>(<span class="keywordtype">int</span> tf)
00062 {
00063 cleanup();
00064 d_buffered = tf;
00065 }
00066
<a name="l00068"></a><a class="code" href="class_qwt_spline.html#a1">00068</a> <a class="code" href="class_qwt_spline.html#a1">QwtSpline::~QwtSpline</a>()
00069 {
00070 cleanup();
00071 }
00072
00074 <span class="keywordtype">int</span> QwtSpline::lookup(<span class="keywordtype">double</span> x)<span class="keyword"> const</span>
00075 <span class="keyword"></span>{
00076 <span class="keywordtype">int</span> i1, i2, i3;
00077
00078 <span class="keywordflow">if</span> (x <= d_x[0])
00079 i1 = 0;
00080 <span class="keywordflow">else</span> <span class="keywordflow">if</span> (x >= d_x[d_size - 2])
00081 i1 = d_size -2;
00082 <span class="keywordflow">else</span>
00083 {
00084 i1 = 0;
00085 i2 = d_size -2;
00086 i3 = 0;
00087
00088 <span class="keywordflow">while</span> ( i2 - i1 > 1 )
00089 {
00090 i3 = i1 + ((i2 - i1) >> 1);
00091
00092 <span class="keywordflow">if</span> (d_x[i3] > x)
00093 i2 = i3;
00094 <span class="keywordflow">else</span>
00095 i1 = i3;
00096
00097 }
00098 }
00099 <span class="keywordflow">return</span> i1;
00100 }
00101
00102
<a name="l00127"></a><a class="code" href="class_qwt_spline.html#a3">00127</a> <span class="keywordtype">int</span> <a class="code" href="class_qwt_spline.html#a3">QwtSpline::recalc</a>(<span class="keywordtype">double</span> *x, <span class="keywordtype">double</span> *y, <span class="keywordtype">int</span> n, <span class="keywordtype">int</span> periodic)
00128 {
00129 <span class="keywordtype">int</span> i, rv = 0;
00130
00131 cleanup();
00132
00133 <span class="keywordflow">if</span> (n > 2)
00134 {
00135 d_size = n;
00136
00137 <span class="keywordflow">if</span> (d_buffered)
00138 {
00139 d_xbuffer = <span class="keyword">new</span> <span class="keywordtype">double</span>[n-1];
00140 d_ybuffer = <span class="keyword">new</span> <span class="keywordtype">double</span>[n-1];
00141
00142 <span class="keywordflow">if</span> ((!d_xbuffer) || (!d_ybuffer))
00143 {
00144 cleanup();
00145 <span class="keywordflow">return</span> Qwt::ErrNoMem;
00146 }
00147 <span class="keywordflow">else</span>
00148 {
00149 <span class="keywordflow">for</span> (i=0;i<n;i++)
00150 {
00151 d_xbuffer[i] = x[i];
00152 d_ybuffer[i] = y[i];
00153 }
00154 d_x = d_xbuffer;
00155 d_y = d_ybuffer;
00156 }
00157 }
00158 <span class="keywordflow">else</span>
00159 {
00160 d_x = x;
00161 d_y = y;
00162 }
00163
00164 d_a = <span class="keyword">new</span> <span class="keywordtype">double</span>[n-1];
00165 d_b = <span class="keyword">new</span> <span class="keywordtype">double</span>[n-1];
00166 d_c = <span class="keyword">new</span> <span class="keywordtype">double</span>[n-1];
00167
00168 <span class="keywordflow">if</span> ( (!d_a) || (!d_b) || (!d_c) )
00169 {
00170 cleanup();
00171 <span class="keywordflow">return</span> Qwt::ErrMono;
00172 }
00173
00174 <span class="keywordflow">if</span>(periodic)
00175 rv = buildPerSpline();
00176 <span class="keywordflow">else</span>
00177 rv = buildNatSpline();
00178
00179 <span class="keywordflow">if</span> (rv) cleanup();
00180 }
00181
00182 <span class="keywordflow">return</span> rv;
00183 }
00184
<a name="l00206"></a><a class="code" href="class_qwt_spline.html#a4">00206</a> <span class="keywordtype">int</span> <a class="code" href="class_qwt_spline.html#a3">QwtSpline::recalc</a>(<span class="keyword">const</span> QwtArray<double> &x, <span class="keyword">const</span> QwtArray<double> &y,
00207 <span class="keywordtype">int</span> periodic)
00208 {
00209 <span class="keywordtype">int</span> n = QMIN(x.size(), y.size());
00210 d_buffered = TRUE;
00211
00212 <span class="keywordflow">return</span> <a class="code" href="class_qwt_spline.html#a3">recalc</a>(x.data(), y.data(), n, periodic);
00213 }
00214
00224 <span class="keywordtype">int</span> QwtSpline::buildNatSpline()
00225 {
00226 <span class="keywordtype">int</span> i;
00227 <span class="keywordtype">double</span> dy1, dy2;
00228
00229 <span class="keywordtype">double</span> *d = <span class="keyword">new</span> <span class="keywordtype">double</span>[d_size-1];
00230 <span class="keywordtype">double</span> *h = <span class="keyword">new</span> <span class="keywordtype">double</span>[d_size-1];
00231 <span class="keywordtype">double</span> *s = <span class="keyword">new</span> <span class="keywordtype">double</span>[d_size];
00232
00233 <span class="keywordflow">if</span> ( (!d) || (!h) || (!s) )
00234 {
00235 cleanup();
00236 <span class="keywordflow">if</span> (h) <span class="keyword">delete</span>[] h;
00237 <span class="keywordflow">if</span> (s) <span class="keyword">delete</span>[] s;
00238 <span class="keywordflow">if</span> (d) <span class="keyword">delete</span>[] d;
00239 <span class="keywordflow">return</span> Qwt::ErrNoMem;
00240 }
00241
00242 <span class="comment">//</span>
00243 <span class="comment">// set up tridiagonal equation system; use coefficient</span>
00244 <span class="comment">// vectors as temporary buffers</span>
00245 <span class="keywordflow">for</span> (i=0; i<d_size - 1; i++)
00246 {
00247 h[i] = d_x[i+1] - d_x[i];
00248 <span class="keywordflow">if</span> (h[i] <= 0)
00249 {
00250 <span class="keyword">delete</span>[] h;
00251 <span class="keyword">delete</span>[] s;
00252 <span class="keyword">delete</span>[] d;
00253 <span class="keywordflow">return</span> Qwt::ErrMono;
00254 }
00255 }
00256
00257 dy1 = (d_y[1] - d_y[0]) / h[0];
00258 <span class="keywordflow">for</span> (i = 1; i < d_size - 1; i++)
00259 {
00260 d_b[i] = d_c[i] = h[i];
00261 d_a[i] = 2.0 * (h[i-1] + h[i]);
00262
00263 dy2 = (d_y[i+1] - d_y[i]) / h[i];
00264 d[i] = 6.0 * ( dy1 - dy2);
00265 dy1 = dy2;
00266 }
00267
00268 <span class="comment">//</span>
00269 <span class="comment">// solve it</span>
00270 <span class="comment">//</span>
00271
00272 <span class="comment">// L-U Factorization</span>
00273 <span class="keywordflow">for</span>(i = 1; i < d_size - 2;i++)
00274 {
00275 d_c[i] /= d_a[i];
00276 d_a[i+1] -= d_b[i] * d_c[i];
00277 }
00278
00279 <span class="comment">// forward elimination</span>
00280 s[1] = d[1];
00281 <span class="keywordflow">for</span>(i=2;i<d_size - 1;i++)
00282 s[i] = d[i] - d_c[i-1] * s[i-1];
00283
00284 <span class="comment">// backward elimination</span>
00285 s[d_size - 2] = - s[d_size - 2] / d_a[d_size - 2];
00286 <span class="keywordflow">for</span> (i= d_size -3; i > 0; i--)
00287 s[i] = - (s[i] + d_b[i] * s[i+1]) / d_a[i];
00288
00289 <span class="comment">//</span>
00290 <span class="comment">// Finally, determine the spline coefficients</span>
00291 <span class="comment">//</span>
00292 s[d_size - 1] = s[0] = 0.0;
00293 <span class="keywordflow">for</span> (i = 0; i < d_size - 1; i++)
00294 {
00295 d_a[i] = ( s[i+1] - s[i] ) / ( 6.0 * h[i]);
00296 d_b[i] = 0.5 * s[i];
00297 d_c[i] = ( d_y[i+1] - d_y[i] )
00298 / h[i] - (s[i+1] + 2.0 * s[i] ) * h[i] / 6.0;
00299 }
00300
00301 <span class="keyword">delete</span>[] d;
00302 <span class="keyword">delete</span>[] s;
00303 <span class="keyword">delete</span>[] h;
00304
00305 <span class="keywordflow">return</span> 0;
00306
00307 }
00308
00318 <span class="keywordtype">int</span> QwtSpline::buildPerSpline()
00319 {
00320 <span class="keywordtype">int</span> i,imax;
00321 <span class="keywordtype">double</span> sum;
00322 <span class="keywordtype">double</span> dy1, dy2,htmp;
00323
00324 <span class="keywordtype">double</span> *d = <span class="keyword">new</span> <span class="keywordtype">double</span>[d_size-1];
00325 <span class="keywordtype">double</span> *h = <span class="keyword">new</span> <span class="keywordtype">double</span>[d_size-1];
00326 <span class="keywordtype">double</span> *s = <span class="keyword">new</span> <span class="keywordtype">double</span>[d_size];
00327
00328 <span class="keywordflow">if</span> ( (!d) || (!h) || (!s) )
00329 {
00330 cleanup();
00331 <span class="keywordflow">if</span> (h) <span class="keyword">delete</span>[] h;
00332 <span class="keywordflow">if</span> (s) <span class="keyword">delete</span>[] s;
00333 <span class="keywordflow">if</span> (d) <span class="keyword">delete</span>[] d;
00334 <span class="keywordflow">return</span> Qwt::ErrNoMem;
00335 }
00336
00337 <span class="comment">//</span>
00338 <span class="comment">// setup equation system; use coefficient</span>
00339 <span class="comment">// vectors as temporary buffers</span>
00340 <span class="comment">//</span>
00341 <span class="keywordflow">for</span> (i=0; i<d_size - 1; i++)
00342 {
00343 h[i] = d_x[i+1] - d_x[i];
00344 <span class="keywordflow">if</span> (h[i] <= 0.0)
00345 {
00346 <span class="keyword">delete</span>[] h;
00347 <span class="keyword">delete</span>[] s;
00348 <span class="keyword">delete</span>[] d;
00349 <span class="keywordflow">return</span> Qwt::ErrMono;
00350 }
00351 }
00352
00353 imax = d_size - 2;
00354 htmp = h[imax];
00355 dy1 = (d_y[0] - d_y[imax]) / htmp;
00356 <span class="keywordflow">for</span> (i=0; i <= imax; i++)
00357 {
00358 d_b[i] = d_c[i] = h[i];
00359 d_a[i] = 2.0 * (htmp + h[i]);
00360 dy2 = (d_y[i+1] - d_y[i]) / h[i];
00361 d[i] = 6.0 * ( dy1 - dy2);
00362 dy1 = dy2;
00363 htmp = h[i];
00364 }
00365
00366 <span class="comment">//</span>
00367 <span class="comment">// solve it</span>
00368 <span class="comment">//</span>
00369
00370 <span class="comment">// L-U Factorization</span>
00371 d_a[0] = sqrt(d_a[0]);
00372 d_c[0] = h[imax] / d_a[0];
00373 sum = 0;
00374
00375 <span class="keywordflow">for</span>(i=0;i<imax-1;i++)
00376 {
00377 d_b[i] /= d_a[i];
00378 <span class="keywordflow">if</span> (i > 0)
00379 d_c[i] = - d_c[i-1] * d_b[i-1] / d_a[i];
00380 d_a[i+1] = sqrt( d_a[i+1] - <a class="code" href="qwt__math_8h.html#a31">qwtSqr</a>(d_b[i]));
00381 sum += <a class="code" href="qwt__math_8h.html#a31">qwtSqr</a>(d_c[i]);
00382 }
00383 d_b[imax-1] = (d_b[imax-1] - d_c[imax-2] * d_b[imax-2]) / d_a[imax-1];
00384 d_a[imax] = sqrt(d_a[imax] - <a class="code" href="qwt__math_8h.html#a31">qwtSqr</a>(d_b[imax-1]) - sum);
00385
00386
00387 <span class="comment">// forward elimination</span>
00388 s[0] = d[0] / d_a[0];
00389 sum = 0;
00390 <span class="keywordflow">for</span>(i=1;i<imax;i++)
00391 {
00392 s[i] = (d[i] - d_b[i-1] * s[i-1]) / d_a[i];
00393 sum += d_c[i-1] * s[i-1];
00394 }
00395 s[imax] = (d[imax] - d_b[imax-1]*s[imax-1] - sum) / d_a[imax];
00396
00397
00398 <span class="comment">// backward elimination</span>
00399 s[imax] = - s[imax] / d_a[imax];
00400 s[imax-1] = -(s[imax-1] + d_b[imax-1] * s[imax]) / d_a[imax-1];
00401 <span class="keywordflow">for</span> (i= imax - 2; i >= 0; i--)
00402 s[i] = - (s[i] + d_b[i] * s[i+1] + d_c[i] * s[imax]) / d_a[i];
00403
00404 <span class="comment">//</span>
00405 <span class="comment">// Finally, determine the spline coefficients</span>
00406 <span class="comment">//</span>
00407 s[d_size-1] = s[0];
00408 <span class="keywordflow">for</span> (i=0;i<d_size-1;i++)
00409 {
00410 d_a[i] = ( s[i+1] - s[i] ) / ( 6.0 * h[i]);
00411 d_b[i] = 0.5 * s[i];
00412 d_c[i] = ( d_y[i+1] - d_y[i] )
00413 / h[i] - (s[i+1] + 2.0 * s[i] ) * h[i] / 6.0;
00414 }
00415
00416 <span class="keyword">delete</span>[] d;
00417 <span class="keyword">delete</span>[] s;
00418 <span class="keyword">delete</span>[] h;
00419
00420 <span class="keywordflow">return</span> 0;
00421 }
00422
00423
00425 <span class="keywordtype">void</span> QwtSpline::cleanup()
00426 {
00427 <span class="keywordflow">if</span> (d_a) <span class="keyword">delete</span>[] d_a;
00428 <span class="keywordflow">if</span> (d_b) <span class="keyword">delete</span>[] d_b;
00429 <span class="keywordflow">if</span> (d_c) <span class="keyword">delete</span>[] d_c;
00430 <span class="keywordflow">if</span> (d_xbuffer) <span class="keyword">delete</span>[] d_xbuffer;
00431 <span class="keywordflow">if</span> (d_ybuffer) <span class="keyword">delete</span>[] d_ybuffer;
00432 d_a = d_b = d_c = 0;
00433 d_xbuffer = d_ybuffer = d_x = d_y = 0;
00434 d_size = 0;
00435 }
</div></pre><hr size="1"><address style="align: right;"><small>Generated on Tue Nov 16 21:12:21 2004 for Qwt User's Guide by
<a href="http://www.doxygen.org/index.html">
<img src="doxygen.png" alt="doxygen" align="middle" border=0 ></a> 1.3.8 </small></address>
</body>
</html>
|
