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