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