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