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