Branch data Line data Source code
1 : : /* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2 : : /*************************************************************************
3 : : *
4 : : * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
5 : : *
6 : : * Copyright 2000, 2010 Oracle and/or its affiliates.
7 : : *
8 : : * OpenOffice.org - a multi-platform office productivity suite
9 : : *
10 : : * This file is part of OpenOffice.org.
11 : : *
12 : : * OpenOffice.org is free software: you can redistribute it and/or modify
13 : : * it under the terms of the GNU Lesser General Public License version 3
14 : : * only, as published by the Free Software Foundation.
15 : : *
16 : : * OpenOffice.org is distributed in the hope that it will be useful,
17 : : * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 : : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 : : * GNU Lesser General Public License version 3 for more details
20 : : * (a copy is included in the LICENSE file that accompanied this code).
21 : : *
22 : : * You should have received a copy of the GNU Lesser General Public License
23 : : * version 3 along with OpenOffice.org. If not, see
24 : : * <http://www.openoffice.org/license.html>
25 : : * for a copy of the LGPLv3 License.
26 : : *
27 : : ************************************************************************/
28 : :
29 : :
30 : : #include <math.h>
31 : :
32 : :
33 : : #include <tools/poly.hxx>
34 : :
35 : : extern "C" {
36 : :
37 : : /*.pn 277 */
38 : : /*.hlAnhang: C - Programme*/
39 : : /*.hrKonstanten- und Macro-Definitionen*/
40 : : /*.fe Die Include-Datei u_const.h ist in das Verzeichnis zu stellen, */
41 : : /*.fe wo der Compiler nach Include-Dateien sucht. */
42 : :
43 : :
44 : : /*----------------------- FILE u_const.h ---------------------------*/
45 : :
46 : : #define IEEE
47 : :
48 : : /* IEEE - Norm fuer die Darstellung von Gleitkommazahlen:
49 : :
50 : : 8 Byte lange Gleitkommazahlen, mit
51 : :
52 : : 53 Bit Mantisse ==> Mantissenbereich: 2 hoch 52 versch. Zahlen
53 : : mit 0.1 <= Zahl < 1.0,
54 : : 1 Vorzeichen-Bit
55 : : 11 Bit Exponent ==> Exponentenbereich: -1024...+1023
56 : :
57 : : Die 1. Zeile ( #define IEEE ) ist zu loeschen, falls die Maschine
58 : : bzw. der Compiler keine Gleitpunktzahlen gemaess der IEEE-Norm
59 : : benutzt. Zusaetzlich muessen die Zahlen MAXEXPON, MINEXPON
60 : : (s.u.) angepasst werden.
61 : : */
62 : :
63 : : #ifdef IEEE /*----------- Falls IEEE Norm --------------------*/
64 : :
65 : : #define MACH_EPS 2.220446049250313e-016 /* Maschinengenauigkeit */
66 : : /* IBM-AT: = 2 hoch -52 */
67 : : /* MACH_EPS ist die kleinste positive, auf der Maschine darstellbare
68 : : Zahl x, die der Bedingung genuegt: 1.0 + x > 1.0 */
69 : :
70 : : #define EPSQUAD 4.930380657631324e-032
71 : : #define EPSROOT 1.490116119384766e-008
72 : :
73 : : #define POSMAX 8.98846567431158e+307 /* groesste positive Zahl */
74 : : #define POSMIN 5.56268464626800e-309 /* kleinste positive Zahl */
75 : : #define MAXROOT 9.48075190810918e+153
76 : :
77 : : #define BASIS 2 /* Basis der Zahlendarst. */
78 : : #ifndef PI
79 : : #define PI 3.141592653589793e+000
80 : : #endif
81 : : #define EXP_1 2.718281828459045e+000
82 : :
83 : : #else /*------------------ sonst -----------------------*/
84 : :
85 : : double exp (double);
86 : : double atan (double);
87 : : double pow (double,double);
88 : : double sqrt (double);
89 : :
90 : : double masch() /* MACH_EPS maschinenunabhaengig bestimmen */
91 : : {
92 : : double eps = 1.0, x = 2.0, y = 1.0;
93 : : while ( y < x )
94 : : { eps *= 0.5;
95 : : x = 1.0 + eps;
96 : : }
97 : : eps *= 2.0; return (eps);
98 : : }
99 : :
100 : : short basis() /* BASIS maschinenunabhaengig bestimmen */
101 : : {
102 : : double x = 1.0, one = 1.0, b = 1.0;
103 : :
104 : : while ( (x + one) - x == one ) x *= 2.0;
105 : : while ( (x + b) == x ) b *= 2.0;
106 : :
107 : : return ( (short) ((x + b) - x) );
108 : : }
109 : :
110 : : #define BASIS basis() /* Basis der Zahlendarst. */
111 : :
112 : : /* Falls die Maschine (der Compiler) keine IEEE-Darstellung fuer
113 : : Gleitkommazahlen nutzt, muessen die folgenden 2 Konstanten an-
114 : : gepasst werden.
115 : : */
116 : :
117 : : #define MAXEXPON 1023.0 /* groesster Exponent */
118 : : #define MINEXPON -1024.0 /* kleinster Exponent */
119 : :
120 : :
121 : : #define MACH_EPS masch()
122 : : #define EPSQUAD MACH_EPS * MACH_EPS
123 : : #define EPSROOT sqrt(MACH_EPS)
124 : :
125 : : #define POSMAX pow ((double) BASIS, MAXEXPON)
126 : : #define POSMIN pow ((double) BASIS, MINEXPON)
127 : : #define MAXROOT sqrt(POSMAX)
128 : :
129 : : #define PI 4.0 * atan (1.0)
130 : : #define EXP_1 exp(1.0)
131 : :
132 : : #endif /*-------------- ENDE ifdef ----------------------*/
133 : :
134 : :
135 : : #define NEGMAX -POSMIN /* groesste negative Zahl */
136 : : #define NEGMIN -POSMAX /* kleinste negative Zahl */
137 : :
138 : : /* Definition von Funktionsmakros:
139 : : */
140 : :
141 : : #define abs(X) ((X) >= 0 ? (X) : -(X)) /* Absolutbetrag von X */
142 : : #define sign(X, Y) (Y < 0 ? -abs(X) : abs(X)) /* Vorzeichen von */
143 : : /* Y mal abs(X) */
144 : : #define sqr(X) ((X) * (X)) /* Quadrat von X */
145 : :
146 : : /*------------------- ENDE FILE u_const.h --------------------------*/
147 : :
148 : :
149 : :
150 : :
151 : :
152 : :
153 : :
154 : :
155 : :
156 : : /*.HL Anhang: C - Programme*/
157 : : /*.HRGleichungssysteme fuer Tridiagonalmatrizen*/
158 : :
159 : : /*.FE P 3.7 TRIDIAGONALE GLEICHUNGSSYSTEME*/
160 : :
161 : :
162 : : /*---------------------- MODUL TRIDIAGONAL ------------------------*/
163 : :
164 : 0 : sal_uInt16 TriDiagGS(sal_Bool rep, sal_uInt16 n, double* lower,
165 : : double* diag, double* upper, double* b)
166 : : /************************/
167 : : /* GAUSS-Verfahren fuer */
168 : : /* Tridiagonalmatrizen */
169 : : /************************/
170 : :
171 : : /*====================================================================*/
172 : : /* */
173 : : /* trdiag bestimmt die Loesung x des linearen Gleichungssystems */
174 : : /* A * x = b mit tridiagonaler n x n Koeffizientenmatrix A, die in */
175 : : /* den 3 Vektoren lower, upper und diag wie folgt abgespeichert ist: */
176 : : /* */
177 : : /* ( diag[0] upper[0] 0 0 . . . 0 ) */
178 : : /* ( lower[1] diag[1] upper[1] 0 . . . ) */
179 : : /* ( 0 lower[2] diag[2] upper[2] 0 . ) */
180 : : /* A = ( . 0 lower[3] . . . ) */
181 : : /* ( . . . . . 0 ) */
182 : : /* ( . . . . . ) */
183 : : /* ( . . . upper[n-2] ) */
184 : : /* ( 0 . . . 0 lower[n-1] diag[n-1] ) */
185 : : /* */
186 : : /*====================================================================*/
187 : : /* */
188 : : /* Anwendung: */
189 : : /* ========= */
190 : : /* Vorwiegend fuer diagonaldominante Tridiagonalmatrizen, wie */
191 : : /* sie bei der Spline-Interpolation auftreten. */
192 : : /* Fuer diagonaldominante Matrizen existiert immer eine LU- */
193 : : /* Zerlegung; fuer nicht diagonaldominante Tridiagonalmatrizen */
194 : : /* sollte die Funktion band vorgezogen werden, da diese mit */
195 : : /* Spaltenpivotsuche arbeitet und daher numerisch stabiler ist. */
196 : : /* */
197 : : /*====================================================================*/
198 : : /* */
199 : : /* Eingabeparameter: */
200 : : /* ================ */
201 : : /* n Dimension der Matrix ( > 1 ) sal_uInt16 n */
202 : : /* */
203 : : /* lower untere Nebendiagonale double lower[n] */
204 : : /* diag Hauptdiagonale double diag[n] */
205 : : /* upper obere Nebendiagonale double upper[n] */
206 : : /* */
207 : : /* bei rep != 0 enthalten lower, diag und upper die */
208 : : /* Dreieckzerlegung der Ausgangsmatrix. */
209 : : /* */
210 : : /* b rechte Seite des Systems double b[n] */
211 : : /* rep = 0 erstmaliger Aufruf sal_Bool rep */
212 : : /* !=0 wiederholter Aufruf */
213 : : /* fuer gleiche Matrix, */
214 : : /* aber verschiedenes b. */
215 : : /* */
216 : : /* Ausgabeparameter: */
217 : : /* ================ */
218 : : /* b Loesungsvektor des Systems; double b[n] */
219 : : /* die urspruengliche rechte Seite wird ueberspeichert */
220 : : /* */
221 : : /* lower ) enthalten bei rep = 0 die Zerlegung der Matrix; */
222 : : /* diag ) die urspruenglichen Werte von lower u. diag werden */
223 : : /* upper ) ueberschrieben */
224 : : /* */
225 : : /* Die Determinante der Matrix ist bei rep = 0 durch */
226 : : /* det A = diag[0] * ... * diag[n-1] bestimmt. */
227 : : /* */
228 : : /* Rueckgabewert: */
229 : : /* ============= */
230 : : /* = 0 alles ok */
231 : : /* = 1 n < 2 gewaehlt */
232 : : /* = 2 Die Dreieckzerlegung der Matrix existiert nicht */
233 : : /* */
234 : : /*====================================================================*/
235 : : /* */
236 : : /* Benutzte Funktionen: */
237 : : /* =================== */
238 : : /* */
239 : : /* Aus der C Bibliothek: fabs() */
240 : : /* */
241 : : /*====================================================================*/
242 : :
243 : : /*.cp 5 */
244 : : {
245 : : sal_uInt16 i;
246 : : short j;
247 : :
248 : : // double fabs(double);
249 : :
250 [ # # ]: 0 : if ( n < 2 ) return(1); /* n mindestens 2 */
251 : :
252 : : /* Wenn rep = 0 ist, */
253 : : /* Dreieckzerlegung der */
254 [ # # ]: 0 : if (rep == 0) /* Matrix u. det be- */
255 : : { /* stimmen */
256 [ # # ]: 0 : for (i = 1; i < n; i++)
257 [ # # ]: 0 : { if ( fabs(diag[i-1]) < MACH_EPS ) /* Wenn ein diag[i] = 0 */
258 : 0 : return(2); /* ist, ex. keine Zerle- */
259 : 0 : lower[i] /= diag[i-1]; /* gung. */
260 : 0 : diag[i] -= lower[i] * upper[i-1];
261 : : }
262 : : }
263 : :
264 [ # # ]: 0 : if ( fabs(diag[n-1]) < MACH_EPS ) return(2);
265 : :
266 [ # # ]: 0 : for (i = 1; i < n; i++) /* Vorwaertselimination */
267 : 0 : b[i] -= lower[i] * b[i-1];
268 : :
269 : 0 : b[n-1] /= diag[n-1]; /* Rueckwaertselimination */
270 [ # # ]: 0 : for (j = n-2; j >= 0; j--) {
271 : 0 : i=j;
272 : 0 : b[i] = ( b[i] - upper[i] * b[i+1] ) / diag[i];
273 : : }
274 : 0 : return(0);
275 : : }
276 : :
277 : : /*----------------------- ENDE TRIDIAGONAL -------------------------*/
278 : :
279 : :
280 : :
281 : :
282 : :
283 : :
284 : :
285 : :
286 : :
287 : : /*.HL Anhang: C - Programme*/
288 : : /*.HRGleichungssysteme mit zyklisch tridiagonalen Matrizen*/
289 : :
290 : : /*.FE P 3.8 SYSTEME MIT ZYKLISCHEN TRIDIAGONALMATRIZEN */
291 : :
292 : :
293 : : /*---------------- MODUL ZYKLISCH TRIDIAGONAL ----------------------*/
294 : :
295 : :
296 : 30 : sal_uInt16 ZyklTriDiagGS(sal_Bool rep, sal_uInt16 n, double* lower, double* diag,
297 : : double* upper, double* lowrow, double* ricol, double* b)
298 : : /******************************/
299 : : /* Systeme mit zyklisch tri- */
300 : : /* diagonalen Matrizen */
301 : : /******************************/
302 : :
303 : : /*====================================================================*/
304 : : /* */
305 : : /* tzdiag bestimmt die Loesung x des linearen Gleichungssystems */
306 : : /* A * x = b mit zyklisch tridiagonaler n x n Koeffizienten- */
307 : : /* matrix A, die in den 5 Vektoren lower, upper, diag, lowrow und */
308 : : /* ricol wie folgt abgespeichert ist: */
309 : : /* */
310 : : /* ( diag[0] upper[0] 0 0 . . 0 ricol[0] ) */
311 : : /* ( lower[1] diag[1] upper[1] 0 . . 0 ) */
312 : : /* ( 0 lower[2] diag[2] upper[2] 0 . ) */
313 : : /* A = ( . 0 lower[3] . . . . ) */
314 : : /* ( . . . . . 0 ) */
315 : : /* ( . . . . . ) */
316 : : /* ( 0 . . . upper[n-2] ) */
317 : : /* ( lowrow[0] 0 . . 0 lower[n-1] diag[n-1] ) */
318 : : /* */
319 : : /* Speicherplatz fuer lowrow[1],..,lowrow[n-3] und ricol[1],..., */
320 : : /* ricol[n-3] muss zusaetzlich bereitgestellt werden, da dieser */
321 : : /* fuer die Aufnahme der Zerlegungsmatrix verfuegbar sein muss, die */
322 : : /* auf die 5 genannten Vektoren ueberspeichert wird. */
323 : : /* */
324 : : /*====================================================================*/
325 : : /* */
326 : : /* Anwendung: */
327 : : /* ========= */
328 : : /* Vorwiegend fuer diagonaldominante zyklische Tridiagonalmatri- */
329 : : /* zen wie sie bei der Spline-Interpolation auftreten. */
330 : : /* Fuer diagonaldominante Matrizen existiert immer eine LU- */
331 : : /* Zerlegung. */
332 : : /* */
333 : : /*====================================================================*/
334 : : /* */
335 : : /* Eingabeparameter: */
336 : : /* ================ */
337 : : /* n Dimension der Matrix ( > 2 ) sal_uInt16 n */
338 : : /* lower untere Nebendiagonale double lower[n] */
339 : : /* diag Hauptdiagonale double diag[n] */
340 : : /* upper obere Nebendiagonale double upper[n] */
341 : : /* b rechte Seite des Systems double b[n] */
342 : : /* rep = 0 erstmaliger Aufruf sal_Bool rep */
343 : : /* !=0 wiederholter Aufruf */
344 : : /* fuer gleiche Matrix, */
345 : : /* aber verschiedenes b. */
346 : : /* */
347 : : /* Ausgabeparameter: */
348 : : /* ================ */
349 : : /* b Loesungsvektor des Systems, double b[n] */
350 : : /* die urspruengliche rechte Seite wird ueberspeichert */
351 : : /* */
352 : : /* lower ) enthalten bei rep = 0 die Zerlegung der Matrix; */
353 : : /* diag ) die urspruenglichen Werte von lower u. diag werden */
354 : : /* upper ) ueberschrieben */
355 : : /* lowrow ) double lowrow[n-2] */
356 : : /* ricol ) double ricol[n-2] */
357 : : /* */
358 : : /* Die Determinante der Matrix ist bei rep = 0 durch */
359 : : /* det A = diag[0] * ... * diag[n-1] bestimmt. */
360 : : /* */
361 : : /* Rueckgabewert: */
362 : : /* ============= */
363 : : /* = 0 alles ok */
364 : : /* = 1 n < 3 gewaehlt */
365 : : /* = 2 Die Zerlegungsmatrix existiert nicht */
366 : : /* */
367 : : /*====================================================================*/
368 : : /* */
369 : : /* Benutzte Funktionen: */
370 : : /* =================== */
371 : : /* */
372 : : /* Aus der C Bibliothek: fabs() */
373 : : /* */
374 : : /*====================================================================*/
375 : :
376 : : /*.cp 5 */
377 : : {
378 : : double temp; // fabs(double);
379 : : sal_uInt16 i;
380 : : short j;
381 : :
382 [ - + ]: 30 : if ( n < 3 ) return(1);
383 : :
384 [ + - ]: 30 : if (rep == 0) /* Wenn rep = 0 ist, */
385 : : { /* Zerlegung der */
386 : 30 : lower[0] = upper[n-1] = 0.0; /* Matrix berechnen. */
387 : :
388 [ - + ]: 30 : if ( fabs (diag[0]) < MACH_EPS ) return(2);
389 : : /* Ist ein Diagonalelement */
390 : 30 : temp = 1.0 / diag[0]; /* betragsmaessig kleiner */
391 : 30 : upper[0] *= temp; /* MACH_EPS, so ex. keine */
392 : 30 : ricol[0] *= temp; /* Zerlegung. */
393 : :
394 [ + + ]: 330 : for (i = 1; i < n-2; i++)
395 : 300 : { diag[i] -= lower[i] * upper[i-1];
396 [ - + ]: 300 : if ( fabs(diag[i]) < MACH_EPS ) return(2);
397 : 300 : temp = 1.0 / diag[i];
398 : 300 : upper[i] *= temp;
399 : 300 : ricol[i] = -lower[i] * ricol[i-1] * temp;
400 : : }
401 : :
402 : 30 : diag[n-2] -= lower[n-2] * upper[n-3];
403 [ - + ]: 30 : if ( fabs(diag[n-2]) < MACH_EPS ) return(2);
404 : :
405 [ + + ]: 330 : for (i = 1; i < n-2; i++)
406 : 300 : lowrow[i] = -lowrow[i-1] * upper[i-1];
407 : :
408 : 30 : lower[n-1] -= lowrow[n-3] * upper[n-3];
409 : 30 : upper[n-2] = ( upper[n-2] - lower[n-2] * ricol[n-3] ) / diag[n-2];
410 : :
411 [ + + ]: 360 : for (temp = 0.0, i = 0; i < n-2; i++)
412 : 330 : temp -= lowrow[i] * ricol[i];
413 : 30 : diag[n-1] += temp - lower[n-1] * upper[n-2];
414 : :
415 [ - + ]: 30 : if ( fabs(diag[n-1]) < MACH_EPS ) return(2);
416 : : } /* end if ( rep == 0 ) */
417 : :
418 : 30 : b[0] /= diag[0]; /* Vorwaertselemination */
419 [ + + ]: 360 : for (i = 1; i < n-1; i++)
420 : 330 : b[i] = ( b[i] - b[i-1] * lower[i] ) / diag[i];
421 : :
422 [ + + ]: 360 : for (temp = 0.0, i = 0; i < n-2; i++)
423 : 330 : temp -= lowrow[i] * b[i];
424 : :
425 : 30 : b[n-1] = ( b[n-1] + temp - lower[n-1] * b[n-2] ) / diag[n-1];
426 : :
427 : 30 : b[n-2] -= b[n-1] * upper[n-2]; /* Rueckwaertselimination */
428 [ + + ]: 360 : for (j = n-3; j >= 0; j--) {
429 : 330 : i=j;
430 : 330 : b[i] -= upper[i] * b[i+1] + ricol[i] * b[n-1];
431 : : }
432 : 30 : return(0);
433 : : }
434 : :
435 : : /*------------------ ENDE ZYKLISCH TRIDIAGONAL ---------------------*/
436 : :
437 : :
438 : : } // extern "C"
439 : :
440 : :
441 : : /*************************************************************************
442 : : |*
443 : : |* NaturalSpline()
444 : : |*
445 : : |* Beschreibung Berechnet die Koeffizienten eines natuerlichen
446 : : |* kubischen Polynomsplines mit n Stuetzstellen.
447 : : |*
448 : : *************************************************************************/
449 : :
450 : 0 : sal_uInt16 NaturalSpline(sal_uInt16 n, double* x, double* y,
451 : : double Marg0, double MargN,
452 : : sal_uInt8 MargCond,
453 : : double* b, double* c, double* d)
454 : : {
455 : : sal_uInt16 i;
456 : : double* a;
457 : : double* h;
458 : : sal_uInt16 error;
459 : :
460 [ # # ]: 0 : if (n<2) return 1;
461 [ # # ]: 0 : if ( (MargCond & ~3) ) return 2;
462 : 0 : a=new double[n+1];
463 : 0 : h=new double[n+1];
464 [ # # ]: 0 : for (i=0;i<n;i++) {
465 : 0 : h[i]=x[i+1]-x[i];
466 [ # # ][ # # ]: 0 : if (h[i]<=0.0) { delete[] a; delete[] h; return 1; }
[ # # ]
467 : : }
468 [ # # ]: 0 : for (i=0;i<n-1;i++) {
469 : 0 : a[i]=3.0*((y[i+2]-y[i+1])/h[i+1]-(y[i+1]-y[i])/h[i]);
470 : 0 : b[i]=h[i];
471 : 0 : c[i]=h[i+1];
472 : 0 : d[i]=2.0*(h[i]+h[i+1]);
473 : : }
474 [ # # # # : 0 : switch (MargCond) {
# ]
475 : : case 0: {
476 [ # # ]: 0 : if (n==2) {
477 : 0 : a[0]=a[0]/3.0;
478 : 0 : d[0]=d[0]*0.5;
479 : : } else {
480 : 0 : a[0] =a[0]*h[1]/(h[0]+h[1]);
481 : 0 : a[n-2]=a[n-2]*h[n-2]/(h[n-1]+h[n-2]);
482 : 0 : d[0] =d[0]-h[0];
483 : 0 : d[n-2]=d[n-2]-h[n-1];
484 : 0 : c[0] =c[0]-h[0];
485 : 0 : b[n-2]=b[n-2]-h[n-1];
486 : : }
487 : : }
488 : : case 1: {
489 : 0 : a[0] =a[0]-1.5*((y[1]-y[0])/h[0]-Marg0);
490 : 0 : a[n-2]=a[n-2]-1.5*(MargN-(y[n]-y[n-1])/h[n-1]);
491 : 0 : d[0] =d[0]-h[0]*0.5;
492 : 0 : d[n-2]=d[n-2]-h[n-1]*0.5;
493 : : }
494 : : case 2: {
495 : 0 : a[0] =a[0]-h[0]*Marg0*0.5;
496 : 0 : a[n-2]=a[n-2]-h[n-1]*MargN*0.5;
497 : : }
498 : : case 3: {
499 : 0 : a[0] =a[0]+Marg0*h[0]*h[0]*0.5;
500 : 0 : a[n-2]=a[n-2]-MargN*h[n-1]*h[n-1]*0.5;
501 : 0 : d[0] =d[0]+h[0];
502 : 0 : d[n-2]=d[n-2]+h[n-1];
503 : : }
504 : : } // switch MargCond
505 [ # # ]: 0 : if (n==2) {
506 : 0 : c[1]=a[0]/d[0];
507 : : } else {
508 : 0 : error=TriDiagGS(sal_False,n-1,b,d,c,a);
509 [ # # ][ # # ]: 0 : if (error!=0) { delete[] a; delete[] h; return error+2; }
[ # # ]
510 [ # # ]: 0 : for (i=0;i<n-1;i++) c[i+1]=a[i];
511 : : }
512 [ # # # # : 0 : switch (MargCond) {
# ]
513 : : case 0: {
514 [ # # ]: 0 : if (n==2) {
515 : 0 : c[2]=c[1];
516 : 0 : c[0]=c[1];
517 : : } else {
518 : 0 : c[0]=c[1]+h[0]*(c[1]-c[2])/h[1];
519 : 0 : c[n]=c[n-1]+h[n-1]*(c[n-1]-c[n-2])/h[n-2];
520 : : }
521 : : }
522 : : case 1: {
523 : 0 : c[0]=1.5*((y[1]-y[0])/h[0]-Marg0);
524 : 0 : c[0]=(c[0]-c[1]*h[0]*0.5)/h[0];
525 : 0 : c[n]=1.5*((y[n]-y[n-1])/h[n-1]-MargN);
526 : 0 : c[n]=(c[n]-c[n-1]*h[n-1]*0.5)/h[n-1];
527 : : }
528 : : case 2: {
529 : 0 : c[0]=Marg0*0.5;
530 : 0 : c[n]=MargN*0.5;
531 : : }
532 : : case 3: {
533 : 0 : c[0]=c[1]-Marg0*h[0]*0.5;
534 : 0 : c[n]=c[n-1]+MargN*h[n-1]*0.5;
535 : : }
536 : : } // switch MargCond
537 [ # # ]: 0 : for (i=0;i<n;i++) {
538 : 0 : b[i]=(y[i+1]-y[i])/h[i]-h[i]*(c[i+1]+2.0*c[i])/3.0;
539 : 0 : d[i]=(c[i+1]-c[i])/(3.0*h[i]);
540 : : }
541 [ # # ]: 0 : delete[] a;
542 [ # # ]: 0 : delete[] h;
543 : 0 : return 0;
544 : : }
545 : :
546 : :
547 : : /*************************************************************************
548 : : |*
549 : : |* PeriodicSpline()
550 : : |*
551 : : |* Beschreibung Berechnet die Koeffizienten eines periodischen
552 : : |* kubischen Polynomsplines mit n Stuetzstellen.
553 : : |*
554 : : *************************************************************************/
555 : :
556 : :
557 : 30 : sal_uInt16 PeriodicSpline(sal_uInt16 n, double* x, double* y,
558 : : double* b, double* c, double* d)
559 : : { // Arrays muessen von [0..n] dimensioniert sein!
560 : : sal_uInt16 Error;
561 : : sal_uInt16 i,im1,nm1; //integer
562 : : double hr,hl;
563 : : double* a;
564 : : double* lowrow;
565 : : double* ricol;
566 : :
567 [ - + ]: 30 : if (n<2) return 4;
568 : 30 : nm1=n-1;
569 [ - + ][ + + ]: 420 : for (i=0;i<=nm1;i++) if (x[i+1]<=x[i]) return 2; // muss streng nonoton fallend sein!
570 [ - + ]: 30 : if (y[n]!=y[0]) return 3; // Anfang muss gleich Ende sein!
571 : :
572 : 30 : a =new double[n+1];
573 : 30 : lowrow=new double[n+1];
574 : 30 : ricol =new double[n+1];
575 : :
576 [ - + ]: 30 : if (n==2) {
577 : 0 : c[1]=3.0*((y[2]-y[1])/(x[2]-x[1]));
578 : 0 : c[1]=c[1]-3.0*((y[i]-y[0])/(x[1]-x[0]));
579 : 0 : c[1]=c[1]/(x[2]-x[0]);
580 : 0 : c[2]=-c[1];
581 : : } else {
582 [ + + ]: 390 : for (i=1;i<=nm1;i++) {
583 : 360 : im1=i-1;
584 : 360 : hl=x[i]-x[im1];
585 : 360 : hr=x[i+1]-x[i];
586 : 360 : b[im1]=hl;
587 : 360 : d[im1]=2.0*(hl+hr);
588 : 360 : c[im1]=hr;
589 : 360 : a[im1]=3.0*((y[i+1]-y[i])/hr-(y[i]-y[im1])/hl);
590 : : }
591 : 30 : hl=x[n]-x[nm1];
592 : 30 : hr=x[1]-x[0];
593 : 30 : b[nm1]=hl;
594 : 30 : d[nm1]=2.0*(hl+hr);
595 : 30 : lowrow[0]=hr;
596 : 30 : ricol[0]=hr;
597 : 30 : a[nm1]=3.0*((y[1]-y[0])/hr-(y[n]-y[nm1])/hl);
598 : 30 : Error=ZyklTriDiagGS(sal_False,n,b,d,c,lowrow,ricol,a);
599 [ - + ]: 30 : if ( Error != 0 )
600 : : {
601 [ # # ]: 0 : delete[] a;
602 [ # # ]: 0 : delete[] lowrow;
603 [ # # ]: 0 : delete[] ricol;
604 : 0 : return(Error+4);
605 : : }
606 [ + + ]: 420 : for (i=0;i<=nm1;i++) c[i+1]=a[i];
607 : : }
608 : 30 : c[0]=c[n];
609 [ + + ]: 420 : for (i=0;i<=nm1;i++) {
610 : 390 : hl=x[i+1]-x[i];
611 : 390 : b[i]=(y[i+1]-y[i])/hl;
612 : 390 : b[i]=b[i]-hl*(c[i+1]+2.0*c[i])/3.0;
613 : 390 : d[i]=(c[i+1]-c[i])/hl/3.0;
614 : : }
615 [ + - ]: 30 : delete[] a;
616 [ + - ]: 30 : delete[] lowrow;
617 [ + - ]: 30 : delete[] ricol;
618 : 30 : return 0;
619 : : }
620 : :
621 : :
622 : :
623 : : /*************************************************************************
624 : : |*
625 : : |* ParaSpline()
626 : : |*
627 : : |* Beschreibung Berechnet die Koeffizienten eines parametrischen
628 : : |* natuerlichen oder periodischen kubischen
629 : : |* Polynomsplines mit n Stuetzstellen.
630 : : |*
631 : : *************************************************************************/
632 : :
633 : 15 : sal_uInt16 ParaSpline(sal_uInt16 n, double* x, double* y, sal_uInt8 MargCond,
634 : : double Marg01, double Marg02,
635 : : double MargN1, double MargN2,
636 : : sal_Bool CondT, double* T,
637 : : double* bx, double* cx, double* dx,
638 : : double* by, double* cy, double* dy)
639 : : {
640 : : sal_uInt16 Error;
641 : : sal_uInt16 i;
642 : : double deltX,deltY,delt,
643 : 15 : alphX = 0,alphY = 0,
644 : 15 : betX = 0,betY = 0;
645 : :
646 [ - + ]: 15 : if (n<2) return 1;
647 [ - + ][ # # ]: 15 : if ((MargCond & ~3) && (MargCond != 4)) return 2; // ungueltige Randbedingung
648 [ + - ]: 15 : if (CondT==sal_False) {
649 : 15 : T[0]=0.0;
650 [ + + ]: 210 : for (i=0;i<n;i++) {
651 : 195 : deltX=x[i+1]-x[i]; deltY=y[i+1]-y[i];
652 : 195 : delt =deltX*deltX+deltY*deltY;
653 [ - + ]: 195 : if (delt<=0.0) return 3; // zwei identische Punkte nacheinander!
654 : 195 : T[i+1]=T[i]+sqrt(delt);
655 : : }
656 : : }
657 [ - - + - : 15 : switch (MargCond) {
- ]
658 : 0 : case 0: break;
659 : : case 1: case 2: {
660 : 0 : alphX=Marg01; betX=MargN1;
661 : 0 : alphY=Marg02; betY=MargN2;
662 : 0 : } break;
663 : : case 3: {
664 [ - + ]: 15 : if (x[n]!=x[0]) return 3;
665 [ - + ]: 15 : if (y[n]!=y[0]) return 4;
666 : 15 : } break;
667 : : case 4: {
668 [ # # ][ # # ]: 0 : if (abs(Marg01)>=MAXROOT) {
669 : 0 : alphX=0.0;
670 [ # # ]: 0 : alphY=sign(1.0,y[1]-y[0]);
671 : : } else {
672 [ # # ][ # # ]: 0 : alphX=sign(sqrt(1.0/(1.0+Marg01*Marg01)),x[1]-x[0]);
[ # # ]
673 : 0 : alphY=alphX*Marg01;
674 : : }
675 [ # # ][ # # ]: 0 : if (abs(MargN1)>=MAXROOT) {
676 : 0 : betX=0.0;
677 [ # # ]: 0 : betY=sign(1.0,y[n]-y[n-1]);
678 : : } else {
679 [ # # ][ # # ]: 0 : betX=sign(sqrt(1.0/(1.0+MargN1*MargN1)),x[n]-x[n-1]);
[ # # ]
680 : 0 : betY=betX*MargN1;
681 : : }
682 : : }
683 : : } // switch MargCond
684 [ + - ]: 15 : if (MargCond==3) {
685 : 15 : Error=PeriodicSpline(n,T,x,bx,cx,dx);
686 [ - + ]: 15 : if (Error!=0) return(Error+4);
687 : 15 : Error=PeriodicSpline(n,T,y,by,cy,dy);
688 [ - + ]: 15 : if (Error!=0) return(Error+10);
689 : : } else {
690 : 0 : Error=NaturalSpline(n,T,x,alphX,betX,MargCond,bx,cx,dx);
691 [ # # ]: 0 : if (Error!=0) return(Error+4);
692 : 0 : Error=NaturalSpline(n,T,y,alphY,betY,MargCond,by,cy,dy);
693 [ # # ]: 0 : if (Error!=0) return(Error+9);
694 : : }
695 : 15 : return 0;
696 : : }
697 : :
698 : :
699 : :
700 : : /*************************************************************************
701 : : |*
702 : : |* CalcSpline()
703 : : |*
704 : : |* Beschreibung Berechnet die Koeffizienten eines parametrischen
705 : : |* natuerlichen oder periodischen kubischen
706 : : |* Polynomsplines. Die Eckpunkte des uebergebenen
707 : : |* Polygons werden als Stuetzstellen angenommen.
708 : : |* n liefert die Anzahl der Teilpolynome.
709 : : |* Ist die Berechnung fehlerfrei verlaufen, so
710 : : |* liefert die Funktion sal_True. Nur in diesem Fall
711 : : |* ist Speicher fuer die Koeffizientenarrays
712 : : |* allokiert, der dann spaeter vom Aufrufer mittels
713 : : |* delete freizugeben ist.
714 : : |*
715 : : *************************************************************************/
716 : :
717 : 15 : sal_Bool CalcSpline(Polygon& rPoly, sal_Bool Periodic, sal_uInt16& n,
718 : : double*& ax, double*& ay, double*& bx, double*& by,
719 : : double*& cx, double*& cy, double*& dx, double*& dy, double*& T)
720 : : {
721 : : sal_uInt8 Marg;
722 : : double Marg01;
723 : : double MargN1,MargN2;
724 : : sal_uInt16 i;
725 : 15 : Point P0(-32768,-32768);
726 : 15 : Point Pt;
727 : :
728 [ + - ]: 15 : n=rPoly.GetSize();
729 [ + - ][ + - ]: 15 : ax=new double[rPoly.GetSize()+2];
730 [ + - ][ + - ]: 15 : ay=new double[rPoly.GetSize()+2];
731 : :
732 : 15 : n=0;
733 [ + - ][ + + ]: 210 : for (i=0;i<rPoly.GetSize();i++) {
734 [ + - ]: 195 : Pt=rPoly.GetPoint(i);
735 [ + + ][ + - ]: 195 : if (i==0 || Pt!=P0) {
[ + - ]
736 : 195 : ax[n]=Pt.X();
737 : 195 : ay[n]=Pt.Y();
738 : 195 : n++;
739 : 195 : P0=Pt;
740 : : }
741 : : }
742 : :
743 [ + - ]: 15 : if (Periodic) {
744 : 15 : Marg=3;
745 : 15 : ax[n]=ax[0];
746 : 15 : ay[n]=ay[0];
747 : 15 : n++;
748 : : } else {
749 : 0 : Marg=2;
750 : : }
751 : :
752 [ + - ]: 15 : bx=new double[n+1];
753 [ + - ]: 15 : by=new double[n+1];
754 [ + - ]: 15 : cx=new double[n+1];
755 [ + - ]: 15 : cy=new double[n+1];
756 [ + - ]: 15 : dx=new double[n+1];
757 [ + - ]: 15 : dy=new double[n+1];
758 [ + - ]: 15 : T =new double[n+1];
759 : :
760 : 15 : Marg01=0.0;
761 : 15 : MargN1=0.0;
762 : 15 : MargN2=0.0;
763 [ + - ]: 15 : if (n>0) n--; // n Korregieren (Anzahl der Teilpolynome)
764 : :
765 : 15 : sal_Bool bRet = sal_False;
766 [ + - ][ - + ]: 15 : if ( ( Marg == 3 && n >= 3 ) || ( Marg == 2 && n >= 2 ) )
[ # # ][ # # ]
767 : : {
768 [ + - ]: 15 : bRet = ParaSpline(n,ax,ay,Marg,Marg01,Marg01,MargN1,MargN2,sal_False,T,bx,cx,dx,by,cy,dy) == 0;
769 : : }
770 [ - + ]: 15 : if ( bRet == sal_False )
771 : : {
772 [ # # ]: 0 : delete[] ax;
773 [ # # ]: 0 : delete[] ay;
774 [ # # ]: 0 : delete[] bx;
775 [ # # ]: 0 : delete[] by;
776 [ # # ]: 0 : delete[] cx;
777 [ # # ]: 0 : delete[] cy;
778 [ # # ]: 0 : delete[] dx;
779 [ # # ]: 0 : delete[] dy;
780 [ # # ]: 0 : delete[] T;
781 : 0 : n=0;
782 : : }
783 : 15 : return bRet;
784 : : }
785 : :
786 : :
787 : : /*************************************************************************
788 : : |*
789 : : |* Spline2Poly()
790 : : |*
791 : : |* Beschreibung Konvertiert einen parametrichen kubischen
792 : : |* Polynomspline Spline (natuerlich oder periodisch)
793 : : |* in ein angenaehertes Polygon.
794 : : |* Die Funktion liefert sal_False, wenn ein Fehler bei
795 : : |* der Koeffizientenberechnung aufgetreten ist oder
796 : : |* das Polygon zu gross wird (>PolyMax=16380). Im 1.
797 : : |* Fall hat das Polygon 0, im 2. Fall PolyMax Punkte.
798 : : |* Um Koordinatenueberlaeufe zu vermeiden werden diese
799 : : |* auf +/-32000 begrenzt.
800 : : |*
801 : : *************************************************************************/
802 : 15 : sal_Bool Spline2Poly(Polygon& rSpln, sal_Bool Periodic, Polygon& rPoly)
803 : : {
804 : 15 : short MinKoord=-32000; // zur Vermeidung
805 : 15 : short MaxKoord=32000; // von Ueberlaeufen
806 : :
807 : : double* ax; // Koeffizienten der Polynome
808 : : double* ay;
809 : : double* bx;
810 : : double* by;
811 : : double* cx;
812 : : double* cy;
813 : : double* dx;
814 : : double* dy;
815 : : double* tv;
816 : :
817 : : double Step; // Schrittweite fuer t
818 : : double dt1,dt2,dt3; // Delta t, y, ^3
819 : : double t;
820 : : sal_Bool bEnde; // Teilpolynom zu Ende?
821 : : sal_uInt16 n; // Anzahl der zu zeichnenden Teilpolynome
822 : : sal_uInt16 i; // aktuelles Teilpolynom
823 : : sal_Bool bOk; // noch alles ok?
824 : 15 : sal_uInt16 PolyMax=16380;// Maximale Anzahl von Polygonpunkten
825 : : long x,y;
826 : :
827 [ + - ]: 15 : bOk=CalcSpline(rSpln,Periodic,n,ax,ay,bx,by,cx,cy,dx,dy,tv);
828 [ + - ]: 15 : if (bOk) {
829 : 15 : Step =10;
830 : :
831 [ + - ]: 15 : rPoly.SetSize(1);
832 [ + - ]: 15 : rPoly.SetPoint(Point(short(ax[0]),short(ay[0])),0); // erster Punkt
833 : 15 : i=0;
834 [ + + ]: 210 : while (i<n) { // n Teilpolynome malen
835 : 195 : t=tv[i]+Step;
836 : 195 : bEnde=sal_False;
837 [ + + ]: 53375 : while (!bEnde) { // ein Teilpolynom interpolieren
838 : 53180 : bEnde=t>=tv[i+1];
839 [ + + ]: 53180 : if (bEnde) t=tv[i+1];
840 : 53180 : dt1=t-tv[i]; dt2=dt1*dt1; dt3=dt2*dt1;
841 : 53180 : x=long(ax[i]+bx[i]*dt1+cx[i]*dt2+dx[i]*dt3);
842 : 53180 : y=long(ay[i]+by[i]*dt1+cy[i]*dt2+dy[i]*dt3);
843 [ - + ][ - + ]: 53180 : if (x<MinKoord) x=MinKoord; if (x>MaxKoord) x=MaxKoord;
844 [ - + ][ - + ]: 53180 : if (y<MinKoord) y=MinKoord; if (y>MaxKoord) y=MaxKoord;
845 [ + - ][ + - ]: 53180 : if (rPoly.GetSize()<PolyMax) {
846 [ + - ][ + - ]: 53180 : rPoly.SetSize(rPoly.GetSize()+1);
847 [ + - ][ + - ]: 53180 : rPoly.SetPoint(Point(short(x),short(y)),rPoly.GetSize()-1);
848 : : } else {
849 : 0 : bOk=sal_False; // Fehler: Polygon wird zu gross
850 : : }
851 : 53180 : t=t+Step;
852 : : } // Ende von Teilpolynom
853 : 195 : i++; // naechstes Teilpolynom
854 : : }
855 [ + - ]: 15 : delete[] ax;
856 [ + - ]: 15 : delete[] ay;
857 [ + - ]: 15 : delete[] bx;
858 [ + - ]: 15 : delete[] by;
859 [ + - ]: 15 : delete[] cx;
860 [ + - ]: 15 : delete[] cy;
861 [ + - ]: 15 : delete[] dx;
862 [ + - ]: 15 : delete[] dy;
863 [ + - ]: 15 : delete[] tv;
864 : 15 : return bOk;
865 : : } // Ende von if (bOk)
866 [ # # ]: 0 : rPoly.SetSize(0);
867 : 15 : return sal_False;
868 : : }
869 : :
870 : : /* vim:set shiftwidth=4 softtabstop=4 expandtab: */
|
