LORENE
op_mult_ct.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_mult_ct_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_mult_ct.C,v 1.8 2014/10/13 08:53:25 j_novak Exp $" ;
24
25/*
26 * $Id: op_mult_ct.C,v 1.8 2014/10/13 08:53:25 j_novak Exp $
27 * $Log: op_mult_ct.C,v $
28 * Revision 1.8 2014/10/13 08:53:25 j_novak
29 * Lorene classes and functions now belong to the namespace Lorene.
30 *
31 * Revision 1.7 2013/04/25 15:46:06 j_novak
32 * Added special treatment in the case np = 1, for type_p = NONSYM.
33 *
34 * Revision 1.6 2009/10/09 14:00:54 j_novak
35 * New bases T_cos and T_SIN.
36 *
37 * Revision 1.5 2007/12/21 10:43:37 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:23 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:41:51 eric
56 * Initialisation a zero du tableau xo avant le calcul.
57 *
58 * Revision 1.2 1999/11/23 16:14:09 novak
59 * *** empty log message ***
60 *
61 * Revision 1.1 1999/11/23 14:31:56 novak
62 * Initial revision
63 *
64 *
65 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_mult_ct.C,v 1.8 2014/10/13 08:53:25 j_novak Exp $
66 *
67 */
68
69/*
70 * Ensemble des routines de base de multiplication par cost theta
71 * (Utilisation interne)
72 *
73 * void _mult_ct_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 _mult_ct_pas_prevu(Tbl * tb, int& base) {
90 cout << "mult_ct 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 _mult_ct_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 // zone de travail
131 double* som = new double [nr] ;
132
133 // pt. sur le tableau de double resultat
134 double* xo = new double[(tb->dim).taille] ;
135
136 // Initialisation a zero :
137 for (int i=0; i<(tb->dim).taille; i++) {
138 xo[i] = 0 ;
139 }
140
141 // On y va...
142 double* xi = tb->t ;
143 double* xci = xi ; // Pointeurs
144 double* xco = xo ; // courants
145
146 // k = 0
147
148 // Dernier j: j = nt-1
149 // Positionnement
150 xci += nr * (nt-1) ;
151 xco += nr * (nt-1) ;
152
153 // Initialisation de som
154 for (int i=0 ; i<nr ; i++) {
155 som[i] = 0.5*xci[i] ;
156 xco[i] = 0. ; // mise a zero du dernier coef en theta
157 }
158
159 // j suivants
160 for (int j=nt-2 ; j > 0 ; j--) {
161 // Positionnement
162 xci -= 2*nr ;
163 xco -= nr ;
164 // Calcul
165 for (int i=0 ; i<nr ; i++ ) {
166 som[i] += 0.5*xci[i] ;
167 xco[i] = som[i] ;
168 } // Fin de la boucle sur r
169 xci += nr ;
170 for (int i=0; i<nr; i++)
171 som[i] = 0.5*xci[i] ;
172 } // Fin de la boucle sur theta
173 // j = 0 : le facteur 2...
174 xci -= nr ;
175 for (int i=0 ; i<nr ; i++ ) {
176 xco[i] += 0.5*xci[i] ;
177 }
178 xco -= nr ;
179 for (int i=0; i<nr; i++)
180 xco[i] = som[i] ;
181
182 // Positionnement phi suivant
183 xci += nr*nt ;
184 xco += nr*nt ;
185
186 // k = 1
187 xci += nr*nt ;
188 xco += nr*nt ;
189
190 // k >= 2
191 for (int k=2 ; k<np+1 ; k++) {
192
193 // Dernier j: j = nt-1
194 // Positionnement
195 xci += nr * (nt-1) ;
196 xco += nr * (nt-1) ;
197
198 // Initialisation de som
199 for (int i=0 ; i<nr ; i++) {
200 som[i] = 0.5*xci[i] ;
201 xco[i] = 0. ; // mise a zero dui dernier coefficient en theta.
202 }
203
204 // j suivants
205 for (int j=nt-2 ; j > 0 ; j--) {
206 // Positionnement
207 xci -= 2*nr ;
208 xco -= nr ;
209 // Calcul
210 for (int i=0 ; i<nr ; i++ ) {
211 som[i] += 0.5*xci[i] ;
212 xco[i] = som[i] ;
213 } // Fin de la boucle sur r
214 xci += nr ;
215 for (int i=0; i<nr; i++)
216 som[i] = 0.5*xci[i] ;
217 } // Fin de la boucle sur theta
218 // j = 0 : le facteur 2...
219 xci -= nr ;
220 for (int i = 0; i<nr; i++) {
221 xco[i] += 0.5*xci[i] ;
222 }
223 xco -= nr ;
224 for (int i=0; i<nr; i++)
225 xco[i] = som[i] ;
226 // Positionnement phi suivant
227 xci += nr*nt ;
228 xco += nr*nt ;
229 } // Fin de la boucle sur phi
230
231 // On remet les choses la ou il faut
232 delete [] tb->t ;
233 tb->t = xo ;
234
235 //Menage
236 delete [] som ;
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 _mult_ct_t_sin(Tbl* tb, int & b)
259{
260 // Peut-etre rien a faire ?
261 if (tb->get_etat() == ETATZERO) {
262 int base_r = b & MSQ_R ;
263 int base_p = b & MSQ_P ;
264 switch(base_r){
265 case(R_CHEBPI_P):
266 b = R_CHEBPI_I | base_p | T_SIN ;
267 break ;
268 case(R_CHEBPI_I):
269 b = R_CHEBPI_P | base_p | T_SIN ;
270 break ;
271 default:
272 b = base_r | base_p | T_SIN ;
273 break;
274 }
275 return ;
276 }
277
278 // Protection
279 assert(tb->get_etat() == ETATQCQ) ;
280
281 // Pour le confort : nbre de points reels.
282 int nr = (tb->dim).dim[0] ;
283 int nt = (tb->dim).dim[1] ;
284 int np = (tb->dim).dim[2] ;
285 np = np - 2 ;
286
287 // zone de travail
288 double* som = new double [nr] ;
289
290 // pt. sur le tableau de double resultat
291 double* xo = new double[(tb->dim).taille] ;
292
293 // Initialisation a zero :
294 for (int i=0; i<(tb->dim).taille; i++) {
295 xo[i] = 0 ;
296 }
297
298 // On y va...
299 double* xi = tb->t ;
300 double* xci = xi ; // Pointeurs
301 double* xco = xo ; // courants
302
303 // k = 0
304
305 // Dernier j: j = nt-1
306 // Positionnement
307 xci += nr * (nt-1) ;
308 xco += nr * (nt-1) ;
309
310 // Initialisation de som
311 for (int i=0 ; i<nr ; i++) {
312 som[i] = 0.5*xci[i] ;
313 xco[i] = 0. ;
314 }
315
316 // j suivants
317 for (int j=nt-2 ; j > 0 ; j--) {
318 // Positionnement
319 xci -= 2*nr ;
320 xco -= nr ;
321 // Calcul
322 for (int i=0 ; i<nr ; i++ ) {
323 som[i] += 0.5*xci[i] ;
324 xco[i] = som[i] ;
325 } // Fin de la boucle sur r
326 xci += nr ;
327 for (int i=0; i<nr; i++) {
328 som[i] = 0.5*xci[i] ;
329 }
330 } // Fin de la boucle sur theta
331 // j = 0 : sin(0*theta)
332 xci -= nr ;
333 xco -= nr ;
334 for (int i =0; i<nr ; i++) {
335 xco[i] = 0. ;
336 }
337 // Positionnement phi suivant
338 xci += nr*nt ;
339 xco += nr*nt ;
340
341 // k = 1
342 xci += nr*nt ;
343 xco += nr*nt ;
344
345 // k >= 2
346 for (int k=2 ; k<np+1 ; k++) {
347
348 // Dernier j: j = nt-1
349 // Positionnement
350 xci += nr * (nt-1) ;
351 xco += nr * (nt-1) ;
352
353 // Initialisation de som
354 for (int i=0 ; i<nr ; i++) {
355 som[i] = 0.5*xci[i] ;
356 xco[i] = 0. ;
357 }
358
359 // j suivants
360 for (int j=nt-2 ; j > 0 ; j--) {
361 // Positionnement
362 xci -= 2*nr ;
363 xco -= nr ;
364 // Calcul
365 for (int i=0 ; i<nr ; i++ ) {
366 som[i] += 0.5*xci[i] ;
367 xco[i] = som[i] ;
368 } // Fin de la boucle sur r
369 xci += nr ;
370 for (int i=0; i<nr; i++) {
371 som[i] = 0.5*xci[i] ;
372 }
373 } // Fin de la boucle sur theta
374 // j = 0 : sin (0*theta)
375 xci -= nr ;
376 xco -= nr ;
377 for (int i=0; i<nr; i++) {
378 xco[i] = 0.0 ;
379 }
380 // Positionnement phi suivant
381 xci += nr*nt ;
382 xco += nr*nt ;
383 } // Fin de la boucle sur phi
384
385 // On remet les choses la ou il faut
386 delete [] tb->t ;
387 tb->t = xo ;
388
389 //Menage
390 delete [] som ;
391
392 // base de developpement
393 int base_r = b & MSQ_R ;
394 int base_p = b & MSQ_P ;
395 switch(base_r){
396 case(R_CHEBPI_P):
397 b = R_CHEBPI_I | base_p | T_SIN ;
398 break ;
399 case(R_CHEBPI_I):
400 b = R_CHEBPI_P | base_p | T_SIN ;
401 break ;
402 default:
403 b = base_r | base_p | T_SIN ;
404 break;
405 }
406}
407 //----------------
408 // cas T_COS_P
409 //----------------
410
411void _mult_ct_t_cos_p(Tbl* tb, int & b)
412{
413
414 // Peut-etre rien a faire ?
415 if (tb->get_etat() == ETATZERO) {
416 int base_r = b & MSQ_R ;
417 int base_p = b & MSQ_P ;
418 b = base_r | base_p | T_COS_I ;
419 return ;
420 }
421
422 // Protection
423 assert(tb->get_etat() == ETATQCQ) ;
424
425 // Pour le confort : nbre de points reels.
426 int nr = (tb->dim).dim[0] ;
427 int nt = (tb->dim).dim[1] ;
428 int np = (tb->dim).dim[2] ;
429 np = np - 2 ;
430
431 // zone de travail
432 double* som = new double [nr] ;
433
434 // pt. sur le tableau de double resultat
435 double* xo = new double[(tb->dim).taille] ;
436
437 // Initialisation a zero :
438 for (int i=0; i<(tb->dim).taille; i++) {
439 xo[i] = 0 ;
440 }
441
442 // On y va...
443 double* xi = tb->t ;
444 double* xci = xi ; // Pointeurs
445 double* xco = xo ; // courants
446
447 // k = 0
448
449 // Dernier j: j = nt-1
450 // Positionnement
451 xci += nr * (nt-1) ;
452 xco += nr * (nt-1) ;
453
454 // Initialisation de som
455 for (int i=0 ; i<nr ; i++) {
456 som[i] = 0.5*xci[i] ;
457 xco[i] = 0. ; // mise a zero du dernier coef en theta
458 }
459
460 // j suivants
461 for (int j=nt-2 ; j > 0 ; j--) {
462 // Positionnement
463 xci -= nr ;
464 xco -= nr ;
465 // Calcul
466 for (int i=0 ; i<nr ; i++ ) {
467 som[i] += 0.5*xci[i] ;
468 xco[i] = som[i] ;
469 som[i] = 0.5*xci[i] ;
470 } // Fin de la boucle sur r
471 } // Fin de la boucle sur theta
472
473 if (nt > 1) {
474 // j = 0
475 xci -= nr ;
476 xco -= nr ;
477 for (int i=0 ; i<nr ; i++ ) {
478 xco[i] = som[i]+xci[i] ;
479 } // Fin de la boucle sur r
480 }
481
482 // Positionnement phi suivant
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-1) ;
496 xco += nr * (nt-1) ;
497
498 // Initialisation de som
499 for (int i=0 ; i<nr ; i++) {
500 som[i] = 0.5*xci[i] ;
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] += 0.5*xci[i] ;
512 xco[i] = som[i] ;
513 som[i] = 0.5*xci[i] ;
514 } // Fin de la boucle sur r
515 } // Fin de la boucle sur theta
516 // j = 0
517 xci -= nr ;
518 xco -= nr ;
519 for (int i = 0; i<nr; i++) {
520 xco[i] = xci[i] + som[i] ;
521 }
522 // Positionnement phi suivant
523 xci += nr*nt ;
524 xco += nr*nt ;
525 } // Fin de la boucle sur phi
526
527 // On remet les choses la ou il faut
528 delete [] tb->t ;
529 tb->t = xo ;
530
531 //Menage
532 delete [] som ;
533
534 // base de developpement
535 int base_r = b & MSQ_R ;
536 int base_p = b & MSQ_P ;
537 b = base_r | base_p | T_COS_I ;
538}
539
540 //--------------
541 // cas T_SIN_P
542 //--------------
543
544void _mult_ct_t_sin_p(Tbl* tb, int & b)
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-1) ;
584 xco += nr * (nt-1) ;
585
586 // Initialisation de som
587 for (int i=0 ; i<nr ; i++) {
588 som[i] = 0.5*xci[i] ;
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] += 0.5*xci[i] ;
600 xco[i] = som[i] ;
601 som[i] = 0.5*xci[i] ;
602 } // Fin de la boucle sur r
603 } // Fin de la boucle sur theta
604 // j = 0
605 if (nt > 1) {
606 xci -= nr ;
607 xco -= nr ;
608 for (int i =0; i<nr ; i++) {
609 xco[i] = som[i] ;
610 }
611 }
612
613 // Positionnement phi suivant
614 xci += nr*nt ;
615 xco += nr*nt ;
616
617 // k = 1
618 xci += nr*nt ;
619 xco += nr*nt ;
620
621 // k >= 2
622 for (int k=2 ; k<np+1 ; k++) {
623
624 // Dernier j: j = nt-1
625 // Positionnement
626 xci += nr * (nt-1) ;
627 xco += nr * (nt-1) ;
628
629 // Initialisation de som
630 for (int i=0 ; i<nr ; i++) {
631 som[i] = 0.5*xci[i] ;
632 xco[i] = 0. ;
633 }
634
635 // j suivants
636 for (int j=nt-2 ; j > 0 ; j--) {
637 // Positionnement
638 xci -= nr ;
639 xco -= nr ;
640 // Calcul
641 for (int i=0 ; i<nr ; i++ ) {
642 som[i] += 0.5*xci[i] ;
643 xco[i] = som[i] ;
644 som[i] = 0.5*xci[i] ;
645 } // Fin de la boucle sur r
646 } // Fin de la boucle sur theta
647 // j = 0
648 xci -= nr ;
649 xco -= nr ;
650 for (int i=0; i<nr; i++) {
651 xco[i] = som[i] ;
652 }
653 // Positionnement phi suivant
654 xci += nr*nt ;
655 xco += nr*nt ;
656 } // Fin de la boucle sur phi
657
658 // On remet les choses la ou il faut
659 delete [] tb->t ;
660 tb->t = xo ;
661
662 //Menage
663 delete [] som ;
664
665 // base de developpement
666 int base_r = b & MSQ_R ;
667 int base_p = b & MSQ_P ;
668 b = base_r | base_p | T_SIN_I ;
669
670}
671 //--------------
672 // cas T_SIN_I
673 //--------------
674
675void _mult_ct_t_sin_i(Tbl* tb, int & b)
676{
677 // Peut-etre rien a faire ?
678 if (tb->get_etat() == ETATZERO) {
679 int base_r = b & MSQ_R ;
680 int base_p = b & MSQ_P ;
681 b = base_r | base_p | T_SIN_P ;
682 return ;
683 }
684
685 // Protection
686 assert(tb->get_etat() == ETATQCQ) ;
687
688 // Pour le confort : nbre de points reels.
689 int nr = (tb->dim).dim[0] ;
690 int nt = (tb->dim).dim[1] ;
691 int np = (tb->dim).dim[2] ;
692 np = np - 2 ;
693
694 // zone de travail
695 double* som = new double [nr] ;
696
697 // pt. sur le tableau de double resultat
698 double* xo = new double[(tb->dim).taille] ;
699
700 // Initialisation a zero :
701 for (int i=0; i<(tb->dim).taille; i++) {
702 xo[i] = 0 ;
703 }
704
705 // On y va...
706 double* xi = tb->t ;
707 double* xci = xi ; // Pointeurs
708 double* xco = xo ; // courants
709
710 // k = 0
711
712 // Dernier j: j = nt-1
713 // Positionnement
714 xci += nr * (nt-1) ;
715 xco += nr * (nt-1) ;
716
717 // Initialisation de som
718 for (int i=0 ; i<nr ; i++) {
719 som[i] = 0. ;
720 }
721
722 // j suivants
723 for (int j=nt-1 ; j > 0 ; j--) {
724 // Positionnement
725 xci -= nr ;
726 // Calcul
727 for (int i=0 ; i<nr ; i++ ) {
728 som[i] += 0.5*xci[i] ;
729 xco[i] = som[i] ;
730 som[i] = 0.5*xci[i] ;
731 } // Fin de la boucle sur r
732 xco -= nr ;
733 } // Fin de la boucle sur theta
734 for (int i=0; i<nr; i++) {
735 xco[i] = 0. ;
736 }
737
738 // Positionnement phi suivant
739 xci += nr*nt ;
740 xco += nr*nt ;
741
742 // k = 1
743 xci += nr*nt ;
744 xco += nr*nt ;
745
746 // k >= 2
747 for (int k=2 ; k<np+1 ; k++) {
748
749 // Dernier j: j = nt-1
750 // Positionnement
751 xci += nr * (nt-1) ;
752 xco += nr * (nt-1) ;
753
754 // Initialisation de som
755 for (int i=0 ; i<nr ; i++) {
756 som[i] = 0. ;
757 }
758
759 // j suivants
760 for (int j=nt-1 ; j > 0 ; j--) {
761 // Positionnement
762 xci -= nr ;
763 // Calcul
764 for (int i=0 ; i<nr ; i++ ) {
765 som[i] += 0.5*xci[i] ;
766 xco[i] = som[i] ;
767 som[i] = 0.5*xci[i] ;
768 } // Fin de la boucle sur r
769 xco -= nr ;
770 } // Fin de la boucle sur theta
771 for (int i=0; i<nr; i++) {
772 xco[i] = 0. ;
773 }
774
775 // Positionnement phi suivant
776 xci += nr*nt ;
777 xco += nr*nt ;
778 } // Fin de la boucle sur phi
779
780 // On remet les choses la ou il faut
781 delete [] tb->t ;
782 tb->t = xo ;
783
784 //Menage
785 delete [] som ;
786
787 // base de developpement
788 int base_r = b & MSQ_R ;
789 int base_p = b & MSQ_P ;
790 b = base_r | base_p | T_SIN_P ;
791
792}
793 //---------------------
794 // cas T_COS_I
795 //---------------------
796void _mult_ct_t_cos_i(Tbl* tb, int & b)
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 }
842
843 // j suivants
844 for (int j=nt-1 ; j > 0 ; j--) {
845 // Positionnement
846 xci -= nr ;
847 // Calcul
848 for (int i=0 ; i<nr ; i++ ) {
849 som[i] += 0.5*xci[i] ;
850 xco[i] = som[i] ;
851 som[i] = 0.5*xci[i] ;
852 } // Fin de la boucle sur r
853 xco -= nr ;
854 } // Fin de la boucle sur theta
855 // premier theta
856 for (int i=0 ; i<nr ; i++) {
857 xco[i] = som[i] ;
858 }
859 // Positionnement phi suivant
860 xci += nr*nt ;
861 xco += nr*nt ;
862
863 // k = 1
864 xci += nr*nt ;
865 xco += nr*nt ;
866
867 // k >= 2
868 for (int k=2 ; k<np+1 ; k++) {
869
870 // Dernier j: j = nt-1
871 // Positionnement
872 xci += nr * (nt-1) ;
873 xco += nr * (nt-1) ;
874
875 // Initialisation de som
876 for (int i=0 ; i<nr ; i++) {
877 som[i] = 0. ;
878 }
879
880 // j suivants
881 for (int j=nt-1 ; j > 0 ; j--) {
882 // Positionnement
883 xci -= nr ;
884 // Calcul
885 for (int i=0 ; i<nr ; i++ ) {
886 som[i] += 0.5*xci[i] ;
887 xco[i] = som[i] ;
888 som[i] = 0.5*xci[i] ;
889 } // Fin de la boucle sur r
890 xco -= nr ;
891 } // Fin de la boucle sur theta
892 // premier theta
893 for (int i=0 ; i<nr ; i++) {
894 xco[i] = som[i] ;
895 }
896 // Positionnement phi suivant
897 xci += nr*nt ;
898 xco += nr*nt ;
899 } // Fin de la boucle sur phi
900
901 // On remet les choses la ou il faut
902 delete [] tb->t ;
903 tb->t = xo ;
904
905 //Menage
906 delete [] som ;
907
908 // base de developpement
909 int base_r = b & MSQ_R ;
910 int base_p = b & MSQ_P ;
911 b = base_r | base_p | T_COS_P ;
912
913}
914 //----------------------
915 // cas T_COSSIN_CP
916 //----------------------
917void _mult_ct_t_cossin_cp(Tbl* tb, int & b)
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 // Dernier j: j = nt-1
955 // Positionnement
956 xci += nr * (nt-1) ;
957 xco += nr * (nt-1) ;
958
959 // Initialisation de som
960 for (int i=0 ; i<nr ; i++) {
961 som[i] = 0.5*xci[i] ;
962 xco[i] = 0. ; // mise a zero dui dernier coefficient en theta.
963 }
964
965 // j suivants
966 for (int j=nt-2 ; j > 0 ; j--) {
967 // Positionnement
968 xci -= nr ;
969 xco -= nr ;
970 // Calcul
971 for (int i=0 ; i<nr ; i++ ) {
972 som[i] += 0.5*xci[i] ;
973 xco[i] = som[i] ;
974 som[i] = 0.5*xci[i] ;
975 } // Fin de la boucle sur r
976 } // Fin de la boucle sur theta
977
978 if (nt > 1 ) {
979 // j = 0
980 xci -= nr ;
981 xco -= nr ;
982 for (int i = 0; i<nr; i++) {
983 xco[i] = xci[i] + som[i] ;
984 }
985 }
986 // Positionnement phi suivant
987 xci += nr*nt ;
988 xco += nr*nt ;
989
990 // k=1
991 xci += nr*nt ;
992 xco += nr*nt ;
993
994 // k >= 2
995 for (int k=2 ; k<np+1 ; k++) {
996 m = (k/2) % 2 ; // Parite de l'harmonique en phi
997
998 switch(m) {
999 case 0: // m pair: cos(pair)
1000 // Dernier j: j = nt-1
1001 // Positionnement
1002 xci += nr * (nt-1) ;
1003 xco += nr * (nt-1) ;
1004
1005 // Initialisation de som
1006 for (int i=0 ; i<nr ; i++) {
1007 som[i] = 0.5*xci[i] ;
1008 xco[i] = 0. ; // mise a zero dui dernier coefficient en theta.
1009 }
1010
1011 // j suivants
1012 for (int j=nt-2 ; j > 0 ; j--) {
1013 // Positionnement
1014 xci -= nr ;
1015 xco -= nr ;
1016 // Calcul
1017 for (int i=0 ; i<nr ; i++ ) {
1018 som[i] += 0.5*xci[i] ;
1019 xco[i] = som[i] ;
1020 som[i] = 0.5*xci[i] ;
1021 } // Fin de la boucle sur r
1022 } // Fin de la boucle sur theta
1023 // j = 0
1024 xci -= nr ;
1025 xco -= nr ;
1026 for (int i = 0; i<nr; i++) {
1027 xco[i] = xci[i] + som[i] ;
1028 }
1029 // Positionnement phi suivant
1030 xci += nr*nt ;
1031 xco += nr*nt ;
1032 break ;
1033
1034 case 1: // m impair: sin(impair)
1035 // Dernier j: j = nt-1
1036 // Positionnement
1037 xci += nr * (nt-1) ;
1038 xco += nr * (nt-1) ;
1039
1040 // Initialisation de som
1041 for (int i=0 ; i<nr ; i++) {
1042 som[i] = 0. ;
1043 }
1044
1045 // j suivants
1046 for (int j=nt-1 ; j > 0 ; j--) {
1047 // Positionnement
1048 xci -= nr ;
1049 // Calcul
1050 for (int i=0 ; i<nr ; i++ ) {
1051 som[i] += 0.5*xci[i] ;
1052 xco[i] = som[i] ;
1053 som[i] = 0.5*xci[i] ;
1054 } // Fin de la boucle sur r
1055 xco -= nr ;
1056 } // Fin de la boucle sur theta
1057 for (int i=0; i<nr; i++) {
1058 xco[i] = 0. ;
1059 }
1060
1061 // Positionnement phi suivant
1062 xci += nr*nt ;
1063 xco += nr*nt ;
1064 break;
1065 } // Fin du switch
1066 } // Fin de la boucle sur phi
1067
1068 // On remet les choses la ou il faut
1069 delete [] tb->t ;
1070 tb->t = xo ;
1071
1072 //Menage
1073 delete [] som ;
1074
1075 // base de developpement
1076 int base_r = b & MSQ_R ;
1077 int base_p = b & MSQ_P ;
1078 b = base_r | base_p | T_COSSIN_CI ;
1079}
1080
1081 //------------------
1082 // cas T_COSSIN_CI
1083 //----------------
1084void _mult_ct_t_cossin_ci(Tbl* tb, int & b)
1085{
1086 // Peut-etre rien a faire ?
1087 if (tb->get_etat() == ETATZERO) {
1088 int base_r = b & MSQ_R ;
1089 int base_p = b & MSQ_P ;
1090 b = base_r | base_p | T_COSSIN_CP ;
1091 return ;
1092 }
1093
1094 // Protection
1095 assert(tb->get_etat() == ETATQCQ) ;
1096
1097 // Pour le confort : nbre de points reels.
1098 int nr = (tb->dim).dim[0] ;
1099 int nt = (tb->dim).dim[1] ;
1100 int np = (tb->dim).dim[2] ;
1101 np = np - 2 ;
1102
1103 // zone de travail
1104 double* som = new double [nr] ;
1105
1106 // pt. sur le tableau de double resultat
1107 double* xo = new double[(tb->dim).taille] ;
1108
1109 // Initialisation a zero :
1110 for (int i=0; i<(tb->dim).taille; i++) {
1111 xo[i] = 0 ;
1112 }
1113
1114 // On y va...
1115 double* xi = tb->t ;
1116 double* xci = xi ; // Pointeurs
1117 double* xco = xo ; // courants
1118
1119 // k = 0
1120 int m = 0 ; // Parite de l'harmonique en phi
1121 // Dernier j: j = nt-1
1122 // Positionnement
1123 xci += nr * (nt-1) ;
1124 xco += nr * (nt-1) ;
1125
1126 // Initialisation de som
1127 for (int i=0 ; i<nr ; i++) {
1128 som[i] = 0. ;
1129 }
1130
1131 // j suivants
1132 for (int j=nt-1 ; j > 0 ; j--) {
1133 // Positionnement
1134 xci -= nr ;
1135 // Calcul
1136 for (int i=0 ; i<nr ; i++ ) {
1137 som[i] += 0.5*xci[i] ;
1138 xco[i] = som[i] ;
1139 som[i] = 0.5*xci[i] ;
1140 } // Fin de la boucle sur r
1141 xco -= nr ;
1142 } // Fin de la boucle sur theta
1143 // premier theta
1144 for (int i=0 ; i<nr ; i++) {
1145 xco[i] = som[i] ;
1146 }
1147 // Positionnement phi suivant
1148 xci += nr*nt ;
1149 xco += nr*nt ;
1150
1151 // k=1
1152 xci += nr*nt ;
1153 xco += nr*nt ;
1154
1155 // k >= 2
1156 for (int k=2 ; k<np+1 ; k++) {
1157 m = (k/2) % 2 ; // Parite de l'harmonique en phi
1158
1159 switch(m) {
1160 case 0: // m pair: cos(impair)
1161 // Dernier j: j = nt-1
1162 // Positionnement
1163 xci += nr * (nt-1) ;
1164 xco += nr * (nt-1) ;
1165
1166 // Initialisation de som
1167 for (int i=0 ; i<nr ; i++) {
1168 som[i] = 0. ;
1169 }
1170
1171 // j suivants
1172 for (int j=nt-1 ; j > 0 ; j--) {
1173 // Positionnement
1174 xci -= nr ;
1175 // Calcul
1176 for (int i=0 ; i<nr ; i++ ) {
1177 som[i] += 0.5*xci[i] ;
1178 xco[i] = som[i] ;
1179 som[i] = 0.5*xci[i] ;
1180 } // Fin de la boucle sur r
1181 xco -= nr ;
1182 } // Fin de la boucle sur theta
1183 // premier theta
1184 for (int i=0 ; i<nr ; i++) {
1185 xco[i] = som[i] ;
1186 }
1187 // Positionnement phi suivant
1188 xci += nr*nt ;
1189 xco += nr*nt ;
1190 break ;
1191
1192 case 1:
1193 // Dernier j: j = nt-1
1194 // Positionnement
1195 xci += nr * (nt-1) ;
1196 xco += nr * (nt-1) ;
1197
1198 // Initialisation de som
1199 for (int i=0 ; i<nr ; i++) {
1200 som[i] = 0.5*xci[i] ;
1201 xco[i] = 0. ;
1202 }
1203
1204 // j suivants
1205 for (int j=nt-2 ; j > 0 ; j--) {
1206 // Positionnement
1207 xci -= nr ;
1208 xco -= nr ;
1209 // Calcul
1210 for (int i=0 ; i<nr ; i++ ) {
1211 som[i] += 0.5*xci[i] ;
1212 xco[i] = som[i] ;
1213 som[i] = 0.5*xci[i] ;
1214 } // Fin de la boucle sur r
1215 } // Fin de la boucle sur theta
1216 // j = 0
1217 xci -= nr ;
1218 xco -= nr ;
1219 for (int i=0; i<nr; i++) {
1220 xco[i] = som[i] ;
1221 }
1222 // Positionnement phi suivant
1223 xci += nr*nt ;
1224 xco += nr*nt ;
1225 break ;
1226 } // Fin du switch
1227 } // Fin de la boucle en phi
1228
1229 // On remet les choses la ou il faut
1230 delete [] tb->t ;
1231 tb->t = xo ;
1232
1233 //Menage
1234 delete [] som ;
1235
1236 // base de developpement
1237 int base_r = b & MSQ_R ;
1238 int base_p = b & MSQ_P ;
1239 b = base_r | base_p | T_COSSIN_CP ;
1240}
1241
1242 //---------------------
1243 // cas T_COSSIN_SI
1244 //----------------
1245void _mult_ct_t_cossin_si(Tbl* tb, int & b)
1246{
1247 // Peut-etre rien a faire ?
1248 if (tb->get_etat() == ETATZERO) {
1249 int base_r = b & MSQ_R ;
1250 int base_p = b & MSQ_P ;
1251 b = base_r | base_p | T_COSSIN_SP ;
1252 return ;
1253 }
1254
1255 // Protection
1256 assert(tb->get_etat() == ETATQCQ) ;
1257
1258 // Pour le confort : nbre de points reels.
1259 int nr = (tb->dim).dim[0] ;
1260 int nt = (tb->dim).dim[1] ;
1261 int np = (tb->dim).dim[2] ;
1262 np = np - 2 ;
1263
1264 // zone de travail
1265 double* som = new double [nr] ;
1266
1267 // pt. sur le tableau de double resultat
1268 double* xo = new double[(tb->dim).taille] ;
1269
1270 // Initialisation a zero :
1271 for (int i=0; i<(tb->dim).taille; i++) {
1272 xo[i] = 0 ;
1273 }
1274
1275 // On y va...
1276 double* xi = tb->t ;
1277 double* xci = xi ; // Pointeurs
1278 double* xco = xo ; // courants
1279
1280 // k = 0
1281 int m = 0 ; // Parite de l'harmonique en phi
1282 // Dernier j: j = nt-1
1283 // Positionnement
1284 xci += nr * (nt-1) ;
1285 xco += nr * (nt-1) ;
1286
1287 // Initialisation de som
1288 for (int i=0 ; i<nr ; i++) {
1289 som[i] = 0. ;
1290 }
1291
1292 // j suivants
1293 for (int j=nt-1 ; j > 0 ; j--) {
1294 // Positionnement
1295 xci -= nr ;
1296 // Calcul
1297 for (int i=0 ; i<nr ; i++ ) {
1298 som[i] += 0.5*xci[i] ;
1299 xco[i] = som[i] ;
1300 som[i] = 0.5*xci[i] ;
1301 } // Fin de la boucle sur r
1302 xco -= nr ;
1303 } // Fin de la boucle sur theta
1304 for (int i=0; i<nr; i++) {
1305 xco[i] = 0. ;
1306 }
1307
1308 // Positionnement phi suivant
1309 xci += nr*nt ;
1310 xco += nr*nt ;
1311
1312 // k=1
1313 xci += nr*nt ;
1314 xco += nr*nt ;
1315
1316 // k >= 2
1317 for (int k=2 ; k<np+1 ; k++) {
1318 m = (k/2) % 2 ; // Parite de l'harmonique en phi
1319
1320 switch(m) {
1321 case 0: // m pair: sin(impair)
1322 // Dernier j: j = nt-1
1323 // Positionnement
1324 xci += nr * (nt-1) ;
1325 xco += nr * (nt-1) ;
1326
1327 // Initialisation de som
1328 for (int i=0 ; i<nr ; i++) {
1329 som[i] = 0. ;
1330 }
1331
1332 // j suivants
1333 for (int j=nt-1 ; j > 0 ; j--) {
1334 // Positionnement
1335 xci -= nr ;
1336 // Calcul
1337 for (int i=0 ; i<nr ; i++ ) {
1338 som[i] += 0.5*xci[i] ;
1339 xco[i] = som[i] ;
1340 som[i] = 0.5*xci[i] ;
1341 } // Fin de la boucle sur r
1342 xco -= nr ;
1343 } // Fin de la boucle sur theta
1344 for (int i=0; i<nr; i++) {
1345 xco[i] = 0. ;
1346 }
1347
1348 // Positionnement phi suivant
1349 xci += nr*nt ;
1350 xco += nr*nt ;
1351 break ;
1352
1353 case 1: // m impair cos(pair)
1354 // Dernier j: j = nt-1
1355 // Positionnement
1356 xci += nr * (nt-1) ;
1357 xco += nr * (nt-1) ;
1358
1359 // Initialisation de som
1360 for (int i=0 ; i<nr ; i++) {
1361 som[i] = 0.5*xci[i] ;
1362 xco[i] = 0. ; // mise a zero dui dernier coefficient en theta.
1363 }
1364
1365 // j suivants
1366 for (int j=nt-2 ; j > 0 ; j--) {
1367 // Positionnement
1368 xci -= nr ;
1369 xco -= nr ;
1370 // Calcul
1371 for (int i=0 ; i<nr ; i++ ) {
1372 som[i] += 0.5*xci[i] ;
1373 xco[i] = som[i] ;
1374 som[i] = 0.5*xci[i] ;
1375 } // Fin de la boucle sur r
1376 } // Fin de la boucle sur theta
1377 // j = 0
1378 xci -= nr ;
1379 xco -= nr ;
1380 for (int i = 0; i<nr; i++) {
1381 xco[i] = xci[i] + som[i] ;
1382 }
1383 // Positionnement phi suivant
1384 xci += nr*nt ;
1385 xco += nr*nt ;
1386 break ;
1387 } // Fin du switch
1388 } // Fin de la boucle en phi
1389
1390 // On remet les choses la ou il faut
1391 delete [] tb->t ;
1392 tb->t = xo ;
1393
1394 //Menage
1395 delete [] som ;
1396
1397 // base de developpement
1398 int base_r = b & MSQ_R ;
1399 int base_p = b & MSQ_P ;
1400 b = base_r | base_p | T_COSSIN_SP ;
1401
1402
1403}
1404 //---------------------
1405 // cas T_COSSIN_SP
1406 //----------------
1407void _mult_ct_t_cossin_sp(Tbl* tb, int & b)
1408{
1409 // Peut-etre rien a faire ?
1410 if (tb->get_etat() == ETATZERO) {
1411 int base_r = b & MSQ_R ;
1412 int base_p = b & MSQ_P ;
1413 b = base_r | base_p | T_COSSIN_SI ;
1414 return ;
1415 }
1416
1417 // Protection
1418 assert(tb->get_etat() == ETATQCQ) ;
1419
1420 // Pour le confort : nbre de points reels.
1421 int nr = (tb->dim).dim[0] ;
1422 int nt = (tb->dim).dim[1] ;
1423 int np = (tb->dim).dim[2] ;
1424 np = np - 2 ;
1425
1426 // zone de travail
1427 double* som = new double [nr] ;
1428
1429 // pt. sur le tableau de double resultat
1430 double* xo = new double[(tb->dim).taille] ;
1431
1432 // Initialisation a zero :
1433 for (int i=0; i<(tb->dim).taille; i++) {
1434 xo[i] = 0 ;
1435 }
1436
1437 // On y va...
1438 double* xi = tb->t ;
1439 double* xci = xi ; // Pointeurs
1440 double* xco = xo ; // courants
1441
1442 // k = 0
1443 int m = 0 ; // Parite de l'harmonique en phi
1444 // Dernier j: j = nt-1
1445 // Positionnement
1446 xci += nr * (nt-1) ;
1447 xco += nr * (nt-1) ;
1448
1449 // Initialisation de som
1450 for (int i=0 ; i<nr ; i++) {
1451 som[i] = 0.5*xci[i] ;
1452 xco[i] = 0. ;
1453 }
1454
1455 // j suivants
1456 for (int j=nt-2 ; j > 0 ; j--) {
1457 // Positionnement
1458 xci -= nr ;
1459 xco -= nr ;
1460 // Calcul
1461 for (int i=0 ; i<nr ; i++ ) {
1462 som[i] += 0.5*xci[i] ;
1463 xco[i] = som[i] ;
1464 som[i] = 0.5*xci[i] ;
1465 } // Fin de la boucle sur r
1466 } // Fin de la boucle sur theta
1467
1468 if (nt > 1 ) {
1469 // j = 0
1470 xci -= nr ;
1471 xco -= nr ;
1472 for (int i=0; i<nr; i++) {
1473 xco[i] = som[i] ;
1474 }
1475 }
1476 // Positionnement phi suivant
1477 xci += nr*nt ;
1478 xco += nr*nt ;
1479
1480 // k=1
1481 xci += nr*nt ;
1482 xco += nr*nt ;
1483
1484 for (int k=2 ; k<np+1 ; k++) {
1485 m = (k/2) % 2 ; // Parite de l'harmonique en phi
1486
1487 switch(m) {
1488 case 1: // m impair: cos(impair)
1489 // Dernier j: j = nt-1
1490 // Positionnement
1491 xci += nr * (nt-1) ;
1492 xco += nr * (nt-1) ;
1493
1494 // Initialisation de som
1495 for (int i=0 ; i<nr ; i++) {
1496 som[i] = 0. ;
1497 }
1498
1499 // j suivants
1500 for (int j=nt-1 ; j > 0 ; j--) {
1501 // Positionnement
1502 xci -= nr ;
1503 // Calcul
1504 for (int i=0 ; i<nr ; i++ ) {
1505 som[i] += 0.5*xci[i] ;
1506 xco[i] = som[i] ;
1507 som[i] = 0.5*xci[i] ;
1508 } // Fin de la boucle sur r
1509 xco -= nr ;
1510 } // Fin de la boucle sur theta
1511 // premier theta
1512 for (int i=0 ; i<nr ; i++) {
1513 xco[i] = som[i] ;
1514 }
1515 // Positionnement phi suivant
1516 xci += nr*nt ;
1517 xco += nr*nt ;
1518 break ;
1519
1520 case 0: // m pair: sin(pair)
1521 // Dernier j: j = nt-1
1522 // Positionnement
1523 xci += nr * (nt-1) ;
1524 xco += nr * (nt-1) ;
1525
1526 // Initialisation de som
1527 for (int i=0 ; i<nr ; i++) {
1528 som[i] = 0.5*xci[i] ;
1529 xco[i] = 0. ;
1530 }
1531
1532 // j suivants
1533 for (int j=nt-2 ; j > 0 ; j--) {
1534 // Positionnement
1535 xci -= nr ;
1536 xco -= nr ;
1537 // Calcul
1538 for (int i=0 ; i<nr ; i++ ) {
1539 som[i] += 0.5*xci[i] ;
1540 xco[i] = som[i] ;
1541 som[i] = 0.5*xci[i] ;
1542 } // Fin de la boucle sur r
1543 } // Fin de la boucle sur theta
1544 // j = 0
1545 xci -= nr ;
1546 xco -= nr ;
1547 for (int i=0; i<nr; i++) {
1548 xco[i] = som[i] ;
1549 }
1550 // Positionnement phi suivant
1551 xci += nr*nt ;
1552 xco += nr*nt ;
1553 break ;
1554 } // Fin du switch
1555 } // Fin de la boucle en phi
1556
1557 // On remet les choses la ou il faut
1558 delete [] tb->t ;
1559 tb->t = xo ;
1560
1561 //Menage
1562 delete [] som ;
1563
1564 // base de developpement
1565 int base_r = b & MSQ_R ;
1566 int base_p = b & MSQ_P ;
1567 b = base_r | base_p | T_COSSIN_SI ;
1568
1569}
1570
1571 //----------------------
1572 // cas T_COSSIN_C
1573 //----------------------
1574void _mult_ct_t_cossin_c(Tbl* tb, int & b)
1575{
1576 // Peut-etre rien a faire ?
1577 if (tb->get_etat() == ETATZERO) {
1578 int base_r = b & MSQ_R ;
1579 int base_p = b & MSQ_P ;
1580 switch(base_r){
1581 case(R_CHEBPI_P):
1582 b = R_CHEBPI_I | base_p | T_COSSIN_C ;
1583 break ;
1584 case(R_CHEBPI_I):
1585 b = R_CHEBPI_P | base_p | T_COSSIN_C ;
1586 break ;
1587 default:
1588 b = base_r | base_p | T_COSSIN_C ;
1589 break;
1590 }
1591 return ;
1592 }
1593
1594 // Protection
1595 assert(tb->get_etat() == ETATQCQ) ;
1596
1597 // Pour le confort : nbre de points reels.
1598 int nr = (tb->dim).dim[0] ;
1599 int nt = (tb->dim).dim[1] ;
1600 int np = (tb->dim).dim[2] ;
1601 np = np - 2 ;
1602
1603 // zone de travail
1604 double* som = new double [nr] ;
1605
1606 // pt. sur le tableau de double resultat
1607 double* xo = new double[(tb->dim).taille] ;
1608
1609 // Initialisation a zero :
1610 for (int i=0; i<(tb->dim).taille; i++) {
1611 xo[i] = 0 ;
1612 }
1613
1614 // On y va...
1615 double* xi = tb->t ;
1616 double* xci = xi ; // Pointeurs
1617 double* xco = xo ; // courants
1618
1619 // k = 0
1620 int m = 0 ; // Parite de l'harmonique en phi
1621 // Dernier j: j = nt-1
1622 // Positionnement
1623 xci += nr * (nt-1) ;
1624 xco += nr * (nt-1) ;
1625
1626 // Initialisation de som
1627 for (int i=0 ; i<nr ; i++) {
1628 som[i] = 0.5*xci[i] ;
1629 xco[i] = 0. ; // mise a zero dui dernier coefficient en theta.
1630 }
1631
1632 // j suivants
1633 for (int j=nt-2 ; j > 0 ; j--) {
1634 // Positionnement
1635 xci -= 2*nr ;
1636 xco -= nr ;
1637 // Calcul
1638 for (int i=0 ; i<nr ; i++ ) {
1639 som[i] += 0.5*xci[i] ;
1640 xco[i] = som[i] ;
1641 } // Fin de la boucle sur r
1642 xci += nr ;
1643 for (int i=0 ; i<nr ; i++ ) {
1644 som[i] = 0.5*xci[i] ;
1645 }
1646 } // Fin de la boucle sur theta
1647 // j = 0 : le facteur 2...
1648 xci -= nr ;
1649 for (int i=0; i<nr; i++) {
1650 xco[i] += 0.5*xci[i] ;
1651 }
1652 xco -= nr ;
1653 for (int i = 0; i<nr; i++) {
1654 xco[i] = som[i] ;
1655 }
1656 // Positionnement phi suivant
1657 xci += nr*nt ;
1658 xco += nr*nt ;
1659
1660 // k=1
1661 xci += nr*nt ;
1662 xco += nr*nt ;
1663
1664 // k >= 2
1665 for (int k=2 ; k<np+1 ; k++) {
1666 m = (k/2) % 2 ; // Parite de l'harmonique en phi
1667
1668 switch(m) {
1669 case 0: // m pair: cos
1670 // Dernier j: j = nt-1
1671 // Positionnement
1672 xci += nr * (nt-1) ;
1673 xco += nr * (nt-1) ;
1674
1675 // Initialisation de som
1676 for (int i=0 ; i<nr ; i++) {
1677 som[i] = 0.5*xci[i] ;
1678 xco[i] = 0. ; // mise a zero dui dernier coefficient en theta.
1679 }
1680
1681 // j suivants
1682 for (int j=nt-2 ; j > 0 ; j--) {
1683 // Positionnement
1684 xci -= 2*nr ;
1685 xco -= nr ;
1686 // Calcul
1687 for (int i=0 ; i<nr ; i++ ) {
1688 som[i] += 0.5*xci[i] ;
1689 xco[i] = som[i] ;
1690 } // Fin de la boucle sur r
1691 xci += nr ;
1692 for (int i=0 ; i<nr ; i++ ) {
1693 som[i] = 0.5*xci[i] ;
1694 }
1695 } // Fin de la boucle sur theta
1696 // j = 0 : le facteur 2...
1697 xci -= nr ;
1698 for (int i=0; i<nr; i++) {
1699 xco[i] += 0.5*xci[i] ;
1700 }
1701 xco -= nr ;
1702 for (int i = 0; i<nr; i++) {
1703 xco[i] = som[i] ;
1704 }
1705 // Positionnement phi suivant
1706 xci += nr*nt ;
1707 xco += nr*nt ;
1708 break ;
1709
1710 case 1: // m impair: sin
1711 // Dernier j: j = nt-1
1712 // Positionnement
1713 xci += nr * (nt-1) ;
1714 xco += nr * (nt-1) ;
1715
1716 // Initialisation de som
1717 for (int i=0 ; i<nr ; i++) {
1718 som[i] = 0.5*xci[i] ;
1719 xco[i] = 0.0 ;
1720 }
1721
1722 // j suivants
1723 for (int j=nt-2 ; j > 0 ; j--) {
1724 // Positionnement
1725 xci -= 2*nr ;
1726 xco -= nr ;
1727 // Calcul
1728 for (int i=0 ; i<nr ; i++ ) {
1729 som[i] += 0.5*xci[i] ;
1730 xco[i] = som[i] ;
1731 } // Fin de la boucle sur r
1732 xci += nr ;
1733 for (int i=0 ; i<nr ; i++ ) {
1734 som[i] = 0.5*xci[i] ;
1735 }
1736 } // Fin de la boucle sur theta
1737 // j = 0 : sin(0*theta)
1738 xci -= nr ;
1739 xco -= nr ;
1740 for (int i=0; i<nr; i++) {
1741 xco[i] = 0. ;
1742 }
1743
1744 // Positionnement phi suivant
1745 xci += nr*nt ;
1746 xco += nr*nt ;
1747 break;
1748 } // Fin du switch
1749 } // Fin de la boucle sur phi
1750
1751 // On remet les choses la ou il faut
1752 delete [] tb->t ;
1753 tb->t = xo ;
1754
1755 //Menage
1756 delete [] som ;
1757
1758 // base de developpement
1759 int base_r = b & MSQ_R ;
1760 int base_p = b & MSQ_P ;
1761 switch(base_r){
1762 case(R_CHEBPI_P):
1763 b = R_CHEBPI_I | base_p | T_COSSIN_C ;
1764 break ;
1765 case(R_CHEBPI_I):
1766 b = R_CHEBPI_P | base_p | T_COSSIN_C ;
1767 break ;
1768 default:
1769 b = base_r | base_p | T_COSSIN_C ;
1770 break;
1771 }
1772}
1773
1774 //---------------------
1775 // cas T_COSSIN_S
1776 //----------------
1777void _mult_ct_t_cossin_s(Tbl* tb, int & b)
1778{
1779 // Peut-etre rien a faire ?
1780 if (tb->get_etat() == ETATZERO) {
1781 int base_r = b & MSQ_R ;
1782 int base_p = b & MSQ_P ;
1783 switch(base_r){
1784 case(R_CHEBPI_P):
1785 b = R_CHEBPI_I | base_p | T_COSSIN_S ;
1786 break ;
1787 case(R_CHEBPI_I):
1788 b = R_CHEBPI_P | base_p | T_COSSIN_S ;
1789 break ;
1790 default:
1791 b = base_r | base_p | T_COSSIN_S ;
1792 break;
1793 }
1794 return ;
1795 }
1796
1797 // Protection
1798 assert(tb->get_etat() == ETATQCQ) ;
1799
1800 // Pour le confort : nbre de points reels.
1801 int nr = (tb->dim).dim[0] ;
1802 int nt = (tb->dim).dim[1] ;
1803 int np = (tb->dim).dim[2] ;
1804 np = np - 2 ;
1805
1806 // zone de travail
1807 double* som = new double [nr] ;
1808
1809 // pt. sur le tableau de double resultat
1810 double* xo = new double[(tb->dim).taille] ;
1811
1812 // Initialisation a zero :
1813 for (int i=0; i<(tb->dim).taille; i++) {
1814 xo[i] = 0 ;
1815 }
1816
1817 // On y va...
1818 double* xi = tb->t ;
1819 double* xci = xi ; // Pointeurs
1820 double* xco = xo ; // courants
1821
1822 // k = 0
1823 int m = 0 ; // Parite de l'harmonique en phi
1824 // Dernier j: j = nt-1
1825 // Positionnement
1826 xci += nr * (nt-1) ;
1827 xco += nr * (nt-1) ;
1828
1829 // Initialisation de som
1830 for (int i=0 ; i<nr ; i++) {
1831 som[i] = 0.5*xci[i] ;
1832 xco[i] = 0. ;
1833 }
1834
1835 // j suivants
1836 for (int j=nt-2 ; j > 0 ; j--) {
1837 // Positionnement
1838 xci -= 2*nr ;
1839 xco -= nr ;
1840 // Calcul
1841 for (int i=0 ; i<nr ; i++ ) {
1842 som[i] += 0.5*xci[i] ;
1843 xco[i] = som[i] ;
1844 } // Fin de la boucle sur r
1845 xci += nr ;
1846 for (int i=0 ; i<nr ; i++ ) {
1847 som[i] = 0.5*xci[i] ;
1848 }
1849 } // Fin de la boucle sur theta
1850 // j = 0 : sin(0*theta)
1851 xci -= nr ;
1852 xco -= nr ;
1853 for (int i=0; i<nr; i++) {
1854 xco[i] = 0.0 ;
1855 }
1856 // Positionnement phi suivant
1857 xci += nr*nt ;
1858 xco += nr*nt ;
1859
1860 // k=1
1861 xci += nr*nt ;
1862 xco += nr*nt ;
1863
1864 for (int k=2 ; k<np+1 ; k++) {
1865 m = (k/2) % 2 ; // Parite de l'harmonique en phi
1866
1867 switch(m) {
1868 case 1: // m impair: cos
1869 // Dernier j: j = nt-1
1870 // Positionnement
1871 xci += nr * (nt-1) ;
1872 xco += nr * (nt-1) ;
1873
1874 // Initialisation de som
1875 for (int i=0 ; i<nr ; i++) {
1876 som[i] = 0.5*xci[i] ;
1877 xco[i] = 0.0 ;
1878 }
1879
1880 // j suivants
1881 for (int j=nt-2 ; j > 0 ; j--) {
1882 // Positionnement
1883 xci -= 2*nr ;
1884 xco -= nr ;
1885 // Calcul
1886 for (int i=0 ; i<nr ; i++ ) {
1887 som[i] += 0.5*xci[i] ;
1888 xco[i] = som[i] ;
1889 } // Fin de la boucle sur r
1890 xci += nr ;
1891 for (int i=0 ; i<nr ; i++ ) {
1892 som[i] = 0.5*xci[i] ;
1893 }
1894 } // Fin de la boucle sur theta
1895 // premier theta
1896 // j=0 : le facteur 2...
1897 xci -= nr ;
1898 for (int i=0; i<nr; i++) {
1899 xco[i] += 0.5*xci[i] ;
1900 }
1901 xco -= nr ;
1902 for (int i=0 ; i<nr ; i++) {
1903 xco[i] = som[i] ;
1904 }
1905 // Positionnement phi suivant
1906 xci += nr*nt ;
1907 xco += nr*nt ;
1908 break ;
1909
1910 case 0: // m pair: sin
1911 // Dernier j: j = nt-1
1912 // Positionnement
1913 xci += nr * (nt-1) ;
1914 xco += nr * (nt-1) ;
1915
1916 // Initialisation de som
1917 for (int i=0 ; i<nr ; i++) {
1918 som[i] = 0.5*xci[i] ;
1919 xco[i] = 0. ;
1920 }
1921
1922 // j suivants
1923 for (int j=nt-2 ; j > 0 ; j--) {
1924 // Positionnement
1925 xci -= 2*nr ;
1926 xco -= nr ;
1927 // Calcul
1928 for (int i=0 ; i<nr ; i++ ) {
1929 som[i] += 0.5*xci[i] ;
1930 xco[i] = som[i] ;
1931 } // Fin de la boucle sur r
1932 xci += nr ;
1933 for (int i=0 ; i<nr ; i++ ) {
1934 som[i] = 0.5*xci[i] ;
1935 }
1936 } // Fin de la boucle sur theta
1937 // j = 0 : sin(0*theta)
1938 xci -= nr ;
1939 xco -= nr ;
1940 for (int i=0; i<nr; i++) {
1941 xco[i] = 0.0 ;
1942 }
1943 // Positionnement phi suivant
1944 xci += nr*nt ;
1945 xco += nr*nt ;
1946 break ;
1947 } // Fin du switch
1948 } // Fin de la boucle en phi
1949
1950 // On remet les choses la ou il faut
1951 delete [] tb->t ;
1952 tb->t = xo ;
1953
1954 //Menage
1955 delete [] som ;
1956
1957 // base de developpement
1958 int base_r = b & MSQ_R ;
1959 int base_p = b & MSQ_P ;
1960 switch(base_r){
1961 case(R_CHEBPI_P):
1962 b = R_CHEBPI_I | base_p | T_COSSIN_S ;
1963 break ;
1964 case(R_CHEBPI_I):
1965 b = R_CHEBPI_P | base_p | T_COSSIN_S ;
1966 break ;
1967 default:
1968 b = base_r | base_p | T_COSSIN_S ;
1969 break;
1970 }
1971}
1972}
#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