1 #include "CurveTest.h"
2 #include "evoral/ControlList.h"
3 #include "evoral/Curve.h"
4 #include <stdlib.h>
5
6 CPPUNIT_TEST_SUITE_REGISTRATION (CurveTest);
7
8 #if defined(PLATFORM_WINDOWS) && defined(COMPILER_MINGW)
9 /* cppunit-1.13.2 uses assertion_traits<double>
10 * sprintf( , "%.*g", precision, x)
11 * to format a double. The actual comparison is performed on a string.
12 * This is problematic with mingw/windows|wine, "%.*g" formatting fails.
13 *
14 * This quick hack compares float, however float compatisons are at most Y.MMMM+eXX,
15 * the max precision needs to be limited. to the last mantissa digit.
16 *
17 * Anyway, actual maths is verified with Linux and OSX unit-tests,
18 * and this needs to go to https://sourceforge.net/p/cppunit/bugs/
19 */
20 #define MAXPREC(P) ((P) < .0005 ? .0005 : (P))
21 #define CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE(M,A,B,P) CPPUNIT_ASSERT_EQUAL_MESSAGE(M, (float)rint ((A) / MAXPREC(P)),(float)rint ((B) / MAXPREC(P)))
22 #define CPPUNIT_ASSERT_DOUBLES_EQUAL(A,B,P) CPPUNIT_ASSERT_EQUAL((float)rint ((A) / MAXPREC(P)),(float)rint ((B) / MAXPREC(P)))
23 #endif
24
25 using namespace Evoral;
26
27 // linear y = Y0 + YS * x ; with x = i * (X1 - X0) + X0; and i = [0..1023]
28 #define VEC1024LINCMP(X0, X1, Y0, YS) \
29 cl->curve ().get_vector ((X0), (X1), vec, 1024); \
30 for (int i = 0; i < 1024; ++i) { \
31 char msg[64]; \
32 snprintf (msg, 64, "at i=%d (x0=%.1f, x1=%.1f, y0=%.1f, ys=%.3f)", \
33 i, X0, X1, Y0, YS); \
34 CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE ( \
35 msg, \
36 (Y0) + i * (YS), vec[i], \
37 1e-24 \
38 ); \
39 }
40
41 void
trivial()42 CurveTest::trivial ()
43 {
44 float vec[1024];
45
46 boost::shared_ptr<Evoral::ControlList> cl = TestCtrlList();
47
48 cl->create_curve ();
49
50 // Empty curve
51 cl->curve().get_vector (1024.0, 2047.0, vec, 1024);
52 for (int i = 0; i < 1024; ++i) {
53 CPPUNIT_ASSERT_EQUAL (0.0f, vec[i]);
54 }
55
56 // Single point curve
57 cl->fast_simple_add(0.0, 42.0);
58 cl->curve().get_vector (1024.0, 2047.0, vec, 1024);
59 for (int i = 0; i < 1024; ++i) {
60 CPPUNIT_ASSERT_EQUAL (42.0f, vec[i]);
61 }
62 }
63
64 void
rtGet()65 CurveTest::rtGet ()
66 {
67 float vec[1024];
68
69 // Create simple control list
70 boost::shared_ptr<Evoral::ControlList> cl = TestCtrlList();
71 cl->create_curve ();
72 cl->fast_simple_add(0.0, 42.0);
73
74 {
75 // Write-lock list
76 Glib::Threads::RWLock::WriterLock lm(cl->lock());
77
78 // Attempt to get vector in RT (expect failure)
79 CPPUNIT_ASSERT (!cl->curve().rt_safe_get_vector (1024.0, 2047.0, vec, 1024));
80 }
81
82 // Attempt to get vector in RT (expect success)
83 CPPUNIT_ASSERT (cl->curve().rt_safe_get_vector (1024.0, 2047.0, vec, 1024));
84 for (int i = 0; i < 1024; ++i) {
85 CPPUNIT_ASSERT_EQUAL (42.0f, vec[i]);
86 }
87 }
88
89 void
twoPointLinear()90 CurveTest::twoPointLinear ()
91 {
92 float vec[1024];
93
94 boost::shared_ptr<Evoral::ControlList> cl = TestCtrlList();
95
96 cl->create_curve ();
97 cl->set_interpolation (ControlList::Linear);
98
99 // add two points to curve
100 cl->fast_simple_add ( 0.0 , 2048.0);
101 cl->fast_simple_add (8192.0 , 4096.0);
102
103 cl->curve ().get_vector (1024.0, 2047.0, vec, 1024);
104
105 VEC1024LINCMP (1024.0, 2047.0, 2304.f, .25f);
106 VEC1024LINCMP (2048.0, 2559.5, 2560.f, .125f);
107 VEC1024LINCMP ( 0.0, 4092.0, 2048.f, 1.f);
108
109 // greetings to tartina
110 cl->curve ().get_vector (2048.0, 2048.0, vec, 1);
111 CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=1 @ 2048..2048", 2560.f, vec[0]);
112
113 /* XXX WHAT DO WE EXPECT WITH veclen=1 AND x1 > x0 ? */
114 #if 0
115 /* .. interpolated value at (x1+x0)/2 */
116 cl->curve ().get_vector (2048.0, 2049.0, vec, 1);
117 CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=1 @ 2048-2049", 2560.125f, vec[0]);
118
119 cl->curve ().get_vector (2048.0, 2056.0, vec, 1);
120 CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=1 @ 2048-2049", 2561.f, vec[0]);
121 #else
122 /* .. value at x0 */
123 cl->curve ().get_vector (2048.0, 2049.0, vec, 1);
124 CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=1 , 2048..2049", 2560.f, vec[0]);
125
126 cl->curve ().get_vector (2048.0, 2056.0, vec, 1);
127 CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=1 , 2048..2049", 2560.f, vec[0]);
128 #endif
129
130 cl->curve ().get_vector (2048.0, 2048.0, vec, 2);
131 CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=2 , 2048..2048 @ 0", 2560.f, vec[0]);
132 CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=2 , 2048..2048 @ 1", 2560.f, vec[1]);
133
134 cl->curve ().get_vector (2048.0, 2056.0, vec, 2);
135 CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=2 , 2048..2056 @ 0", 2560.f, vec[0]);
136 CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=2 , 2048..2056 @ 0", 2562.f, vec[1]);
137
138 cl->curve ().get_vector (2048.0, 2056.0, vec, 3);
139 CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=3 , 2048..2056 @ 0", 2560.f, vec[0]);
140 CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=3 , 2048..2056 @ 1", 2561.f, vec[1]);
141 CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=3 , 2048..2056 @ 2", 2562.f, vec[2]);
142
143 /* check out-of range..
144 * we expect the first and last value - no interpolation
145 */
146 cl->curve ().get_vector (-1, -1, vec, 1);
147 CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=1 @ -1", 2048.f, vec[0]);
148
149 cl->curve ().get_vector (9999.0, 9999.0, vec, 1);
150 CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=1 @ 9999", 4096.f, vec[0]);
151
152 cl->curve ().get_vector (-999.0, 0, vec, 13);
153 for (int i = 0; i < 13; ++i) {
154 CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=13 @ -999..0", 2048.f, vec[i]);
155 }
156
157 cl->curve ().get_vector (9998.0, 9999.0, vec, 8);
158 for (int i = 0; i < 8; ++i) {
159 CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=8 @ 9998..9999", 4096.f, vec[i]);
160 }
161 }
162
163 void
threePointLinear()164 CurveTest::threePointLinear ()
165 {
166 float vec[4];
167
168 boost::shared_ptr<Evoral::ControlList> cl = TestCtrlList();
169
170 cl->create_curve ();
171 cl->set_interpolation (ControlList::Linear);
172
173 // add 3 points to curve
174 cl->fast_simple_add ( 0.0 , 2.0);
175 cl->fast_simple_add ( 100.0 , 4.0);
176 cl->fast_simple_add ( 200.0 , 0.0);
177
178 cl->curve ().get_vector (50.0, 60.0, vec, 1);
179 CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=1 @ 50", 3.f, vec[0]);
180
181 cl->curve ().get_vector (100.0, 100.0, vec, 1);
182 CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=1 @ 100", 4.f, vec[0]);
183
184 cl->curve ().get_vector (150.0, 150.0, vec, 1);
185 CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=1 @ 150", 2.f, vec[0]);
186
187 cl->curve ().get_vector (130.0, 150.0, vec, 3);
188 CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=3 130..150 @ 0", 2.8f, vec[0]);
189 CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=3 130..150 @ 2", 2.4f, vec[1]);
190 CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=3 130..150 @ 3", 2.0f, vec[2]);
191
192 cl->curve ().get_vector (80.0, 160.0, vec, 3);
193 CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=3 80..160 @ 0", 3.6f, vec[0]);
194 CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=3 80..160 @ 2", 3.2f, vec[1]);
195 CPPUNIT_ASSERT_EQUAL_MESSAGE ("veclen=3 80..160 @ 3", 1.6f, vec[2]);
196 }
197
198 void
threePointDiscete()199 CurveTest::threePointDiscete ()
200 {
201 boost::shared_ptr<Evoral::ControlList> cl = TestCtrlList();
202 cl->set_interpolation (ControlList::Discrete);
203
204 // add 3 points to curve
205 cl->fast_simple_add ( 0.0 , 2.0);
206 cl->fast_simple_add ( 100.0 , 4.0);
207 cl->fast_simple_add ( 200.0 , 0.0);
208
209 CPPUNIT_ASSERT_EQUAL(2.0, cl->unlocked_eval(80.));
210 CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(120.));
211 CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(160.));
212
213 cl->set_interpolation (ControlList::Linear);
214
215 CPPUNIT_ASSERT_EQUAL(3.6, cl->unlocked_eval(80.));
216 CPPUNIT_ASSERT_EQUAL(3.2, cl->unlocked_eval(120.));
217 CPPUNIT_ASSERT_EQUAL(1.6, cl->unlocked_eval(160.));
218 }
219
220 void
ctrlListEval()221 CurveTest::ctrlListEval ()
222 {
223 boost::shared_ptr<Evoral::ControlList> cl = TestCtrlList();
224
225 cl->fast_simple_add ( 0.0 , 2.0);
226
227 cl->set_interpolation (ControlList::Discrete);
228 CPPUNIT_ASSERT_EQUAL(2.0, cl->unlocked_eval(80.));
229 CPPUNIT_ASSERT_EQUAL(2.0, cl->unlocked_eval(120.));
230 CPPUNIT_ASSERT_EQUAL(2.0, cl->unlocked_eval(160.));
231
232 cl->set_interpolation (ControlList::Linear);
233 CPPUNIT_ASSERT_EQUAL(2.0, cl->unlocked_eval(80.));
234 CPPUNIT_ASSERT_EQUAL(2.0, cl->unlocked_eval(120.));
235 CPPUNIT_ASSERT_EQUAL(2.0, cl->unlocked_eval(160.));
236
237 cl->fast_simple_add ( 100.0 , 4.0);
238
239 cl->set_interpolation (ControlList::Discrete);
240 CPPUNIT_ASSERT_EQUAL(2.0, cl->unlocked_eval(80.));
241 CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(120.));
242 CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(160.));
243
244 cl->set_interpolation (ControlList::Linear);
245 CPPUNIT_ASSERT_EQUAL(3.6, cl->unlocked_eval(80.));
246 CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(120.));
247 CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(160.));
248
249 cl->fast_simple_add ( 200.0 , 0.0);
250
251 cl->set_interpolation (ControlList::Discrete);
252 CPPUNIT_ASSERT_EQUAL(2.0, cl->unlocked_eval(80.));
253 CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(120.));
254 CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(160.));
255
256 cl->set_interpolation (ControlList::Linear);
257 CPPUNIT_ASSERT_EQUAL(3.6, cl->unlocked_eval(80.));
258 CPPUNIT_ASSERT_EQUAL(3.2, cl->unlocked_eval(120.));
259 CPPUNIT_ASSERT_EQUAL(1.6, cl->unlocked_eval(160.));
260
261 cl->fast_simple_add ( 300.0 , 8.0);
262
263 cl->set_interpolation (ControlList::Discrete);
264 CPPUNIT_ASSERT_EQUAL(2.0, cl->unlocked_eval(80.));
265 CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(120.));
266 CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(160.));
267 CPPUNIT_ASSERT_EQUAL(0.0, cl->unlocked_eval(250.));
268 CPPUNIT_ASSERT_EQUAL(8.0, cl->unlocked_eval(999.));
269
270 cl->set_interpolation (ControlList::Linear);
271 CPPUNIT_ASSERT_EQUAL(3.6, cl->unlocked_eval(80.));
272 CPPUNIT_ASSERT_EQUAL(3.2, cl->unlocked_eval(120.));
273 CPPUNIT_ASSERT_EQUAL(1.6, cl->unlocked_eval(160.));
274 CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(250.));
275 CPPUNIT_ASSERT_EQUAL(8.0, cl->unlocked_eval(999.));
276
277 cl->fast_simple_add ( 400.0 , 9.0);
278
279 cl->set_interpolation (ControlList::Discrete);
280 CPPUNIT_ASSERT_EQUAL(2.0, cl->unlocked_eval(80.));
281 CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(120.));
282 CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(160.));
283 CPPUNIT_ASSERT_EQUAL(0.0, cl->unlocked_eval(250.));
284 CPPUNIT_ASSERT_EQUAL(8.0, cl->unlocked_eval(350.));
285 CPPUNIT_ASSERT_EQUAL(9.0, cl->unlocked_eval(999.));
286
287 cl->set_interpolation (ControlList::Linear);
288 CPPUNIT_ASSERT_EQUAL(3.6, cl->unlocked_eval(80.));
289 CPPUNIT_ASSERT_EQUAL(3.2, cl->unlocked_eval(120.));
290 CPPUNIT_ASSERT_EQUAL(1.6, cl->unlocked_eval(160.));
291 CPPUNIT_ASSERT_EQUAL(4.0, cl->unlocked_eval(250.));
292 CPPUNIT_ASSERT_EQUAL(8.5, cl->unlocked_eval(350.));
293 CPPUNIT_ASSERT_EQUAL(9.0, cl->unlocked_eval(999.));
294 }
295
296 void
constrainedCubic()297 CurveTest::constrainedCubic ()
298 {
299
300 struct point {
301 int x, y;
302 };
303
304 static const struct point data[] = {
305 /* values from worked example in www.korf.co.uk/spline.pdf */
306 { 0, 30 },
307 { 10, 130 },
308 { 30, 150 },
309 { 50, 150 },
310 { 70, 170 },
311 { 90, 220 },
312 { 100, 320 },
313 };
314
315 int32_t type = 0;
316 Evoral::Parameter p(type);
317 Evoral::ParameterDescriptor pd;
318 pd.lower = 5;
319 pd.upper = 325;
320 Evoral::ControlList l(p,pd);
321
322 size_t i;
323 l.set_interpolation(Evoral::ControlList::Curved);
324
325 for (i=0; i<sizeof(data)/sizeof(data[0]); i++) {
326 l.add (data[i].x, data[i].y);
327 }
328
329 Evoral::Curve curve(l);
330
331 float f[121];
332 curve.get_vector(-10, 110, f, 121);
333
334 const float *g = &f[10]; /* so g starts at x==0 */
335
336 /* given points - should be exactly equal */
337 CPPUNIT_ASSERT_EQUAL( 30.0f, g[-10]);
338 CPPUNIT_ASSERT_EQUAL( 30.0f, g[ 0]);
339 CPPUNIT_ASSERT_EQUAL(130.0f, g[ 10]);
340 CPPUNIT_ASSERT_EQUAL(150.0f, g[ 30]);
341 CPPUNIT_ASSERT_EQUAL(150.0f, g[ 40]);
342 CPPUNIT_ASSERT_EQUAL(150.0f, g[ 50]);
343 CPPUNIT_ASSERT_EQUAL(320.0f, g[100]);
344 CPPUNIT_ASSERT_EQUAL(320.0f, g[110]);
345
346 /*
347 First segment, i=1, for 0 <= x <= 10
348 f'1(x1) = 2/((x2 – x1)/(y2 – y1) + (x1 – x0)/(y1 – y0))
349 = 2/((30 – 10)/(150 – 130) + (10 – 0)/(130 – 30))
350 = 1.8181
351 f'1(x0) = 3/2*(y1 – y0)/(x1 – x0) - f'1(x1)/2
352 = 3/2*(130 – 30)/(10 – 0) – 1.818/2
353 = 14.0909
354 f"1(x0) = -2*(f'1(x1) + 2* f'1(x0))/(x1 – x0) + 6*(y1 – y0)/ (x1 – x0)^2
355 = -2*(1.8181 + 2*14.0909)/(10 – 0) + 6*(130 – 30)/(10 – 0)^2
356 = 0
357 f"1(x1) = 2*(2*f'1(x1) + f'1(x0))/(x1 – x0) - 6*(y1 – y0)/ (x1 – x0)^2
358 = 2*(2*1.818 + 14.0909)/(10 – 0) – 6*(130 – 30)/(10 – 0)^2
359 = -2.4545
360 d1 = 1/6 * (f"1(x1) - f"1(x0))/(x1 – x0)
361 = 1/6 * (-2.4545 – 0)/(10 – 0)
362 = -0.0409
363 c1 = 1/2 * (x1*f"1(x0) – x0*f"1(x1))/(x1 – x0)
364 = 1/2 * (10*0 – 0*1.8181)/(10 – 0)
365 = 0
366 b1 = ((y1 – y0) – c1*(x21 – x20) – d1*( x31 – x30))/(x1 – x0)
367 = ((130 – 30) – 0*(102 – 02) + 0.0409*(103 – 03))/(10 – 0)
368 = 14.09
369 a1 = y0 – b1*x0 – c1*x20 – d1*x30
370 = 30
371 y1 = 30 + 14.09x - 0.0409x3 for 0 <= x <= 10
372 */
373 /*
374 Second segment, i=2, for 10 <= x <= 30
375 f'2(x2) = 2/((x3 – x2)/(y3 – y2) + (x2 – x1)/(y2 – y1))
376 = 2/((50 – 30)/(150 – 150) + (30 – 10)/(150 – 130))
377 = 0
378 f'2(x1) = 2/((x2 – x1)/(y2 – y1) + (x1 – x0)/(y1 – y0))
379 = 1.8181
380
381 f"2(x1) = -2*(f'2(x2) + 2* f'2(x1))/(x2 – x1) + 6*(y2 – y1)/ (x2 – x1)^2
382 = -2*(0 + 2*1.8181)/(30 – 10) + 6*(150 – 130)/(30 – 10)2
383 = -0.063636
384 f"2(x2) = 2*(2*f'2(x2) + f'2(x1))/(x2 – x1) - 6*(y2 – y1)/ (x2 – x1)^2
385 = 2*(2*0 + 1.8181)/(30 – 10) – 6*(150 – 130)/(30 – 10)^2
386 = -0.11818
387
388 d2 = 1/6 * (f"2(x2) - f"2(x1))/(x2 – x1)
389 = 1/6 * (-0.11818 + 0.063636)/(30 – 10)
390 = -0.0004545
391 c2 = 1/2 * (x2*f"2(x1) – x1*f"2(x2))/(x2 – x1)
392 = 1/2 * (-30*0.063636 + 10*0.11818)/(30 – 10)
393 = -0.01818
394 b2 = ((y2 – y1) – c2*(x2^2 – x1^2) – d2*( x2^3 – x1^3))/(x2 – x1)
395 = ((150 – 130) + 0.01818*(302 – 102) + 0.0004545*(303 – 103))/(30 – 10)
396 = 2.31818
397 a2 = y1 – b2*x1 – c2*x1^2 – d2*x1^3
398 = 130 – 2.31818*10 + 0.01818*102 + 0.0004545*103
399 = 109.09
400 y2 = 109.09 + 2.31818x - 0.01818x^2 - 0.0004545x^3 for 10 <= x <= 30
401 */
402
403
404 int x;
405 long double a1, b1, c1, d1, a2, b2, c2, d2, fdx0, fddx0, fdx1, fdx2, fddx1, fddx2;
406 double x0 = data[0].x;
407 double y0 = data[0].y;
408 double x1 = data[1].x;
409 double y1 = data[1].y;
410 double x2 = data[2].x;
411 double y2 = data[2].y;
412 double x3 = data[3].x;
413 double y3 = data[3].y;
414
415 double dx0 = x1 - x0;
416 double dy0 = y1 - y0;
417 double dx1 = x2 - x1;
418 double dy1 = y2 - y1;
419 double dx2 = x3 - x2;
420 double dy2 = y3 - y2;
421
422 // First (leftmost) segment
423 fdx1 = 2.0 / ( dx1 / dy1 + dx0 / dy0 );
424 fdx0 = 3.0 / 2.0 * dy0 / dx0 - fdx1 / 2.0;
425
426 fddx0 = -2.0 * (fdx1 + 2.0 * fdx0) / dx0 + 6.0 * dy0 / (dx0*dx0);
427 fddx1 = 2.0 * (2.0 * fdx1 + fdx0) / dx0 - 6.0 * dy0 / (dx0*dx0);
428 d1 = 1.0 / 6.0 * (fddx1 - fddx0) / dx0;
429 c1 = 1.0 / 2.0 * (x1 * fddx0 - x0 * fddx1) / dx0;
430 b1 = (dy0 - c1 * (x1* x1 - x0*x0) - d1 * (x1*x1*x1 - x0*x0*x0)) / dx0;
431 a1 = y0 - b1*x0 - c1*x0*x0 - d1*x0*x0*x0;
432
433 // printf("dx0=%f, dy0=%f, dx1=%f, dy1=%f\n", dx0, dy0, dx1, dy1);
434 // printf("fdx0=%Lf, fdx1=%Lf, fddx0=%Lf, fddx1=%Lf\n", fdx0, fdx1, fddx0, fddx1);
435 // printf("a1=%Lf, b1=%Lf, c1=%Lf, d1=%Lf\n", a1, b1, c1, d1);
436
437 // values from worked example: deltas rather arbitrary, I'm afraid
438 CPPUNIT_ASSERT_DOUBLES_EQUAL(30.0, a1, 0.1);
439 CPPUNIT_ASSERT_DOUBLES_EQUAL(14.09, b1, 0.01);
440 CPPUNIT_ASSERT_DOUBLES_EQUAL(0.0, c1, 0.1);
441 CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.0409, d1, 0.0001);
442
443 for (x = 0; x <= 10; x++) {
444 double v = a1 + b1*x + c1*x*x + d1*x*x*x;
445 char msg[64];
446 snprintf(msg, 64, "interpolating %d: v=%f, x=%f...\n", x, v, g[x]);
447 CPPUNIT_ASSERT_DOUBLES_EQUAL(v, g[x], 0.000004);
448 }
449
450 // Second segment
451 fdx2 = 2.0 / ( dx2 / dy2 + dx1 / dy1 );
452
453 fddx1 = -2.0 * (fdx2 + 2.0 * fdx1) / dx1 + 6.0 * dy1 / (dx1*dx1);
454 fddx2 = 2.0 * (2.0 * fdx2 + fdx1) / dx1 - 6.0 * dy1 / (dx1*dx1);
455 d2 = 1.0 / 6.0 * (fddx2 - fddx1) / dx1;
456 c2 = 1.0 / 2.0 * (x2 * fddx1 - x1 * fddx2) / dx1;
457 b2 = (dy1 - c2 * (x2*x2 - x1*x1) - d2 * (x2*x2*x2 - x1*x1*x1)) / dx1;
458 a2 = y1 - b2*x1 - c2*x1*x1 - d2*x1*x1*x1;
459
460 // printf("dx0=%f, dy0=%f, dx1=%f, dy1=%f dx2=%f, dy2=%f\n", dx0, dy0, dx1, dy1, dx2, dy2);
461 // printf("fdx1=%Lf, fdx2=%Lf, fddx1=%Lf, fddx2=%Lf\n", fdx1, fdx2, fddx1, fddx2);
462 // printf("a2=%Lf, b2=%Lf, c2=%Lf, d2=%Lf\n", a2, b2, c2, d2);
463
464 // values from worked example: deltas rather arbitrary, I'm afraid
465 CPPUNIT_ASSERT_DOUBLES_EQUAL(109.09, a2, 0.01);
466 CPPUNIT_ASSERT_DOUBLES_EQUAL(2.31818, b2, 0.00001);
467 CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.01818, c2, 0.00001);
468 CPPUNIT_ASSERT_DOUBLES_EQUAL(-0.0004545, d2, 0.0000001);
469
470 for (x = 10; x <= 30; x++) {
471 double v = a2 + b2*x + c2*x*x + d2*x*x*x;
472 char msg[64];
473 snprintf(msg, 64, "interpolating %d: v=%f, x=%f...\n", x, v, g[x]);
474 CPPUNIT_ASSERT_DOUBLES_EQUAL(v, g[x], 0.000008);
475 }
476 }
477