LORENE
op_scost.C
1/*
2 * Copyright (c) 1999-2001 Jerome Novak
3 *
4 * This file is part of LORENE.
5 *
6 * LORENE is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * LORENE is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with LORENE; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 *
20 */
21
22
23
24
25/*
26 * $Id: op_scost.C,v 1.8 2016/12/05 16:18:08 j_novak Exp $
27 * $Log: op_scost.C,v $
28 * Revision 1.8 2016/12/05 16:18:08 j_novak
29 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
30 *
31 * Revision 1.7 2014/10/13 08:53:26 j_novak
32 * Lorene classes and functions now belong to the namespace Lorene.
33 *
34 * Revision 1.6 2009/10/10 18:28:11 j_novak
35 * New bases T_COS and T_SIN.
36 *
37 * Revision 1.5 2007/12/21 10:43:38 j_novak
38 * Corrected some bugs in the case nt=1 (1 point in theta).
39 *
40 * Revision 1.4 2007/10/05 12:37:20 j_novak
41 * Corrected a few errors in the theta-nonsymmetric case (bases T_COSSIN_C and
42 * T_COSSIN_S).
43 *
44 * Revision 1.3 2005/02/16 15:29:40 m_forot
45 * Correct T_COSSIN_S and T_COSSIN_C cases
46 *
47 * Revision 1.2 2004/11/23 15:16:01 m_forot
48 *
49 * Added the bases for the cases without any equatorial symmetry
50 * (T_COSSIN_C, T_COSSIN_S, T_LEG, R_CHEBPI_P, R_CHEBPI_I).
51 *
52 * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
53 * LORENE
54 *
55 * Revision 1.3 2000/02/24 16:42:36 eric
56 * Initialisation a zero du tableau xo avant le calcul.
57 *
58 * Revision 1.2 1999/11/22 14:34:16 novak
59 * Suppression des variables inutiles
60 *
61 * Revision 1.1 1999/11/16 13:38:00 novak
62 * Initial revision
63 *
64 *
65 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_scost.C,v 1.8 2016/12/05 16:18:08 j_novak Exp $
66 *
67 */
68
69/*
70 * Ensemble des routines de base de division par rapport a cost theta
71 * (Utilisation interne)
72 *
73 * void _scost_XXXX(Tbl * t, int & b)
74 * t pointeur sur le Tbl d'entree-sortie
75 * b base spectrale
76 *
77 */
78
79
80// Fichier includes
81#include "tbl.h"
82
83
84 //-----------------------------------
85 // Routine pour les cas non prevus --
86 //-----------------------------------
87
88namespace Lorene {
89void _scost_pas_prevu(Tbl * tb, int& base) {
90 cout << "scost pas prevu..." << endl ;
91 cout << "Tbl: " << tb << " base: " << base << endl ;
92 abort () ;
93 exit(-1) ;
94}
95
96 //--------------
97 // cas T_COS
98 //--------------
99
100void _scost_t_cos(Tbl* tb, int & b)
101{
102
103 // Peut-etre rien a faire ?
104 if (tb->get_etat() == ETATZERO) {
105 int base_r = b & MSQ_R ;
106 int base_p = b & MSQ_P ;
107 switch(base_r){
108 case(R_CHEBPI_P):
109 b = R_CHEBPI_I | base_p | T_COS ;
110 break ;
111 case(R_CHEBPI_I):
112 b = R_CHEBPI_P | base_p | T_COS ;
113 break ;
114 default:
115 b = base_r | base_p | T_COS ;
116 break;
117 }
118 return ;
119 }
120
121 // Protection
122 assert(tb->get_etat() == ETATQCQ) ;
123
124 // Pour le confort : nbre de points reels.
125 int nr = (tb->dim).dim[0] ;
126 int nt = (tb->dim).dim[1] ;
127 int np = (tb->dim).dim[2] ;
128 np = np - 2 ;
129
130 // pt. sur le tableau de double resultat
131 double* somP = new double [nr] ;
132 double* somN = new double [nr] ;
133
134 // pt. sur le tableau de double resultat
135 double* xo = new double[(tb->dim).taille] ;
136
137 // Initialisation a zero :
138 for (int i=0; i<(tb->dim).taille; i++) {
139 xo[i] = 0 ;
140 }
141
142 // On y va...
143 double* xi = tb->t ;
144 double* xci = xi ; // Pointeurs
145 double* xco = xo ; // courants
146
147 // k = 0
148
149 // Dernier j: j = nt-1
150 // Positionnement
151 xci += nr * (nt-1) ;
152 xco += nr * (nt-1) ;
153
154 // Initialisation de som
155 for (int i=0 ; i<nr ; i++) {
156 somP[i] = 0. ;
157 somN[i] = 0. ;
158 xco[i] = 0. ; // mise a zero du dernier coef en theta
159 }
160
161 // j suivants
162 for (int j=nt-2 ; j >= 0 ; j--) {
163 // Positionnement
164 xci -= nr ;
165 xco -= nr ;
166 // Calcul
167 for (int i=0 ; i<nr ; i++ ) {
168 if((j%2) == 1) {
169 somN[i] = -somN[i] ;
170 somN[i] += 2*xci[i] ;
171 xco[i] = somP[i] ;
172 }
173 else {
174 somP[i] = -somP[i] ;
175 somP[i] += 2*xci[i] ;
176 xco[i] = somN[i] ;
177 }
178 } // Fin de la boucle sur r
179 } // Fin de la boucle sur theta
180 // j=0 : le facteur 2...
181 for (int i=0 ; i<nr ; i++) xco[i] *= .5 ;
182
183 // Positionnement phi suivant
184 xci += nr*nt ;
185 xco += nr*nt ;
186
187 // k = 1
188 xci += nr*nt ;
189 xco += nr*nt ;
190
191 // k >= 2
192 for (int k=2 ; k<np+1 ; k++) {
193
194 // Dernier j: j = nt-1
195 // Positionnement
196 xci += nr * (nt-1) ;
197 xco += nr * (nt-1) ;
198
199 // Initialisation de som
200 for (int i=0 ; i<nr ; i++) {
201 somP[i] = 0. ;
202 somN[i] = 0. ;
203 xco[i] = 0. ; // mise a zero dui dernier coefficient en theta.
204 }
205
206 // j suivants
207 for (int j=nt-2 ; j >= 0 ; j--) {
208 // Positionnement
209 xci -= nr ;
210 xco -= nr ;
211 // Calcul
212 for (int i=0 ; i<nr ; i++ ) {
213 if((j%2) == 1) {
214 somN[i] = -somN[i];
215 somN[i] += 2 * xci[i] ;
216 xco[i] = somP[i] ;
217 }
218 else {
219 somP[i] = -somP[i];
220 somP[i] += 2 * xci[i] ;
221 xco[i] = somN[i];
222 }
223 } // Fin de la boucle sur r
224 } // Fin de la boucle sur theta
225 // j=0 : facteur 2 ...
226 for (int i=0 ; i<nr ; i++) xco[i] *= 0.5 ;
227
228 // Positionnement phi suivant
229 xci += nr*nt ;
230 xco += nr*nt ;
231 } // Fin de la boucle sur phi
232
233 // On remet les choses la ou il faut
234 delete [] tb->t ;
235 tb->t = xo ;
236
237 //Menage
238 delete [] somP ;
239 delete [] somN ;
240
241 // base de developpement
242 int base_r = b & MSQ_R ;
243 int base_p = b & MSQ_P ;
244 switch(base_r){
245 case(R_CHEBPI_P):
246 b = R_CHEBPI_I | base_p | T_COS ;
247 break ;
248 case(R_CHEBPI_I):
249 b = R_CHEBPI_P | base_p | T_COS ;
250 break ;
251 default:
252 b = base_r | base_p | T_COS ;
253 break;
254 }
255}
256
257 //------------
258 // cas T_SIN
259 //------------
260
261void _scost_t_sin(Tbl* tb, int & b)
262{
263
264
265 // Peut-etre rien a faire ?
266 if (tb->get_etat() == ETATZERO) {
267 int base_r = b & MSQ_R ;
268 int base_p = b & MSQ_P ;
269 switch(base_r){
270 case(R_CHEBPI_P):
271 b = R_CHEBPI_I | base_p | T_SIN ;
272 break ;
273 case(R_CHEBPI_I):
274 b = R_CHEBPI_P | base_p | T_SIN ;
275 break ;
276 default:
277 b = base_r | base_p | T_SIN ;
278 break;
279 }
280 return ;
281 }
282
283 // Protection
284 assert(tb->get_etat() == ETATQCQ) ;
285
286 // Pour le confort : nbre de points reels.
287 int nr = (tb->dim).dim[0] ;
288 int nt = (tb->dim).dim[1] ;
289 int np = (tb->dim).dim[2] ;
290 np = np - 2 ;
291
292 // zone de travail
293 double* somP = new double [nr] ;
294 double* somN = new double [nr] ;
295
296 // pt. sur le tableau de double resultat
297 double* xo = new double[(tb->dim).taille] ;
298
299 // Initialisation a zero :
300 for (int i=0; i<(tb->dim).taille; i++) {
301 xo[i] = 0 ;
302 }
303
304 // On y va...
305 double* xi = tb->t ;
306 double* xci = xi ; // Pointeurs
307 double* xco = xo ; // courants
308
309 // k = 0
310
311 // Dernier j: j = nt-1
312 // Positionnement
313 xci += nr * (nt-1) ;
314 xco += nr * (nt-1) ;
315
316 // Initialisation de som
317 for (int i=0 ; i<nr ; i++) {
318 somP[i] = 0. ;
319 somN[i] = 0. ;
320 xco[i] = 0. ;
321 }
322
323 // j suivants
324 for (int j=nt-2 ; j >= 0 ; j--) {
325 // Positionnement
326 xci -= nr ;
327 xco -= nr ;
328 // Calcul
329 for (int i=0 ; i<nr ; i++ ) {
330 if((j%2) == 1) {
331 somN[i] = -somN[i] ;
332 somN[i] += 2*xci[i] ;
333 xco[i] = somP[i] ;
334 }
335 else {
336 somP[i] = -somP[i] ;
337 somP[i] += 2*xci[i] ;
338 xco[i] = somN[i] ;
339 }
340 } // Fin de la boucle sur r
341 } // Fin de la boucle sur theta
342 // j=0 : sin(0*theta)
343 for (int i=0 ; i<nr ; i++) xco[i] = 0. ;
344
345 // Positionnement phi suivant
346 xci += nr*nt ;
347 xco += nr*nt ;
348
349 // k = 1
350 xci += nr*nt ;
351 xco += nr*nt ;
352
353 // k >= 2
354 for (int k=2 ; k<np+1 ; k++) {
355
356 // Dernier j: j = nt-1
357 // Positionnement
358 xci += nr * (nt-1) ;
359 xco += nr * (nt-1) ;
360
361 // Initialisation de som
362 for (int i=0 ; i<nr ; i++) {
363 somP[i] = 0. ;
364 somN[i] = 0. ;
365 xco[i] = 0. ;
366 }
367
368 // j suivants
369 for (int j=nt-2 ; j >= 0 ; j--) {
370 // Positionnement
371 xci -= nr ;
372 xco -= nr ;
373 // Calcul
374 for (int i=0 ; i<nr ; i++ ) {
375 if((j%2) == 1) {
376 somN[i] = -somN[i] ;
377 somN[i] += 2*xci[i] ;
378 xco[i] = somP[i] ;
379 }
380 else {
381 somP[i] = -somP[i] ;
382 somP[i] += 2*xci[i] ;
383 xco[i] = somN[i] ;
384 }
385 } // Fin de la boucle sur r
386 } // Fin de la boucle sur theta
387 // j=0 : sin(0*theta)
388 for (int i=0 ; i<nr ; i++) xco[i] = 0. ;
389
390 // Positionnement phi suivant
391 xci += nr*nt ;
392 xco += nr*nt ;
393
394 } // Fin de la boucle sur phi
395
396 // On remet les choses la ou il faut
397 delete [] tb->t ;
398 tb->t = xo ;
399
400 //Menage
401 delete [] somP ;
402 delete [] somN ;
403
404 // base de developpement
405 int base_r = b & MSQ_R ;
406 int base_p = b & MSQ_P ;
407 switch(base_r){
408 case(R_CHEBPI_P):
409 b = R_CHEBPI_I | base_p | T_SIN ;
410 break ;
411 case(R_CHEBPI_I):
412 b = R_CHEBPI_P | base_p | T_SIN ;
413 break ;
414 default:
415 b = base_r | base_p | T_SIN ;
416 break;
417 }
418}
419 //----------------
420 // cas T_COS_P
421 //----------------
422
423void _scost_t_cos_p(Tbl* tb, int & b)
424{
425
426 // Peut-etre rien a faire ?
427 if (tb->get_etat() == ETATZERO) {
428 int base_r = b & MSQ_R ;
429 int base_p = b & MSQ_P ;
430 b = base_r | base_p | T_COS_I ;
431 return ;
432 }
433
434 // Protection
435 assert(tb->get_etat() == ETATQCQ) ;
436
437 // Pour le confort : nbre de points reels.
438 int nr = (tb->dim).dim[0] ;
439 int nt = (tb->dim).dim[1] ;
440 int np = (tb->dim).dim[2] ;
441 np = np - 2 ;
442
443 // zone de travail
444 double* som = new double [nr] ;
445
446 // pt. sur le tableau de double resultat
447 double* xo = new double[(tb->dim).taille] ;
448
449 // Initialisation a zero :
450 for (int i=0; i<(tb->dim).taille; i++) {
451 xo[i] = 0 ;
452 }
453
454 // On y va...
455 double* xi = tb->t ;
456 double* xci = xi ; // Pointeurs
457 double* xco = xo ; // courants
458
459 // k = 0
460
461 // Dernier j: j = nt-1
462 // Positionnement
463 xci += nr * (nt) ;
464 xco += nr * (nt-1) ;
465
466 // Initialisation de som
467 for (int i=0 ; i<nr ; i++) {
468 som[i] = 0. ;
469 xco[i] = 0. ; // mise a zero du dernier coef en theta
470 }
471
472 // j suivants
473 for (int j=nt-2 ; j >= 0 ; j--) {
474 // Positionnement
475 xci -= nr ;
476 xco -= nr ;
477 // Calcul
478 for (int i=0 ; i<nr ; i++ ) {
479 som[i] += 2*xci[i] ;
480 xco[i] = som[i] ;
481 som[i] = -som[i] ;
482 } // Fin de la boucle sur r
483 } // Fin de la boucle sur theta
484 // Positionnement phi suivant
485 xci -= nr ;
486 xci += nr*nt ;
487 xco += nr*nt ;
488
489 // k = 1
490 xci += nr*nt ;
491 xco += nr*nt ;
492
493 // k >= 2
494 for (int k=2 ; k<np+1 ; k++) {
495
496 // Dernier j: j = nt-1
497 // Positionnement
498 xci += nr * (nt) ;
499 xco += nr * (nt-1) ;
500
501 // Initialisation de som
502 for (int i=0 ; i<nr ; i++) {
503 som[i] = 0. ;
504 xco[i] = 0. ; // mise a zero dui dernier coefficient en theta.
505 }
506
507 // j suivants
508 for (int j=nt-2 ; j >= 0 ; j--) {
509 // Positionnement
510 xci -= nr ;
511 xco -= nr ;
512 // Calcul
513 for (int i=0 ; i<nr ; i++ ) {
514 som[i] += 2*xci[i] ;
515 xco[i] = som[i] ;
516 som[i] = - som[i] ;
517 } // Fin de la boucle sur r
518 } // Fin de la boucle sur theta
519 // Positionnement phi suivant
520 xci -= nr ;
521 xci += nr*nt ;
522 xco += nr*nt ;
523 } // Fin de la boucle sur phi
524
525 // On remet les choses la ou il faut
526 delete [] tb->t ;
527 tb->t = xo ;
528
529 //Menage
530 delete [] som ;
531
532 // base de developpement
533 int base_r = b & MSQ_R ;
534 int base_p = b & MSQ_P ;
535 b = base_r | base_p | T_COS_I ;
536}
537
538 //--------------
539 // cas T_SIN_P
540 //--------------
541
542void _scost_t_sin_p(Tbl* tb, int & b)
543{
544
545
546 // Peut-etre rien a faire ?
547 if (tb->get_etat() == ETATZERO) {
548 int base_r = b & MSQ_R ;
549 int base_p = b & MSQ_P ;
550 b = base_r | base_p | T_SIN_I ;
551 return ;
552 }
553
554 // Protection
555 assert(tb->get_etat() == ETATQCQ) ;
556
557 // Pour le confort : nbre de points reels.
558 int nr = (tb->dim).dim[0] ;
559 int nt = (tb->dim).dim[1] ;
560 int np = (tb->dim).dim[2] ;
561 np = np - 2 ;
562
563 // zone de travail
564 double* som = new double [nr] ;
565
566 // pt. sur le tableau de double resultat
567 double* xo = new double[(tb->dim).taille] ;
568
569 // Initialisation a zero :
570 for (int i=0; i<(tb->dim).taille; i++) {
571 xo[i] = 0 ;
572 }
573
574 // On y va...
575 double* xi = tb->t ;
576 double* xci = xi ; // Pointeurs
577 double* xco = xo ; // courants
578
579 // k = 0
580
581 // Dernier j: j = nt-1
582 // Positionnement
583 xci += nr * nt ;
584 xco += nr * (nt-1) ;
585
586 // Initialisation de som
587 for (int i=0 ; i<nr ; i++) {
588 som[i] = 0. ;
589 xco[i] = 0. ;
590 }
591
592 // j suivants
593 for (int j=nt-2 ; j >= 0 ; j--) {
594 // Positionnement
595 xci -= nr ;
596 xco -= nr ;
597 // Calcul
598 for (int i=0 ; i<nr ; i++ ) {
599 som[i] += 2*xci[i] ;
600 xco[i] = som[i] ;
601 som[i] = -som[i] ;
602 } // Fin de la boucle sur r
603 } // Fin de la boucle sur theta
604 // Positionnement phi suivant
605 xci -= nr ;
606 xci += nr*nt ;
607 xco += nr*nt ;
608
609 // k = 1
610 xci += nr*nt ;
611 xco += nr*nt ;
612
613 // k >= 2
614 for (int k=2 ; k<np+1 ; k++) {
615
616 // Dernier j: j = nt-1
617 // Positionnement
618 xci += nr * nt ;
619 xco += nr * (nt-1) ;
620
621 // Initialisation de som
622 for (int i=0 ; i<nr ; i++) {
623 som[i] = 0. ;
624 xco[i] = 0. ;
625 }
626
627 // j suivants
628 for (int j=nt-2 ; j >= 0 ; j--) {
629 // Positionnement
630 xci -= nr ;
631 xco -= nr ;
632 // Calcul
633 for (int i=0 ; i<nr ; i++ ) {
634 som[i] += 2*xci[i] ;
635 xco[i] = som[i] ;
636 som[i] = -som[i] ;
637 } // Fin de la boucle sur r
638 } // Fin de la boucle sur theta
639 // Positionnement phi suivant
640 xci -= nr ;
641 xci += nr*nt ;
642 xco += nr*nt ;
643 } // Fin de la boucle sur phi
644
645 // On remet les choses la ou il faut
646 delete [] tb->t ;
647 tb->t = xo ;
648
649 //Menage
650 delete [] som ;
651
652 // base de developpement
653 int base_r = b & MSQ_R ;
654 int base_p = b & MSQ_P ;
655 b = base_r | base_p | T_SIN_I ;
656
657}
658 //--------------
659 // cas T_SIN_I
660 //--------------
661
662void _scost_t_sin_i(Tbl* tb, int & b)
663{
664
665
666 // Peut-etre rien a faire ?
667 if (tb->get_etat() == ETATZERO) {
668 int base_r = b & MSQ_R ;
669 int base_p = b & MSQ_P ;
670 b = base_r | base_p | T_SIN_P ;
671 return ;
672 }
673
674 // Protection
675 assert(tb->get_etat() == ETATQCQ) ;
676
677 // Pour le confort : nbre de points reels.
678 int nr = (tb->dim).dim[0] ;
679 int nt = (tb->dim).dim[1] ;
680 int np = (tb->dim).dim[2] ;
681 np = np - 2 ;
682
683 // zone de travail
684 double* som = new double [nr] ;
685
686 // pt. sur le tableau de double resultat
687 double* xo = new double[(tb->dim).taille] ;
688
689 // Initialisation a zero :
690 for (int i=0; i<(tb->dim).taille; i++) {
691 xo[i] = 0 ;
692 }
693
694 // On y va...
695 double* xi = tb->t ;
696 double* xci = xi ; // Pointeurs
697 double* xco = xo ; // courants
698
699
700 // k = 0
701
702 // Dernier j: j = nt-1
703 // Positionnement
704 xci += nr * (nt-1) ;
705 xco += nr * (nt-1) ;
706
707 // Initialisation de som
708 for (int i=0 ; i<nr ; i++) {
709 xco[i] = 0. ;
710 som[i] = 0. ;
711 }
712
713 // j suivants
714 for (int j=nt-2 ; j >= 1 ; j--) {
715 // Positionnement
716 xci -= nr ;
717 xco -= nr ;
718 // Calcul
719 for (int i=0 ; i<nr ; i++ ) {
720 som[i] += 2*xci[i] ;
721 xco[i] = som[i] ;
722 som[i] = -som[i] ;
723 } // Fin de la boucle sur r
724 } // Fin de la boucle sur theta
725 if (nt > 1) {
726 xci -= nr ;
727 xco -= nr ;
728 // Calcul pour le premier theta
729 for (int i=0 ; i<nr ; i++ ) {
730 xco[i] =0. ;
731 }
732 }
733 // Positionnement phi suivant
734 xci += nr*nt ;
735 xco += nr*nt ;
736
737 // k = 1
738 xci += nr*nt ;
739 xco += nr*nt ;
740
741 // k >= 2
742 for (int k=2 ; k<np+1 ; k++) {
743
744 // Dernier j: j = nt-1
745 // Positionnement
746 xci += nr * (nt-1) ;
747 xco += nr * (nt-1) ;
748
749 // Initialisation de som
750 for (int i=0 ; i<nr ; i++) {
751 xco[i] = 0. ;
752 som[i] = 0. ;
753 }
754
755 // j suivants
756 for (int j=nt-2 ; j >= 1 ; j--) {
757 // Positionnement
758 xci -= nr ;
759 xco -= nr ;
760 // Calcul
761 for (int i=0 ; i<nr ; i++ ) {
762 som[i] += 2*xci[i] ;
763 xco[i] = som[i] ;
764 som[i] = -som[i] ;
765 } // Fin de la boucle sur r
766 } // Fin de la boucle sur theta
767 xci -= nr ;
768 xco -= nr ;
769 // Calcul pour le premier theta
770 for (int i=0 ; i<nr ; i++ ) {
771 xco[i] =0. ;
772 }
773 // Positionnement phi suivant
774 xci += nr*nt ;
775 xco += nr*nt ;
776 } // Fin de la boucle sur phi
777
778
779 // On remet les choses la ou il faut
780 delete [] tb->t ;
781 tb->t = xo ;
782
783 //Menage
784 delete [] som ;
785
786 // base de developpement
787 int base_r = b & MSQ_R ;
788 int base_p = b & MSQ_P ;
789 b = base_r | base_p | T_SIN_P ;
790
791}
792 //---------------------
793 // cas T_COS_I
794 //---------------------
795void _scost_t_cos_i(Tbl* tb, int & b)
796{
797
798 // Peut-etre rien a faire ?
799 if (tb->get_etat() == ETATZERO) {
800 int base_r = b & MSQ_R ;
801 int base_p = b & MSQ_P ;
802 b = base_r | base_p | T_COS_P ;
803 return ;
804 }
805
806 // Protection
807 assert(tb->get_etat() == ETATQCQ) ;
808
809 // Pour le confort : nbre de points reels.
810 int nr = (tb->dim).dim[0] ;
811 int nt = (tb->dim).dim[1] ;
812 int np = (tb->dim).dim[2] ;
813 np = np - 2 ;
814
815 // zone de travail
816 double* som = new double [nr] ;
817
818 // pt. sur le tableau de double resultat
819 double* xo = new double[(tb->dim).taille] ;
820
821 // Initialisation a zero :
822 for (int i=0; i<(tb->dim).taille; i++) {
823 xo[i] = 0 ;
824 }
825
826 // On y va...
827 double* xi = tb->t ;
828 double* xci = xi ; // Pointeurs
829 double* xco = xo ; // courants
830
831 // k = 0
832
833 // Dernier j: j = nt-1
834 // Positionnement
835 xci += nr * (nt-1) ;
836 xco += nr * (nt-1) ;
837
838 // Initialisation de som
839 for (int i=0 ; i<nr ; i++) {
840 som[i] = 0. ;
841 xco[i] = 0. ; // mise a zero du dernier coef en theta
842 }
843
844 // j suivants
845 for (int j=nt-2 ; j >= 0 ; j--) {
846 // Positionnement
847 xci -= nr ;
848 xco -= nr ;
849 // Calcul
850 for (int i=0 ; i<nr ; i++ ) {
851 som[i] += 2*xci[i] ;
852 xco[i] = som[i] ;
853 som[i] = - som[i] ;
854 } // Fin de la boucle sur r
855 } // Fin de la boucle sur theta
856 // Normalisation du premier theta
857 for (int i=0 ; i<nr ; i++) {
858 xco[i] *= .5 ;
859 }
860 // Positionnement phi suivant
861 xci += nr*nt ;
862 xco += nr*nt ;
863
864 // k = 1
865 xci += nr*nt ;
866 xco += nr*nt ;
867
868 // k >= 2
869 for (int k=2 ; k<np+1 ; k++) {
870
871 // Dernier j: j = nt-1
872 // Positionnement
873 xci += nr * (nt-1) ;
874 xco += nr * (nt-1) ;
875
876 // Initialisation de som
877 for (int i=0 ; i<nr ; i++) {
878 som[i] = 0. ;
879 xco[i] = 0. ; // mise a zero du dernier coef en theta
880 }
881
882 // j suivants
883 for (int j=nt-2 ; j >= 0 ; j--) {
884 // Positionnement
885 xci -= nr ;
886 xco -= nr ;
887 // Calcul
888 for (int i=0 ; i<nr ; i++ ) {
889 som[i] += 2*xci[i] ;
890 xco[i] = som[i] ;
891 som[i] = -som[i] ;
892 } // Fin de la boucle sur r
893 } // Fin de la boucle sur theta
894 // Normalisation du premier theta
895 for (int i=0 ; i<nr ; i++) {
896 xco[i] *= .5 ;
897 }
898 // Positionnement phi suivant
899 xci += nr*nt ;
900 xco += nr*nt ;
901 } // Fin de la boucle sur phi
902
903 // On remet les choses la ou il faut
904 delete [] tb->t ;
905 tb->t = xo ;
906
907 //Menage
908 delete [] som ;
909
910 // base de developpement
911 int base_r = b & MSQ_R ;
912 int base_p = b & MSQ_P ;
913 b = base_r | base_p | T_COS_P ;
914
915}
916 //----------------------
917 // cas T_COSSIN_CP
918 //----------------------
919void _scost_t_cossin_cp(Tbl* tb, int & b)
920{
921
922 // Peut-etre rien a faire ?
923 if (tb->get_etat() == ETATZERO) {
924 int base_r = b & MSQ_R ;
925 int base_p = b & MSQ_P ;
926 b = base_r | base_p | T_COSSIN_CI ;
927 return ;
928 }
929
930 // Protection
931 assert(tb->get_etat() == ETATQCQ) ;
932
933 // Pour le confort : nbre de points reels.
934 int nr = (tb->dim).dim[0] ;
935 int nt = (tb->dim).dim[1] ;
936 int np = (tb->dim).dim[2] ;
937 np = np - 2 ;
938
939 // zone de travail
940 double* som = new double [nr] ;
941
942 // pt. sur le tableau de double resultat
943 double* xo = new double[(tb->dim).taille] ;
944
945 // Initialisation a zero :
946 for (int i=0; i<(tb->dim).taille; i++) {
947 xo[i] = 0 ;
948 }
949
950 // On y va...
951 double* xi = tb->t ;
952 double* xci = xi ; // Pointeurs
953 double* xco = xo ; // courants
954
955 // k = 0
956 int m = 0 ; // Parite de l'harmonique en phi
957
958 // Dernier j: j = nt-1
959 // Positionnement
960 xci += nr * (nt) ;
961 xco += nr * (nt-1) ;
962
963 // Initialisation de som
964 for (int i=0 ; i<nr ; i++) {
965 som[i] = 0. ;
966 xco[i] = 0. ; // mise a zero du dernier coef en theta
967 }
968
969 // j suivants
970 for (int j=nt-2 ; j >= 0 ; j--) {
971 // Positionnement
972 xci -= nr ;
973 xco -= nr ;
974 // Calcul
975 for (int i=0 ; i<nr ; i++ ) {
976 som[i] += 2*xci[i] ;
977 xco[i] = som[i] ;
978 som[i] = -som[i] ;
979 } // Fin de la boucle sur r
980 } // Fin de la boucle sur theta
981 // Positionnement phi suivant
982 xci -= nr ;
983 xci += nr*nt ;
984 xco += nr*nt ;
985
986 // k=1
987 xci += nr*nt ;
988 xco += nr*nt ;
989
990 // k >= 2
991 for (int k=2 ; k<np+1 ; k++) {
992 m = (k/2) % 2 ; // Parite de l'harmonique en phi
993
994 switch(m) {
995 case 0: // m pair: cos(pair)
996 // Dernier j: j = nt-1
997 // Positionnement
998 xci += nr * (nt) ;
999 xco += nr * (nt-1) ;
1000
1001 // Initialisation de som
1002 for (int i=0 ; i<nr ; i++) {
1003 som[i] = 0. ;
1004 xco[i] = 0. ; // mise a zero du dernier coef en theta
1005 }
1006
1007 // j suivants
1008 for (int j=nt-2 ; j >= 0 ; j--) {
1009 // Positionnement
1010 xci -= nr ;
1011 xco -= nr ;
1012 // Calcul
1013 for (int i=0 ; i<nr ; i++ ) {
1014 som[i] += 2*xci[i] ;
1015 xco[i] = som[i] ;
1016 som[i] = -som[i] ;
1017 } // Fin de la boucle sur r
1018 } // Fin de la boucle sur theta
1019 // Positionnement phi suivant
1020 xci -= nr ;
1021 xci += nr*nt ;
1022 xco += nr*nt ;
1023 break ;
1024
1025 case 1: // m impair: sin(impair)
1026 // Dernier j: j = nt-1
1027 // Positionnement
1028 xci += nr * (nt-1) ;
1029 xco += nr * (nt-1) ;
1030
1031 // Initialisation de som
1032 for (int i=0 ; i<nr ; i++) {
1033 xco[i] = 0. ;
1034 som[i] = 0. ;
1035 }
1036
1037 // j suivants
1038 for (int j=nt-2 ; j >= 1 ; j--) {
1039 // Positionnement
1040 xci -= nr ;
1041 xco -= nr ;
1042 // Calcul
1043 for (int i=0 ; i<nr ; i++ ) {
1044 som[i] += 2*xci[i] ;
1045 xco[i] = som[i] ;
1046 som[i] = -som[i] ;
1047 } // Fin de la boucle sur r
1048 } // Fin de la boucle sur theta
1049 xci -= nr ;
1050 xco -= nr ;
1051 // Calcul pour le premier theta
1052 for (int i=0 ; i<nr ; i++ ) {
1053 xco[i] =0. ;
1054 }
1055 // Positionnement phi suivant
1056 xci += nr*nt ;
1057 xco += nr*nt ;
1058 break;
1059 } // Fin du switch
1060 } // Fin de la boucle sur phi
1061
1062 // On remet les choses la ou il faut
1063 delete [] tb->t ;
1064 tb->t = xo ;
1065
1066 //Menage
1067 delete [] som ;
1068
1069 // base de developpement
1070 int base_r = b & MSQ_R ;
1071 int base_p = b & MSQ_P ;
1072 b = base_r | base_p | T_COSSIN_CI ;
1073}
1074
1075 //------------------
1076 // cas T_COSSIN_CI
1077 //----------------
1078void _scost_t_cossin_ci(Tbl* tb, int & b)
1079{
1080 // Peut-etre rien a faire ?
1081 if (tb->get_etat() == ETATZERO) {
1082 int base_r = b & MSQ_R ;
1083 int base_p = b & MSQ_P ;
1084 b = base_r | base_p | T_COSSIN_CP ;
1085 return ;
1086 }
1087
1088 // Protection
1089 assert(tb->get_etat() == ETATQCQ) ;
1090
1091 // Pour le confort : nbre de points reels.
1092 int nr = (tb->dim).dim[0] ;
1093 int nt = (tb->dim).dim[1] ;
1094 int np = (tb->dim).dim[2] ;
1095 np = np - 2 ;
1096
1097 // zone de travail
1098 double* som = new double [nr] ;
1099
1100 // pt. sur le tableau de double resultat
1101 double* xo = new double[(tb->dim).taille] ;
1102
1103 // Initialisation a zero :
1104 for (int i=0; i<(tb->dim).taille; i++) {
1105 xo[i] = 0 ;
1106 }
1107
1108 // On y va...
1109 double* xi = tb->t ;
1110 double* xci = xi ; // Pointeurs
1111 double* xco = xo ; // courants
1112
1113 // k = 0
1114 int m = 0 ; // Parite de l'harmonique en phi
1115
1116
1117 // Dernier j: j = nt-1
1118 // Positionnement
1119 xci += nr * (nt-1) ;
1120 xco += nr * (nt-1) ;
1121
1122 // Initialisation de som
1123 for (int i=0 ; i<nr ; i++) {
1124 som[i] = 0. ;
1125 xco[i] = 0. ; // mise a zero du dernier coef en theta
1126 }
1127
1128 // j suivants
1129 for (int j=nt-2 ; j >= 0 ; j--) {
1130 // Positionnement
1131 xci -= nr ;
1132 xco -= nr ;
1133 // Calcul
1134 for (int i=0 ; i<nr ; i++ ) {
1135 som[i] += 2*xci[i] ;
1136 xco[i] = som[i] ;
1137 som[i] = - som[i] ;
1138 } // Fin de la boucle sur r
1139 } // Fin de la boucle sur theta
1140 // Normalisation du premier theta
1141 for (int i=0 ; i<nr ; i++) {
1142 xco[i] *= .5 ;
1143 }
1144 // Positionnement phi suivant
1145 xci += nr*nt ;
1146 xco += nr*nt ;
1147
1148 // k=1
1149 xci += nr*nt ;
1150 xco += nr*nt ;
1151
1152 // k >= 2
1153 for (int k=2 ; k<np+1 ; k++) {
1154 m = (k/2) % 2 ; // Parite de l'harmonique en phi
1155
1156 switch(m) {
1157 case 0: // m pair: cos(impair)
1158 // Dernier j: j = nt-1
1159 // Positionnement
1160 xci += nr * (nt-1) ;
1161 xco += nr * (nt-1) ;
1162
1163 // Initialisation de som
1164 for (int i=0 ; i<nr ; i++) {
1165 som[i] = 0. ;
1166 xco[i] = 0. ; // mise a zero du dernier coef en theta
1167 }
1168
1169 // j suivants
1170 for (int j=nt-2 ; j >= 0 ; j--) {
1171 // Positionnement
1172 xci -= nr ;
1173 xco -= nr ;
1174 // Calcul
1175 for (int i=0 ; i<nr ; i++ ) {
1176 som[i] += 2*xci[i] ;
1177 xco[i] = som[i] ;
1178 som[i] = - som[i] ;
1179 } // Fin de la boucle sur r
1180 } // Fin de la boucle sur theta
1181 // Normalisation du premier theta
1182 for (int i=0 ; i<nr ; i++) {
1183 xco[i] *= .5 ;
1184 }
1185 // Positionnement phi suivant
1186 xci += nr*nt ;
1187 xco += nr*nt ;
1188 break ;
1189
1190 case 1:
1191 // m impair: sin(pair)
1192 // Dernier j: j = nt-1
1193 // Positionnement
1194 xci += nr * nt ;
1195 xco += nr * (nt-1) ;
1196
1197 // Initialisation de som
1198 for (int i=0 ; i<nr ; i++) {
1199 som[i] = 0. ;
1200 xco[i] = 0. ;
1201 }
1202
1203 // j suivants
1204 for (int j=nt-2 ; j >= 0 ; j--) {
1205 // Positionnement
1206 xci -= nr ;
1207 xco -= nr ;
1208 // Calcul
1209 for (int i=0 ; i<nr ; i++ ) {
1210 som[i] += 2*xci[i] ;
1211 xco[i] = som[i] ;
1212 som[i] = -som[i] ;
1213 } // Fin de la boucle sur r
1214 } // Fin de la boucle sur theta
1215 // Positionnement phi suivant
1216 xci -= nr ;
1217 xci += nr*nt ;
1218 xco += nr*nt ;
1219 break ;
1220 } // Fin du switch
1221 } // Fin de la boucle en phi
1222
1223 // On remet les choses la ou il faut
1224 delete [] tb->t ;
1225 tb->t = xo ;
1226
1227 //Menage
1228 delete [] som ;
1229
1230 // base de developpement
1231 int base_r = b & MSQ_R ;
1232 int base_p = b & MSQ_P ;
1233 b = base_r | base_p | T_COSSIN_CP ;
1234}
1235
1236 //---------------------
1237 // cas T_COSSIN_SI
1238 //----------------
1239void _scost_t_cossin_si(Tbl* tb, int & b)
1240{
1241 // Peut-etre rien a faire ?
1242 if (tb->get_etat() == ETATZERO) {
1243 int base_r = b & MSQ_R ;
1244 int base_p = b & MSQ_P ;
1245 b = base_r | base_p | T_COSSIN_SP ;
1246 return ;
1247 }
1248
1249 // Protection
1250 assert(tb->get_etat() == ETATQCQ) ;
1251
1252 // Pour le confort : nbre de points reels.
1253 int nr = (tb->dim).dim[0] ;
1254 int nt = (tb->dim).dim[1] ;
1255 int np = (tb->dim).dim[2] ;
1256 np = np - 2 ;
1257
1258 // zone de travail
1259 double* som = new double [nr] ;
1260
1261 // pt. sur le tableau de double resultat
1262 double* xo = new double[(tb->dim).taille] ;
1263
1264 // Initialisation a zero :
1265 for (int i=0; i<(tb->dim).taille; i++) {
1266 xo[i] = 0 ;
1267 }
1268
1269 // On y va...
1270 double* xi = tb->t ;
1271 double* xci = xi ; // Pointeurs
1272 double* xco = xo ; // courants
1273
1274 // k = 0
1275 int m = 0 ; // Parite de l'harmonique en phi
1276
1277 // Dernier j: j = nt-1
1278 // Positionnement
1279 xci += nr * (nt-1) ;
1280 xco += nr * (nt-1) ;
1281
1282 // Initialisation de som
1283 for (int i=0 ; i<nr ; i++) {
1284 xco[i] = 0. ;
1285 som[i] = 0. ;
1286 }
1287
1288 // j suivants
1289 for (int j=nt-2 ; j >= 1 ; j--) {
1290 // Positionnement
1291 xci -= nr ;
1292 xco -= nr ;
1293 // Calcul
1294 for (int i=0 ; i<nr ; i++ ) {
1295 som[i] += 2*xci[i] ;
1296 xco[i] = som[i] ;
1297 som[i] = -som[i] ;
1298 } // Fin de la boucle sur r
1299 } // Fin de la boucle sur theta
1300 xci -= nr ;
1301 xco -= nr ;
1302 // Calcul pour le premier theta
1303 for (int i=0 ; i<nr ; i++ ) {
1304 xco[i] =0. ;
1305 }
1306 // Positionnement phi suivant
1307 xci += nr*nt ;
1308 xco += nr*nt ;
1309
1310 // k=1
1311 xci += nr*nt ;
1312 xco += nr*nt ;
1313
1314 // k >= 2
1315 for (int k=2 ; k<np+1 ; k++) {
1316 m = (k/2) % 2 ; // Parite de l'harmonique en phi
1317
1318 switch(m) {
1319 case 0: // m pair: sin(impair)
1320 // Dernier j: j = nt-1
1321 // Positionnement
1322 xci += nr * (nt-1) ;
1323 xco += nr * (nt-1) ;
1324
1325 // Initialisation de som
1326 for (int i=0 ; i<nr ; i++) {
1327 xco[i] = 0. ;
1328 som[i] = 0. ;
1329 }
1330
1331 // j suivants
1332 for (int j=nt-2 ; j >= 1 ; j--) {
1333 // Positionnement
1334 xci -= nr ;
1335 xco -= nr ;
1336 // Calcul
1337 for (int i=0 ; i<nr ; i++ ) {
1338 som[i] += 2*xci[i] ;
1339 xco[i] = som[i] ;
1340 som[i] = -som[i] ;
1341 } // Fin de la boucle sur r
1342 } // Fin de la boucle sur theta
1343 xci -= nr ;
1344 xco -= nr ;
1345 // Calcul pour le premier theta
1346 for (int i=0 ; i<nr ; i++ ) {
1347 xco[i] =0. ;
1348 }
1349 // Positionnement phi suivant
1350 xci += nr*nt ;
1351 xco += nr*nt ;
1352 break ;
1353
1354 case 1: // m impair cos(pair)
1355 // Dernier j: j = nt-1
1356 // Positionnement
1357 xci += nr * (nt) ;
1358 xco += nr * (nt-1) ;
1359
1360 // Initialisation de som
1361 for (int i=0 ; i<nr ; i++) {
1362 som[i] = 0. ;
1363 xco[i] = 0. ; // mise a zero dui dernier coefficient en theta.
1364 }
1365
1366 // j suivants
1367 for (int j=nt-2 ; j >= 0 ; j--) {
1368 // Positionnement
1369 xci -= nr ;
1370 xco -= nr ;
1371 // Calcul
1372 for (int i=0 ; i<nr ; i++ ) {
1373 som[i] += 2*xci[i] ;
1374 xco[i] = som[i] ;
1375 som[i] = - som[i] ;
1376 } // Fin de la boucle sur r
1377 } // Fin de la boucle sur theta
1378 // Positionnement phi suivant
1379 xci -= nr ;
1380 xci += nr*nt ;
1381 xco += nr*nt ;
1382
1383 break ;
1384 } // Fin du switch
1385 } // Fin de la boucle en phi
1386
1387 // On remet les choses la ou il faut
1388 delete [] tb->t ;
1389 tb->t = xo ;
1390
1391 //Menage
1392 delete [] som ;
1393
1394 // base de developpement
1395 int base_r = b & MSQ_R ;
1396 int base_p = b & MSQ_P ;
1397 b = base_r | base_p | T_COSSIN_SP ;
1398
1399
1400}
1401 //---------------------
1402 // cas T_COSSIN_SP
1403 //----------------
1404void _scost_t_cossin_sp(Tbl* tb, int & b)
1405{
1406 // Peut-etre rien a faire ?
1407 if (tb->get_etat() == ETATZERO) {
1408 int base_r = b & MSQ_R ;
1409 int base_p = b & MSQ_P ;
1410 b = base_r | base_p | T_COSSIN_SI ;
1411 return ;
1412 }
1413
1414 // Protection
1415 assert(tb->get_etat() == ETATQCQ) ;
1416
1417 // Pour le confort : nbre de points reels.
1418 int nr = (tb->dim).dim[0] ;
1419 int nt = (tb->dim).dim[1] ;
1420 int np = (tb->dim).dim[2] ;
1421 np = np - 2 ;
1422
1423 // zone de travail
1424 double* som = new double [nr] ;
1425
1426 // pt. sur le tableau de double resultat
1427 double* xo = new double[(tb->dim).taille] ;
1428
1429 // Initialisation a zero :
1430 for (int i=0; i<(tb->dim).taille; i++) {
1431 xo[i] = 0 ;
1432 }
1433
1434 // On y va...
1435 double* xi = tb->t ;
1436 double* xci = xi ; // Pointeurs
1437 double* xco = xo ; // courants
1438
1439 // k = 0
1440 int m = 0 ; // Parite de l'harmonique en phi
1441
1442 // Dernier j: j = nt-1
1443 // Positionnement
1444 xci += nr * nt ;
1445 xco += nr * (nt-1) ;
1446
1447 // Initialisation de som
1448 for (int i=0 ; i<nr ; i++) {
1449 som[i] = 0. ;
1450 xco[i] = 0. ;
1451 }
1452
1453 // j suivants
1454 for (int j=nt-2 ; j >= 0 ; j--) {
1455 // Positionnement
1456 xci -= nr ;
1457 xco -= nr ;
1458 // Calcul
1459 for (int i=0 ; i<nr ; i++ ) {
1460 som[i] += 2*xci[i] ;
1461 xco[i] = som[i] ;
1462 som[i] = -som[i] ;
1463 } // Fin de la boucle sur r
1464 } // Fin de la boucle sur theta
1465 // Positionnement phi suivant
1466 xci -= nr ;
1467 xci += nr*nt ;
1468 xco += nr*nt ;
1469
1470 // k=1
1471 xci += nr*nt ;
1472 xco += nr*nt ;
1473
1474 for (int k=2 ; k<np+1 ; k++) {
1475 m = (k/2) % 2 ; // Parite de l'harmonique en phi
1476
1477 switch(m) {
1478 case 1: // m impair: cos(impair)
1479 // Dernier j: j = nt-1
1480 // Positionnement
1481 xci += nr * (nt-1) ;
1482 xco += nr * (nt-1) ;
1483
1484 // Initialisation de som
1485 for (int i=0 ; i<nr ; i++) {
1486 som[i] = 0. ;
1487 xco[i] = 0. ; // mise a zero du dernier coef en theta
1488 }
1489
1490 // j suivants
1491 for (int j=nt-2 ; j >= 0 ; j--) {
1492 // Positionnement
1493 xci -= nr ;
1494 xco -= nr ;
1495 // Calcul
1496 for (int i=0 ; i<nr ; i++ ) {
1497 som[i] += 2*xci[i] ;
1498 xco[i] = som[i] ;
1499 som[i] = -som[i] ;
1500 } // Fin de la boucle sur r
1501 } // Fin de la boucle sur theta
1502 // Normalisation du premier theta
1503 for (int i=0 ; i<nr ; i++) {
1504 xco[i] *= .5 ;
1505 }
1506 // Positionnement phi suivant
1507 xci += nr*nt ;
1508 xco += nr*nt ;
1509 break ;
1510
1511 case 0: // m pair: sin(pair)
1512 // Dernier j: j = nt-1
1513 // Positionnement
1514 xci += nr * nt ;
1515 xco += nr * (nt-1) ;
1516
1517 // Initialisation de som
1518 for (int i=0 ; i<nr ; i++) {
1519 som[i] = 0. ;
1520 xco[i] = 0. ;
1521 }
1522
1523 // j suivants
1524 for (int j=nt-2 ; j >= 0 ; j--) {
1525 // Positionnement
1526 xci -= nr ;
1527 xco -= nr ;
1528 // Calcul
1529 for (int i=0 ; i<nr ; i++ ) {
1530 som[i] += 2*xci[i] ;
1531 xco[i] = som[i] ;
1532 som[i] = -som[i] ;
1533 } // Fin de la boucle sur r
1534 } // Fin de la boucle sur theta
1535 // Positionnement phi suivant
1536 xci -= nr ;
1537 xci += nr*nt ;
1538 xco += nr*nt ;
1539 break ;
1540 } // Fin du switch
1541 } // Fin de la boucle en phi
1542
1543 // On remet les choses la ou il faut
1544 delete [] tb->t ;
1545 tb->t = xo ;
1546
1547 //Menage
1548 delete [] som ;
1549
1550 // base de developpement
1551 int base_r = b & MSQ_R ;
1552 int base_p = b & MSQ_P ;
1553 b = base_r | base_p | T_COSSIN_SI ;
1554
1555}
1556
1557 //----------------------
1558 // cas T_COSSIN_C
1559 //----------------------
1560void _scost_t_cossin_c(Tbl* tb, int & b) {
1561
1562 // Peut-etre rien a faire ?
1563 if (tb->get_etat() == ETATZERO) {
1564 int base_r = b & MSQ_R ;
1565 int base_p = b & MSQ_P ;
1566 switch(base_r){
1567 case(R_CHEBPI_P):
1568 b = R_CHEBPI_I | base_p | T_COSSIN_C ;
1569 break ;
1570 case(R_CHEBPI_I):
1571 b = R_CHEBPI_P | base_p | T_COSSIN_C ;
1572 break ;
1573 default:
1574 b = base_r | base_p | T_COSSIN_C ;
1575 break;
1576 }
1577 return ;
1578 }
1579
1580 // Protection
1581 assert(tb->get_etat() == ETATQCQ) ;
1582
1583 // Pour le confort : nbre de points reels.
1584 int nr = (tb->dim).dim[0] ;
1585 int nt = (tb->dim).dim[1] ;
1586 int np = (tb->dim).dim[2] ;
1587 np = np - 2 ;
1588
1589 // zone de travail
1590 double* somP = new double [nr] ;
1591 double* somN = new double [nr] ;
1592
1593 // pt. sur le tableau de double resultat
1594 double* xo = new double[(tb->dim).taille] ;
1595
1596 // Initialisation a zero :
1597 for (int i=0; i<(tb->dim).taille; i++) xo[i] = 0 ;
1598
1599 // On y va...
1600 double* xi = tb->t ;
1601 double* xci = xi ; // Pointeurs
1602 double* xco = xo ; // courants
1603
1604 // k = 0
1605
1606 // Dernier j: j = nt-1
1607 // Positionnement
1608 xci += nr * (nt-1) ;
1609 xco += nr * (nt-1) ;
1610
1611 // Initialisation de som
1612 for (int i=0 ; i<nr ; i++) {
1613 somP[i] = 0. ;
1614 somN[i] = 0. ;
1615 xco[i] = 0. ; // mise a zero du dernier coef en theta
1616 }
1617
1618 // j suivants
1619 for (int j=nt-2 ; j >= 0 ; j--) {
1620 // Positionnement
1621 xci -= nr ;
1622 xco -= nr ;
1623 // Calcul
1624 for (int i=0 ; i<nr ; i++ ) {
1625 if((j%2) == 1) {
1626 somN[i] = -somN[i] ;
1627 somN[i] += 2*xci[i] ;
1628 xco[i] = somP[i] ;
1629 }
1630 else {
1631 somP[i] = -somP[i] ;
1632 somP[i] += 2*xci[i] ;
1633 xco[i] = somN[i] ;
1634 }
1635 } // Fin de la boucle sur r
1636 } // Fin de la boucle sur theta
1637 // j=0 : le facteur 2...
1638 for (int i=0 ; i<nr ; i++) xco[i] *= .5 ;
1639
1640 // Positionnement phi suivant
1641 xci += nr*nt ;
1642 xco += nr*nt ;
1643
1644 // k=1
1645 xci += nr*nt ;
1646 xco += nr*nt ;
1647
1648 // k >= 2
1649 for (int k=2 ; k<np+1 ; k++) {
1650
1651 // Dernier j: j = nt-1
1652 // Positionnement
1653 xco += nr * (nt-1) ;
1654 xci += nr * (nt-1) ;
1655
1656 // Initialisation de som
1657 for (int i=0 ; i<nr ; i++) {
1658 somP[i] = 0. ;
1659 somN[i] = 0. ;
1660 xco[i] = 0. ;
1661 }
1662
1663 // j suivants
1664 for (int j=nt-2 ; j >= 0 ; j--) {
1665 // Positionnement
1666 xci -= nr ;
1667 xco -= nr ;
1668 // Calcul
1669 for (int i=0 ; i<nr ; i++ ) {
1670 if((j%2) == 1) {
1671 somN[i] = -somN[i];
1672 somN[i] += 2 * xci[i] ;
1673 xco[i] = somP[i] ;
1674 }
1675 else {
1676 somP[i] = -somP[i];
1677 somP[i] += 2 * xci[i] ;
1678 xco[i] = somN[i];
1679 }
1680 } // Fin de la boucle sur r
1681 } // Fin de la boucle sur theta
1682 double fac_m = ( (k/2)%2 == 1 ? 0. : 0.5) ;
1683 // j=0 : sin(0*theta) ou facteur 2, suivant la parite de m
1684 for (int i=0 ; i<nr ; i++) xco[i] *= fac_m ;
1685
1686 // Positionnement phi suivant
1687 xci += nr*nt ;
1688 xco += nr*nt ;
1689 } // Fin de la boucle sur phi
1690
1691 // On remet les choses la ou il faut
1692 delete [] tb->t ;
1693 tb->t = xo ;
1694
1695 //Menage
1696 delete [] somP ;
1697 delete [] somN ;
1698
1699 // base de developpement
1700 int base_r = b & MSQ_R ;
1701 int base_p = b & MSQ_P ;
1702 switch(base_r){
1703 case(R_CHEBPI_P):
1704 b = R_CHEBPI_I | base_p | T_COSSIN_C ;
1705 break ;
1706 case(R_CHEBPI_I):
1707 b = R_CHEBPI_P | base_p | T_COSSIN_C ;
1708 break ;
1709 default:
1710 b = base_r | base_p | T_COSSIN_C ;
1711 break;
1712 }
1713}
1714
1715 //---------------------
1716 // cas T_COSSIN_S
1717 //----------------
1718void _scost_t_cossin_s(Tbl* tb, int & b) {
1719
1720 // Peut-etre rien a faire ?
1721 if (tb->get_etat() == ETATZERO) {
1722 int base_r = b & MSQ_R ;
1723 int base_p = b & MSQ_P ;
1724 switch(base_r){
1725 case(R_CHEBPI_P):
1726 b = R_CHEBPI_I | base_p | T_COSSIN_S ;
1727 break ;
1728 case(R_CHEBPI_I):
1729 b = R_CHEBPI_P | base_p | T_COSSIN_S ;
1730 break ;
1731 default:
1732 b = base_r | base_p | T_COSSIN_S ;
1733 break;
1734 }
1735 return ;
1736 }
1737
1738 // Protection
1739 assert(tb->get_etat() == ETATQCQ) ;
1740
1741 // Pour le confort : nbre de points reels.
1742 int nr = (tb->dim).dim[0] ;
1743 int nt = (tb->dim).dim[1] ;
1744 int np = (tb->dim).dim[2] ;
1745 np = np - 2 ;
1746
1747 // zone de travail
1748 double* somP = new double [nr] ;
1749 double* somN = new double [nr] ;
1750
1751 // pt. sur le tableau de double resultat
1752 double* xo = new double[(tb->dim).taille] ;
1753
1754 // Initialisation a zero :
1755 for (int i=0; i<(tb->dim).taille; i++) xo[i] = 0 ;
1756
1757 // On y va...
1758 double* xi = tb->t ;
1759 double* xci = xi ; // Pointeurs
1760 double* xco = xo ; // courants
1761
1762 // k = 0
1763
1764 // Dernier j: j = nt-1
1765 // Positionnement
1766 xci += nr * (nt-1) ;
1767 xco += nr * (nt-1) ;
1768
1769 // Initialisation de som
1770 for (int i=0 ; i<nr ; i++) {
1771 somP[i] = 0. ;
1772 somN[i] = 0. ;
1773 xco[i] = 0. ; // mise a zero du dernier coef en theta
1774 }
1775
1776 // j suivants
1777 for (int j=nt-2 ; j >= 0 ; j--) {
1778 // Positionnement
1779 xci -= nr ;
1780 xco -= nr ;
1781 // Calcul
1782 for (int i=0 ; i<nr ; i++ ) {
1783 if((j%2) == 1) {
1784 somN[i] = -somN[i] ;
1785 somN[i] += 2*xci[i] ;
1786 xco[i] = somP[i] ;
1787 }
1788 else {
1789 somP[i] = -somP[i] ;
1790 somP[i] += 2*xci[i] ;
1791 xco[i] = somN[i] ;
1792 }
1793 } // Fin de la boucle sur r
1794 } // Fin de la boucle sur theta
1795 // j=0 : sin(0*theta)
1796 for (int i=0 ; i<nr ; i++) xco[i] = 0. ;
1797
1798 // Positionnement phi suivant
1799 xci += nr*nt ;
1800 xco += nr*nt ;
1801
1802 // k=1
1803 xci += nr*nt ;
1804 xco += nr*nt ;
1805
1806 // k >= 2
1807 for (int k=2 ; k<np+1 ; k++) {
1808
1809 // Dernier j: j = nt-1
1810 // Positionnement
1811 xco += nr * (nt-1) ;
1812 xci += nr * (nt-1) ;
1813
1814 // Initialisation de som
1815 for (int i=0 ; i<nr ; i++) {
1816 somP[i] = 0. ;
1817 somN[i] = 0. ;
1818 xco[i] = 0. ;
1819 }
1820
1821 // j suivants
1822 for (int j=nt-2 ; j >= 0 ; j--) {
1823 // Positionnement
1824 xci -= nr ;
1825 xco -= nr ;
1826 // Calcul
1827 for (int i=0 ; i<nr ; i++ ) {
1828 if((j%2) == 1) {
1829 somN[i] = -somN[i] ;
1830 somN[i] += 2*xci[i] ;
1831 xco[i] = somP[i] ;
1832 }
1833 else {
1834 somP[i] = -somP[i] ;
1835 somP[i] += 2*xci[i] ;
1836 xco[i] = somN[i] ;
1837 }
1838 } // Fin de la boucle sur r
1839 } // Fin de la boucle sur theta
1840
1841 double fac_m = ( (k/2)%2 == 0 ? 0. : 0.5) ;
1842 // j=0 : sin(0*theta) ou facteur 2, suivant la parite de m
1843 for (int i=0 ; i<nr ; i++) xco[i] *= fac_m ;
1844
1845 // Positionnement phi suivant
1846 xci += nr*nt ;
1847 xco += nr*nt ;
1848
1849 } // Fin de la boucle sur phi
1850
1851 // On remet les choses la ou il faut
1852 delete [] tb->t ;
1853 tb->t = xo ;
1854
1855 //Menage
1856 delete [] somP ;
1857 delete [] somN ;
1858
1859 // base de developpement
1860 int base_r = b & MSQ_R ;
1861 int base_p = b & MSQ_P ;
1862 switch(base_r){
1863 case(R_CHEBPI_P):
1864 b = R_CHEBPI_I | base_p | T_COSSIN_S ;
1865 break ;
1866 case(R_CHEBPI_I):
1867 b = R_CHEBPI_P | base_p | T_COSSIN_S ;
1868 break ;
1869 default:
1870 b = base_r | base_p | T_COSSIN_S ;
1871 break;
1872 }
1873}
1874}
Lorene::Tbl
Basic array class.
Definition tbl.h:161
T_COSSIN_SP
#define T_COSSIN_SP
sin pair-cos impair alternes, sin pour m=0
Definition type_parite.h:210
T_SIN_P
#define T_SIN_P
dev. sin seulement, harmoniques paires
Definition type_parite.h:202
MSQ_R
#define MSQ_R
Extraction de l'info sur R.
Definition type_parite.h:152
T_COSSIN_S
#define T_COSSIN_S
dev. cos-sin alternes, sin pour m=0
Definition type_parite.h:194
T_COSSIN_SI
#define T_COSSIN_SI
sin impair-cos pair alternes, sin pour m=0
Definition type_parite.h:214
R_CHEBPI_I
#define R_CHEBPI_I
Cheb. pair-impair suivant l impair pour l=0.
Definition type_parite.h:174
T_COS_P
#define T_COS_P
dev. cos seulement, harmoniques paires
Definition type_parite.h:200
T_COSSIN_CI
#define T_COSSIN_CI
cos impair-sin pair alternes, cos pour m=0
Definition type_parite.h:212
T_COSSIN_CP
#define T_COSSIN_CP
cos pair-sin impair alternes, cos pour m=0
Definition type_parite.h:208
T_SIN_I
#define T_SIN_I
dev. sin seulement, harmoniques impaires
Definition type_parite.h:206
T_COS
#define T_COS
dev. cos seulement
Definition type_parite.h:196
MSQ_P
#define MSQ_P
Extraction de l'info sur Phi.
Definition type_parite.h:156
T_SIN
#define T_SIN
dev. sin seulement
Definition type_parite.h:198
T_COS_I
#define T_COS_I
dev. cos seulement, harmoniques impaires
Definition type_parite.h:204
T_COSSIN_C
#define T_COSSIN_C
dev. cos-sin alternes, cos pour m=0
Definition type_parite.h:192
R_CHEBPI_P
#define R_CHEBPI_P
Cheb. pair-impair suivant l pair pour l=0.
Definition type_parite.h:172
Lorene
Lorene prototypes.
Definition app_hor.h:67