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