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