LORENE
op_mult_st.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_st_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_mult_st.C,v 1.8 2014/10/13 08:53:25 j_novak Exp $" ;
24
25/*
26 * $Id: op_mult_st.C,v 1.8 2014/10/13 08:53:25 j_novak Exp $
27 * $Log: op_mult_st.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:28 m_forot
45 * Correct T_COSSIN_S and T_COSSIN_C cases
46 *
47 * Revision 1.2 2004/11/23 15:16:01 m_forot
48 *
49 * Added the bases for the cases without any equatorial symmetry
50 * (T_COSSIN_C, T_COSSIN_S, T_LEG, R_CHEBPI_P, R_CHEBPI_I).
51 *
52 * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
53 * LORENE
54 *
55 * Revision 1.3 2000/02/24 16:42:06 eric
56 * Initialisation a zero du tableau xo avant le calcul.
57 *
58 * Revision 1.2 1999/11/23 16:13:53 novak
59 * *** empty log message ***
60 *
61 * Revision 1.1 1999/11/23 14:31:50 novak
62 * Initial revision
63 *
64 *
65 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_mult_st.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 sin theta
71 * (Utilisation interne)
72 *
73 * void _mult_st_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_st_pas_prevu(Tbl * tb, int& base) {
90 cout << "mult_st 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_st_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_SIN ;
110 break ;
111 case(R_CHEBPI_I):
112 b = R_CHEBPI_P | base_p | T_SIN ;
113 break ;
114 default:
115 b = base_r | base_p | T_SIN ;
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 }
173 } // Fin de la boucle sur theta
174 // j = 0
175 xci -= nr ;
176 for (int i=0 ; i<nr ; i++ ) {
177 xco[i] += 0.5*xci[i] ;
178 }
179 xco -= nr ;
180 for (int i=0; i<nr; i++) {
181 xco[i] = 0 ;
182 }
183 // Positionnement phi suivant
184 xci += nr*nt ;
185 xco += nr*nt ;
186
187 // k = 1
188 xci += nr*nt ;
189 xco += nr*nt ;
190
191 // k >= 2
192 for (int k=2 ; k<np+1 ; k++) {
193
194 // Dernier j: j = nt-1
195 // Positionnement
196 xci += nr * (nt-1) ;
197 xco += nr * (nt-1) ;
198
199 // Initialisation de som
200 for (int i=0 ; i<nr ; i++) {
201 som[i] = -0.5*xci[i] ;
202 xco[i] = 0. ; // mise a zero dui dernier coefficient en theta.
203 }
204
205 // j suivants
206 for (int j=nt-2 ; j > 0 ; j--) {
207 // Positionnement
208 xci -= 2*nr ;
209 xco -= nr ;
210 // Calcul
211 for (int i=0 ; i<nr ; i++ ) {
212 som[i] += 0.5*xci[i] ;
213 xco[i] = som[i] ;
214 } // Fin de la boucle sur r
215 xci += nr ;
216 for (int i=0; i<nr; i++) {
217 som[i] = -0.5*xci[i] ;
218 }
219 } // Fin de la boucle sur theta
220 // j = 0
221 xci -= nr ;
222 for (int i = 0; i<nr; i++) {
223 xco[i] += 0.5*xci[i] ;
224 }
225 xco -= nr ;
226 for (int i=0; i<nr; i++) {
227 xco[i] = 0 ;
228 }
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_SIN ;
247 break ;
248 case(R_CHEBPI_I):
249 b = R_CHEBPI_P | base_p | T_SIN ;
250 break ;
251 default:
252 b = base_r | base_p | T_SIN ;
253 break;
254 }
255}
256
257 //------------
258 // cas T_SIN
259 //------------
260
261void _mult_st_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_COS ;
270 break ;
271 case(R_CHEBPI_I):
272 b = R_CHEBPI_P | base_p | T_COS ;
273 break ;
274 default:
275 b = base_r | base_p | T_COS ;
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
335 xci -= nr ;
336 xco -= nr ;
337 for (int i =0; i<nr ; i++) {
338 xco[i] = som[i] ;
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
378 xci -= nr ;
379 xco -= nr ;
380 for (int i=0; i<nr; i++) {
381 xco[i] = som[i] ;
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_COS ;
401 break ;
402 case(R_CHEBPI_I):
403 b = R_CHEBPI_P | base_p | T_COS ;
404 break ;
405 default:
406 b = base_r | base_p | T_COS ;
407 break;
408 }
409}
410 //----------------
411 // cas T_COS_P
412 //----------------
413
414void _mult_st_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_SIN_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 if (nt > 1) {
476 // j = 0
477 xci -= nr ;
478 xco -= nr ;
479 for (int i=0 ; i<nr ; i++ ) {
480 xco[i] = som[i]+xci[i] ;
481 } // Fin de la boucle sur r
482 }
483 // Positionnement phi suivant
484 xci += nr*nt ;
485 xco += nr*nt ;
486
487 // k = 1
488 xci += nr*nt ;
489 xco += nr*nt ;
490
491 // k >= 2
492 for (int k=2 ; k<np+1 ; k++) {
493
494 // Dernier j: j = nt-1
495 // Positionnement
496 xci += nr * (nt-1) ;
497 xco += nr * (nt-1) ;
498
499 // Initialisation de som
500 for (int i=0 ; i<nr ; i++) {
501 som[i] = -0.5*xci[i] ;
502 xco[i] = 0. ; // mise a zero dui dernier coefficient en theta.
503 }
504
505 // j suivants
506 for (int j=nt-2 ; j > 0 ; j--) {
507 // Positionnement
508 xci -= nr ;
509 xco -= nr ;
510 // Calcul
511 for (int i=0 ; i<nr ; i++ ) {
512 som[i] += 0.5*xci[i] ;
513 xco[i] = som[i] ;
514 som[i] = -0.5*xci[i] ;
515 } // Fin de la boucle sur r
516 } // Fin de la boucle sur theta
517 // j = 0
518 xci -= nr ;
519 xco -= nr ;
520 for (int i = 0; i<nr; i++) {
521 xco[i] = xci[i] + som[i] ;
522 }
523 // Positionnement phi suivant
524 xci += nr*nt ;
525 xco += nr*nt ;
526 } // Fin de la boucle sur phi
527
528 // On remet les choses la ou il faut
529 delete [] tb->t ;
530 tb->t = xo ;
531
532 //Menage
533 delete [] som ;
534
535 // base de developpement
536 int base_r = b & MSQ_R ;
537 int base_p = b & MSQ_P ;
538 b = base_r | base_p | T_SIN_I ;
539}
540
541 //--------------
542 // cas T_SIN_P
543 //--------------
544
545void _mult_st_t_sin_p(Tbl* tb, int & b)
546{
547 // Peut-etre rien a faire ?
548 if (tb->get_etat() == ETATZERO) {
549 int base_r = b & MSQ_R ;
550 int base_p = b & MSQ_P ;
551 b = base_r | base_p | T_COS_I ;
552 return ;
553 }
554
555 // Protection
556 assert(tb->get_etat() == ETATQCQ) ;
557
558 // Pour le confort : nbre de points reels.
559 int nr = (tb->dim).dim[0] ;
560 int nt = (tb->dim).dim[1] ;
561 int np = (tb->dim).dim[2] ;
562 np = np - 2 ;
563
564 // zone de travail
565 double* som = new double [nr] ;
566
567 // pt. sur le tableau de double resultat
568 double* xo = new double[(tb->dim).taille] ;
569
570 // Initialisation a zero :
571 for (int i=0; i<(tb->dim).taille; i++) {
572 xo[i] = 0 ;
573 }
574
575 // On y va...
576 double* xi = tb->t ;
577 double* xci = xi ; // Pointeurs
578 double* xco = xo ; // courants
579
580 // k = 0
581
582 // Dernier j: j = nt-1
583 // Positionnement
584 xci += nr * (nt-1) ;
585 xco += nr * (nt-1) ;
586
587 // Initialisation de som
588 for (int i=0 ; i<nr ; i++) {
589 som[i] = 0.5*xci[i] ;
590 xco[i] = 0. ;
591 }
592
593 // j suivants
594 for (int j=nt-2 ; j > 0 ; j--) {
595 // Positionnement
596 xci -= nr ;
597 xco -= nr ;
598 // Calcul
599 for (int i=0 ; i<nr ; i++ ) {
600 som[i] -= 0.5*xci[i] ;
601 xco[i] = som[i] ;
602 som[i] = 0.5*xci[i] ;
603 } // Fin de la boucle sur r
604 } // Fin de la boucle sur theta
605 if (nt > 1) {
606 // j = 0
607 xci -= nr ;
608 xco -= nr ;
609 for (int i =0; i<nr ; i++) {
610 xco[i] = som[i] ;
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_COS_I ;
669
670}
671 //--------------
672 // cas T_SIN_I
673 //--------------
674
675void _mult_st_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_COS_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 // premier theta
735 for (int i=0 ; i<nr ; i++) {
736 xco[i] = som[i] ;
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 // premier theta
772 for (int i=0 ; i<nr ; i++) {
773 xco[i] = som[i] ;
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_COS_P ;
791
792}
793 //---------------------
794 // cas T_COS_I
795 //---------------------
796void _mult_st_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_SIN_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 for (int i=0; i<nr; i++) {
856 xco[i] = 0. ;
857 }
858 // Positionnement phi suivant
859 xci += nr*nt ;
860 xco += nr*nt ;
861
862 // k = 1
863 xci += nr*nt ;
864 xco += nr*nt ;
865
866 // k >= 2
867 for (int k=2 ; k<np+1 ; k++) {
868
869 // Dernier j: j = nt-1
870 // Positionnement
871 xci += nr * (nt-1) ;
872 xco += nr * (nt-1) ;
873
874 // Initialisation de som
875 for (int i=0 ; i<nr ; i++) {
876 som[i] = 0. ;
877 }
878
879 // j suivants
880 for (int j=nt-1 ; j > 0 ; j--) {
881 // Positionnement
882 xci -= nr ;
883 // Calcul
884 for (int i=0 ; i<nr ; i++ ) {
885 som[i] += 0.5*xci[i] ;
886 xco[i] = som[i] ;
887 som[i] = -0.5*xci[i] ;
888 } // Fin de la boucle sur r
889 xco -= nr ;
890 } // Fin de la boucle sur theta
891 for (int i=0; i<nr; i++) {
892 xco[i] = 0. ;
893 }
894
895 // Positionnement phi suivant
896 xci += nr*nt ;
897 xco += nr*nt ;
898 } // Fin de la boucle sur phi
899
900 // On remet les choses la ou il faut
901 delete [] tb->t ;
902 tb->t = xo ;
903
904 //Menage
905 delete [] som ;
906
907 // base de developpement
908 int base_r = b & MSQ_R ;
909 int base_p = b & MSQ_P ;
910 b = base_r | base_p | T_SIN_P ;
911
912}
913 //----------------------
914 // cas T_COSSIN_CP
915 //----------------------
916void _mult_st_t_cossin_cp(Tbl* tb, int & b)
917{
918 // Peut-etre rien a faire ?
919 if (tb->get_etat() == ETATZERO) {
920 int base_r = b & MSQ_R ;
921 int base_p = b & MSQ_P ;
922 b = base_r | base_p | T_COSSIN_SI ;
923 return ;
924 }
925
926 // Protection
927 assert(tb->get_etat() == ETATQCQ) ;
928
929 // Pour le confort : nbre de points reels.
930 int nr = (tb->dim).dim[0] ;
931 int nt = (tb->dim).dim[1] ;
932 int np = (tb->dim).dim[2] ;
933 np = np - 2 ;
934
935 // zone de travail
936 double* som = new double [nr] ;
937
938 // pt. sur le tableau de double resultat
939 double* xo = new double[(tb->dim).taille] ;
940
941 // Initialisation a zero :
942 for (int i=0; i<(tb->dim).taille; i++) {
943 xo[i] = 0 ;
944 }
945
946 // On y va...
947 double* xi = tb->t ;
948 double* xci = xi ; // Pointeurs
949 double* xco = xo ; // courants
950
951 // k = 0
952 int m = 0 ; // Parite de l'harmonique en phi
953 // Dernier j: j = nt-1
954 // Positionnement
955 xci += nr * (nt-1) ;
956 xco += nr * (nt-1) ;
957
958 // Initialisation de som
959 for (int i=0 ; i<nr ; i++) {
960 som[i] = -0.5*xci[i] ;
961 xco[i] = 0. ; // mise a zero dui dernier coefficient en theta.
962 }
963
964 // j suivants
965 for (int j=nt-2 ; j > 0 ; j--) {
966 // Positionnement
967 xci -= nr ;
968 xco -= nr ;
969 // Calcul
970 for (int i=0 ; i<nr ; i++ ) {
971 som[i] += 0.5*xci[i] ;
972 xco[i] = som[i] ;
973 som[i] = -0.5*xci[i] ;
974 } // Fin de la boucle sur r
975 } // Fin de la boucle sur theta
976
977 if (nt > 1 ) {
978 // j = 0
979 xci -= nr ;
980 xco -= nr ;
981 for (int i = 0; i<nr; i++) {
982 xco[i] = xci[i] + som[i] ;
983 }
984 }
985 // Positionnement phi suivant
986 xci += nr*nt ;
987 xco += nr*nt ;
988
989 // k=1
990 xci += nr*nt ;
991 xco += nr*nt ;
992
993 // k >= 2
994 for (int k=2 ; k<np+1 ; k++) {
995 m = (k/2) % 2 ; // Parite de l'harmonique en phi
996
997 switch(m) {
998 case 0: // m pair: cos(pair)
999 // Dernier j: j = nt-1
1000 // Positionnement
1001 xci += nr * (nt-1) ;
1002 xco += nr * (nt-1) ;
1003
1004 // Initialisation de som
1005 for (int i=0 ; i<nr ; i++) {
1006 som[i] = -0.5*xci[i] ;
1007 xco[i] = 0. ; // mise a zero dui dernier coefficient en theta.
1008 }
1009
1010 // j suivants
1011 for (int j=nt-2 ; j > 0 ; j--) {
1012 // Positionnement
1013 xci -= nr ;
1014 xco -= nr ;
1015 // Calcul
1016 for (int i=0 ; i<nr ; i++ ) {
1017 som[i] += 0.5*xci[i] ;
1018 xco[i] = som[i] ;
1019 som[i] = -0.5*xci[i] ;
1020 } // Fin de la boucle sur r
1021 } // Fin de la boucle sur theta
1022 // j = 0
1023 xci -= nr ;
1024 xco -= nr ;
1025 for (int i = 0; i<nr; i++) {
1026 xco[i] = xci[i] + som[i] ;
1027 }
1028 // Positionnement phi suivant
1029 xci += nr*nt ;
1030 xco += nr*nt ;
1031 break ;
1032
1033 case 1: // m impair: sin(impair)
1034 // Dernier j: j = nt-1
1035 // Positionnement
1036 xci += nr * (nt-1) ;
1037 xco += nr * (nt-1) ;
1038
1039 // Initialisation de som
1040 for (int i=0 ; i<nr ; i++) {
1041 som[i] = 0. ;
1042 }
1043
1044 // j suivants
1045 for (int j=nt-1 ; j > 0 ; j--) {
1046 // Positionnement
1047 xci -= nr ;
1048 // Calcul
1049 for (int i=0 ; i<nr ; i++ ) {
1050 som[i] -= 0.5*xci[i] ;
1051 xco[i] = som[i] ;
1052 som[i] = 0.5*xci[i] ;
1053 } // Fin de la boucle sur r
1054 xco -= nr ;
1055 } // Fin de la boucle sur theta
1056 // premier theta
1057 for (int i=0 ; i<nr ; i++) {
1058 xco[i] = som[i] ;
1059 }
1060 // Positionnement phi suivant
1061 xci += nr*nt ;
1062 xco += nr*nt ;
1063 break;
1064 } // Fin du switch
1065 } // Fin de la boucle sur phi
1066
1067 // On remet les choses la ou il faut
1068 delete [] tb->t ;
1069 tb->t = xo ;
1070
1071 //Menage
1072 delete [] som ;
1073
1074 // base de developpement
1075 int base_r = b & MSQ_R ;
1076 int base_p = b & MSQ_P ;
1077 b = base_r | base_p | T_COSSIN_SI ;
1078}
1079
1080 //------------------
1081 // cas T_COSSIN_CI
1082 //----------------
1083void _mult_st_t_cossin_ci(Tbl* tb, int & b)
1084{
1085 // Peut-etre rien a faire ?
1086 if (tb->get_etat() == ETATZERO) {
1087 int base_r = b & MSQ_R ;
1088 int base_p = b & MSQ_P ;
1089 b = base_r | base_p | T_COSSIN_SP ;
1090 return ;
1091 }
1092
1093 // Protection
1094 assert(tb->get_etat() == ETATQCQ) ;
1095
1096 // Pour le confort : nbre de points reels.
1097 int nr = (tb->dim).dim[0] ;
1098 int nt = (tb->dim).dim[1] ;
1099 int np = (tb->dim).dim[2] ;
1100 np = np - 2 ;
1101
1102 // zone de travail
1103 double* som = new double [nr] ;
1104
1105 // pt. sur le tableau de double resultat
1106 double* xo = new double[(tb->dim).taille] ;
1107
1108 // Initialisation a zero :
1109 for (int i=0; i<(tb->dim).taille; i++) {
1110 xo[i] = 0 ;
1111 }
1112
1113 // On y va...
1114 double* xi = tb->t ;
1115 double* xci = xi ; // Pointeurs
1116 double* xco = xo ; // courants
1117
1118 // k = 0
1119 int m = 0 ; // Parite de l'harmonique en phi
1120 // Dernier j: j = nt-1
1121 // Positionnement
1122 xci += nr * (nt-1) ;
1123 xco += nr * (nt-1) ;
1124
1125 // Initialisation de som
1126 for (int i=0 ; i<nr ; i++) {
1127 som[i] = 0. ;
1128 }
1129
1130 // j suivants
1131 for (int j=nt-1 ; j > 0 ; j--) {
1132 // Positionnement
1133 xci -= nr ;
1134 // Calcul
1135 for (int i=0 ; i<nr ; i++ ) {
1136 som[i] += 0.5*xci[i] ;
1137 xco[i] = som[i] ;
1138 som[i] = -0.5*xci[i] ;
1139 } // Fin de la boucle sur r
1140 xco -= nr ;
1141 } // Fin de la boucle sur theta
1142 for (int i=0; i<nr; i++) {
1143 xco[i] = 0. ;
1144 }
1145
1146 // Positionnement phi suivant
1147 xci += nr*nt ;
1148 xco += nr*nt ;
1149
1150 // k=1
1151 xci += nr*nt ;
1152 xco += nr*nt ;
1153
1154 // k >= 2
1155 for (int k=2 ; k<np+1 ; k++) {
1156 m = (k/2) % 2 ; // Parite de l'harmonique en phi
1157
1158 switch(m) {
1159 case 0: // m pair: cos(impair)
1160 // Dernier j: j = nt-1
1161 // Positionnement
1162 xci += nr * (nt-1) ;
1163 xco += nr * (nt-1) ;
1164
1165 // Initialisation de som
1166 for (int i=0 ; i<nr ; i++) {
1167 som[i] = 0. ;
1168 }
1169
1170 // j suivants
1171 for (int j=nt-1 ; j > 0 ; j--) {
1172 // Positionnement
1173 xci -= nr ;
1174 // Calcul
1175 for (int i=0 ; i<nr ; i++ ) {
1176 som[i] += 0.5*xci[i] ;
1177 xco[i] = som[i] ;
1178 som[i] = -0.5*xci[i] ;
1179 } // Fin de la boucle sur r
1180 xco -= nr ;
1181 } // Fin de la boucle sur theta
1182 for (int i=0; i<nr; i++) {
1183 xco[i] = 0. ;
1184 }
1185
1186 // Positionnement phi suivant
1187 xci += nr*nt ;
1188 xco += nr*nt ;
1189 break ;
1190
1191 case 1:
1192 // Dernier j: j = nt-1
1193 // Positionnement
1194 xci += nr * (nt-1) ;
1195 xco += nr * (nt-1) ;
1196
1197 // Initialisation de som
1198 for (int i=0 ; i<nr ; i++) {
1199 som[i] = 0.5*xci[i] ;
1200 xco[i] = 0. ;
1201 }
1202
1203 // j suivants
1204 for (int j=nt-2 ; j > 0 ; j--) {
1205 // Positionnement
1206 xci -= nr ;
1207 xco -= nr ;
1208 // Calcul
1209 for (int i=0 ; i<nr ; i++ ) {
1210 som[i] -= 0.5*xci[i] ;
1211 xco[i] = som[i] ;
1212 som[i] = 0.5*xci[i] ;
1213 } // Fin de la boucle sur r
1214 } // Fin de la boucle sur theta
1215 // j = 0
1216 xci -= nr ;
1217 xco -= nr ;
1218 for (int i=0; i<nr; i++) {
1219 xco[i] = som[i] ;
1220 }
1221 // Positionnement phi suivant
1222 xci += nr*nt ;
1223 xco += nr*nt ;
1224 break ;
1225 } // Fin du switch
1226 } // Fin de la boucle en phi
1227
1228 // On remet les choses la ou il faut
1229 delete [] tb->t ;
1230 tb->t = xo ;
1231
1232 //Menage
1233 delete [] som ;
1234
1235 // base de developpement
1236 int base_r = b & MSQ_R ;
1237 int base_p = b & MSQ_P ;
1238 b = base_r | base_p | T_COSSIN_SP ;
1239}
1240
1241 //---------------------
1242 // cas T_COSSIN_SI
1243 //----------------
1244void _mult_st_t_cossin_si(Tbl* tb, int & b)
1245{
1246 // Peut-etre rien a faire ?
1247 if (tb->get_etat() == ETATZERO) {
1248 int base_r = b & MSQ_R ;
1249 int base_p = b & MSQ_P ;
1250 b = base_r | base_p | T_COSSIN_CP ;
1251 return ;
1252 }
1253
1254 // Protection
1255 assert(tb->get_etat() == ETATQCQ) ;
1256
1257 // Pour le confort : nbre de points reels.
1258 int nr = (tb->dim).dim[0] ;
1259 int nt = (tb->dim).dim[1] ;
1260 int np = (tb->dim).dim[2] ;
1261 np = np - 2 ;
1262
1263 // zone de travail
1264 double* som = new double [nr] ;
1265
1266 // pt. sur le tableau de double resultat
1267 double* xo = new double[(tb->dim).taille] ;
1268
1269 // Initialisation a zero :
1270 for (int i=0; i<(tb->dim).taille; i++) {
1271 xo[i] = 0 ;
1272 }
1273
1274 // On y va...
1275 double* xi = tb->t ;
1276 double* xci = xi ; // Pointeurs
1277 double* xco = xo ; // courants
1278
1279 // k = 0
1280 int m = 0 ; // Parite de l'harmonique en phi
1281 // Dernier j: j = nt-1
1282 // Positionnement
1283 xci += nr * (nt-1) ;
1284 xco += nr * (nt-1) ;
1285
1286 // Initialisation de som
1287 for (int i=0 ; i<nr ; i++) {
1288 som[i] = 0. ;
1289 }
1290
1291 // j suivants
1292 for (int j=nt-1 ; j > 0 ; j--) {
1293 // Positionnement
1294 xci -= nr ;
1295 // Calcul
1296 for (int i=0 ; i<nr ; i++ ) {
1297 som[i] -= 0.5*xci[i] ;
1298 xco[i] = som[i] ;
1299 som[i] = 0.5*xci[i] ;
1300 } // Fin de la boucle sur r
1301 xco -= nr ;
1302 } // Fin de la boucle sur theta
1303 // premier theta
1304 for (int i=0 ; i<nr ; i++) {
1305 xco[i] = som[i] ;
1306 }
1307 // Positionnement phi suivant
1308 xci += nr*nt ;
1309 xco += nr*nt ;
1310
1311 // k=1
1312 xci += nr*nt ;
1313 xco += nr*nt ;
1314
1315 // k >= 2
1316 for (int k=2 ; k<np+1 ; k++) {
1317 m = (k/2) % 2 ; // Parite de l'harmonique en phi
1318
1319 switch(m) {
1320 case 0: // m pair: sin(impair)
1321 // Dernier j: j = nt-1
1322 // Positionnement
1323 xci += nr * (nt-1) ;
1324 xco += nr * (nt-1) ;
1325
1326 // Initialisation de som
1327 for (int i=0 ; i<nr ; i++) {
1328 som[i] = 0. ;
1329 }
1330
1331 // j suivants
1332 for (int j=nt-1 ; j > 0 ; j--) {
1333 // Positionnement
1334 xci -= nr ;
1335 // Calcul
1336 for (int i=0 ; i<nr ; i++ ) {
1337 som[i] -= 0.5*xci[i] ;
1338 xco[i] = som[i] ;
1339 som[i] = 0.5*xci[i] ;
1340 } // Fin de la boucle sur r
1341 xco -= nr ;
1342 } // Fin de la boucle sur theta
1343 // premier theta
1344 for (int i=0 ; i<nr ; i++) {
1345 xco[i] = som[i] ;
1346 }
1347 // Positionnement phi suivant
1348 xci += nr*nt ;
1349 xco += nr*nt ;
1350 break ;
1351
1352 case 1: // m impair cos(pair)
1353 // Dernier j: j = nt-1
1354 // Positionnement
1355 xci += nr * (nt-1) ;
1356 xco += nr * (nt-1) ;
1357
1358 // Initialisation de som
1359 for (int i=0 ; i<nr ; i++) {
1360 som[i] = -0.5*xci[i] ;
1361 xco[i] = 0. ; // mise a zero dui dernier coefficient en theta.
1362 }
1363
1364 // j suivants
1365 for (int j=nt-2 ; j > 0 ; j--) {
1366 // Positionnement
1367 xci -= nr ;
1368 xco -= nr ;
1369 // Calcul
1370 for (int i=0 ; i<nr ; i++ ) {
1371 som[i] += 0.5*xci[i] ;
1372 xco[i] = som[i] ;
1373 som[i] = -0.5*xci[i] ;
1374 } // Fin de la boucle sur r
1375 } // Fin de la boucle sur theta
1376 // j = 0
1377 xci -= nr ;
1378 xco -= nr ;
1379 for (int i = 0; i<nr; i++) {
1380 xco[i] = xci[i] + som[i] ;
1381 }
1382 // Positionnement phi suivant
1383 xci += nr*nt ;
1384 xco += nr*nt ;
1385 break ;
1386 } // Fin du switch
1387 } // Fin de la boucle en phi
1388
1389 // On remet les choses la ou il faut
1390 delete [] tb->t ;
1391 tb->t = xo ;
1392
1393 //Menage
1394 delete [] som ;
1395
1396 // base de developpement
1397 int base_r = b & MSQ_R ;
1398 int base_p = b & MSQ_P ;
1399 b = base_r | base_p | T_COSSIN_CP ;
1400
1401
1402}
1403 //---------------------
1404 // cas T_COSSIN_SP
1405 //----------------
1406void _mult_st_t_cossin_sp(Tbl* tb, int & b)
1407{
1408 // Peut-etre rien a faire ?
1409 if (tb->get_etat() == ETATZERO) {
1410 int base_r = b & MSQ_R ;
1411 int base_p = b & MSQ_P ;
1412 b = base_r | base_p | T_COSSIN_CI ;
1413 return ;
1414 }
1415
1416 // Protection
1417 assert(tb->get_etat() == ETATQCQ) ;
1418
1419 // Pour le confort : nbre de points reels.
1420 int nr = (tb->dim).dim[0] ;
1421 int nt = (tb->dim).dim[1] ;
1422 int np = (tb->dim).dim[2] ;
1423 np = np - 2 ;
1424
1425 // zone de travail
1426 double* som = new double [nr] ;
1427
1428 // pt. sur le tableau de double resultat
1429 double* xo = new double[(tb->dim).taille] ;
1430
1431 // Initialisation a zero :
1432 for (int i=0; i<(tb->dim).taille; i++) {
1433 xo[i] = 0 ;
1434 }
1435
1436 // On y va...
1437 double* xi = tb->t ;
1438 double* xci = xi ; // Pointeurs
1439 double* xco = xo ; // courants
1440
1441 // k = 0
1442 int m = 0 ; // Parite de l'harmonique en phi
1443 // Dernier j: j = nt-1
1444 // Positionnement
1445 xci += nr * (nt-1) ;
1446 xco += nr * (nt-1) ;
1447
1448 // Initialisation de som
1449 for (int i=0 ; i<nr ; i++) {
1450 som[i] = 0.5*xci[i] ;
1451 xco[i] = 0. ;
1452 }
1453
1454 // j suivants
1455 for (int j=nt-2 ; j > 0 ; j--) {
1456 // Positionnement
1457 xci -= nr ;
1458 xco -= nr ;
1459 // Calcul
1460 for (int i=0 ; i<nr ; i++ ) {
1461 som[i] -= 0.5*xci[i] ;
1462 xco[i] = som[i] ;
1463 som[i] = 0.5*xci[i] ;
1464 } // Fin de la boucle sur r
1465 } // Fin de la boucle sur theta
1466
1467 if (nt > 1 ) {
1468 // j = 0
1469 xci -= nr ;
1470 xco -= nr ;
1471 for (int i=0; i<nr; i++) {
1472 xco[i] = som[i] ;
1473 }
1474 }
1475 // Positionnement phi suivant
1476 xci += nr*nt ;
1477 xco += nr*nt ;
1478
1479 // k=1
1480 xci += nr*nt ;
1481 xco += nr*nt ;
1482
1483 for (int k=2 ; k<np+1 ; k++) {
1484 m = (k/2) % 2 ; // Parite de l'harmonique en phi
1485
1486 switch(m) {
1487 case 1: // m impair: cos(impair)
1488 // Dernier j: j = nt-1
1489 // Positionnement
1490 xci += nr * (nt-1) ;
1491 xco += nr * (nt-1) ;
1492
1493 // Initialisation de som
1494 for (int i=0 ; i<nr ; i++) {
1495 som[i] = 0. ;
1496 }
1497
1498 // j suivants
1499 for (int j=nt-1 ; j > 0 ; j--) {
1500 // Positionnement
1501 xci -= nr ;
1502 // Calcul
1503 for (int i=0 ; i<nr ; i++ ) {
1504 som[i] += 0.5*xci[i] ;
1505 xco[i] = som[i] ;
1506 som[i] = -0.5*xci[i] ;
1507 } // Fin de la boucle sur r
1508 xco -= nr ;
1509 } // Fin de la boucle sur theta
1510 for (int i=0; i<nr; i++) {
1511 xco[i] = 0. ;
1512 }
1513
1514 // Positionnement phi suivant
1515 xci += nr*nt ;
1516 xco += nr*nt ;
1517 break ;
1518
1519 case 0: // m pair: sin(pair)
1520 // Dernier j: j = nt-1
1521 // Positionnement
1522 xci += nr * (nt-1) ;
1523 xco += nr * (nt-1) ;
1524
1525 // Initialisation de som
1526 for (int i=0 ; i<nr ; i++) {
1527 som[i] = 0.5*xci[i] ;
1528 xco[i] = 0. ;
1529 }
1530
1531 // j suivants
1532 for (int j=nt-2 ; j > 0 ; j--) {
1533 // Positionnement
1534 xci -= nr ;
1535 xco -= nr ;
1536 // Calcul
1537 for (int i=0 ; i<nr ; i++ ) {
1538 som[i] -= 0.5*xci[i] ;
1539 xco[i] = som[i] ;
1540 som[i] = 0.5*xci[i] ;
1541 } // Fin de la boucle sur r
1542 } // Fin de la boucle sur theta
1543 // j = 0
1544 xci -= nr ;
1545 xco -= nr ;
1546 for (int i=0; i<nr; i++) {
1547 xco[i] = som[i] ;
1548 }
1549 // Positionnement phi suivant
1550 xci += nr*nt ;
1551 xco += nr*nt ;
1552 break ;
1553 } // Fin du switch
1554 } // Fin de la boucle en phi
1555
1556 // On remet les choses la ou il faut
1557 delete [] tb->t ;
1558 tb->t = xo ;
1559
1560 //Menage
1561 delete [] som ;
1562
1563 // base de developpement
1564 int base_r = b & MSQ_R ;
1565 int base_p = b & MSQ_P ;
1566 b = base_r | base_p | T_COSSIN_CI ;
1567
1568}
1569
1570 //----------------------
1571 // cas T_COSSIN_C
1572 //----------------------
1573void _mult_st_t_cossin_c(Tbl* tb, int & b)
1574{
1575 // Peut-etre rien a faire ?
1576 if (tb->get_etat() == ETATZERO) {
1577 int base_r = b & MSQ_R ;
1578 int base_p = b & MSQ_P ;
1579 switch(base_r){
1580 case(R_CHEBPI_P):
1581 b = R_CHEBPI_I | base_p | T_COSSIN_S ;
1582 break ;
1583 case(R_CHEBPI_I):
1584 b = R_CHEBPI_P | base_p | T_COSSIN_S ;
1585 break ;
1586 default:
1587 b = base_r | base_p | T_COSSIN_S ;
1588 break;
1589 }
1590 return ;
1591 }
1592
1593 // Protection
1594 assert(tb->get_etat() == ETATQCQ) ;
1595
1596 // Pour le confort : nbre de points reels.
1597 int nr = (tb->dim).dim[0] ;
1598 int nt = (tb->dim).dim[1] ;
1599 int np = (tb->dim).dim[2] ;
1600 np = np - 2 ;
1601
1602 // zone de travail
1603 double* som = new double [nr] ;
1604
1605 // pt. sur le tableau de double resultat
1606 double* xo = new double[(tb->dim).taille] ;
1607
1608 // Initialisation a zero :
1609 for (int i=0; i<(tb->dim).taille; i++) {
1610 xo[i] = 0 ;
1611 }
1612
1613 // On y va...
1614 double* xi = tb->t ;
1615 double* xci = xi ; // Pointeurs
1616 double* xco = xo ; // courants
1617
1618 // k = 0
1619 int m = 0 ; // Parite de l'harmonique en phi
1620 // Dernier j: j = nt-1
1621 // Positionnement
1622 xci += nr * (nt-1) ;
1623 xco += nr * (nt-1) ;
1624
1625 // Initialisation de som
1626 for (int i=0 ; i<nr ; i++) {
1627 som[i] = -0.5*xci[i] ;
1628 xco[i] = 0. ; // mise a zero dui dernier coefficient en theta.
1629 }
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
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] = 0. ;
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
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] = 0. ;
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. ; // mise a zero dui dernier coefficient en theta.
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 xci -= nr;
1738 xco -= nr;
1739 // premier theta
1740 for (int i=0 ; i<nr ; i++) {
1741 xco[i] = som[i] ;
1742 }
1743 // Positionnement phi suivant
1744 xci += nr*nt ;
1745 xco += nr*nt ;
1746 break;
1747 } // Fin du switch
1748 } // Fin de la boucle sur phi
1749
1750 // On remet les choses la ou il faut
1751 delete [] tb->t ;
1752 tb->t = xo ;
1753
1754 //Menage
1755 delete [] som ;
1756
1757 // base de developpement
1758 int base_r = b & MSQ_R ;
1759 int base_p = b & MSQ_P ;
1760 switch(base_r){
1761 case(R_CHEBPI_P):
1762 b = R_CHEBPI_I | base_p | T_COSSIN_S ;
1763 break ;
1764 case(R_CHEBPI_I):
1765 b = R_CHEBPI_P | base_p | T_COSSIN_S ;
1766 break ;
1767 default:
1768 b = base_r | base_p | T_COSSIN_S ;
1769 break;
1770 }
1771}
1772
1773 //---------------------
1774 // cas T_COSSIN_S
1775 //----------------
1776void _mult_st_t_cossin_s(Tbl* tb, int & b)
1777{
1778 // Peut-etre rien a faire ?
1779 if (tb->get_etat() == ETATZERO) {
1780 int base_r = b & MSQ_R ;
1781 int base_p = b & MSQ_P ;
1782 switch(base_r){
1783 case(R_CHEBPI_P):
1784 b = R_CHEBPI_I | base_p | T_COSSIN_C ;
1785 break ;
1786 case(R_CHEBPI_I):
1787 b = R_CHEBPI_P | base_p | T_COSSIN_C ;
1788 break ;
1789 default:
1790 b = base_r | base_p | T_COSSIN_C ;
1791 break;
1792 }
1793 return ;
1794 }
1795
1796 // Protection
1797 assert(tb->get_etat() == ETATQCQ) ;
1798
1799 // Pour le confort : nbre de points reels.
1800 int nr = (tb->dim).dim[0] ;
1801 int nt = (tb->dim).dim[1] ;
1802 int np = (tb->dim).dim[2] ;
1803 np = np - 2 ;
1804
1805 // zone de travail
1806 double* som = new double [nr] ;
1807
1808 // pt. sur le tableau de double resultat
1809 double* xo = new double[(tb->dim).taille] ;
1810
1811 // Initialisation a zero :
1812 for (int i=0; i<(tb->dim).taille; i++) {
1813 xo[i] = 0 ;
1814 }
1815
1816 // On y va...
1817 double* xi = tb->t ;
1818 double* xci = xi ; // Pointeurs
1819 double* xco = xo ; // courants
1820
1821 // k = 0
1822 int m = 0 ; // Parite de l'harmonique en phi
1823 // Dernier j: j = nt-1
1824 // Positionnement
1825 xci += nr * (nt-1) ;
1826 xco += nr * (nt-1) ;
1827
1828 // Initialisation de som
1829 for (int i=0 ; i<nr ; i++) {
1830 som[i] = 0.5*xci[i] ;
1831 xco[i] = 0. ;
1832 }
1833
1834 // j suivants
1835 for (int j=nt-2 ; j > 0 ; j--) {
1836 // Positionnement
1837 xci -= 2*nr ;
1838 xco -= nr ;
1839 // Calcul
1840 for (int i=0 ; i<nr ; i++ ) {
1841 som[i] -= 0.5*xci[i] ;
1842 xco[i] = som[i] ;
1843 } // Fin de la boucle sur r
1844 xci += nr ;
1845 for (int i=0 ; i<nr ; i++ ) {
1846 som[i] = 0.5*xci[i] ;
1847 }
1848 } // Fin de la boucle sur theta
1849 // j = 0
1850 xci -= nr ;
1851 xco -= nr ;
1852 for (int i=0; i<nr; i++) {
1853 xco[i] = som[i] ;
1854 }
1855 // Positionnement phi suivant
1856 xci += nr*nt ;
1857 xco += nr*nt ;
1858
1859 // k=1
1860 xci += nr*nt ;
1861 xco += nr*nt ;
1862
1863 for (int k=2 ; k<np+1 ; k++) {
1864 m = (k/2) % 2 ; // Parite de l'harmonique en phi
1865
1866 switch(m) {
1867 case 1: // m impair: cos(impair)
1868 // Dernier j: j = nt-1
1869 // Positionnement
1870 xci += nr * (nt-1) ;
1871 xco += nr * (nt-1) ;
1872
1873 // Initialisation de som
1874 for (int i=0 ; i<nr ; i++) {
1875 som[i] = -0.5*xci[i] ;
1876 xco[i] = 0.0;
1877 }
1878
1879 // j suivants
1880 for (int j=nt-2 ; j > 0 ; j--) {
1881 // Positionnement
1882 xci -= 2*nr ;
1883 xco -= nr ;
1884 // Calcul
1885 for (int i=0 ; i<nr ; i++ ) {
1886 som[i] += 0.5*xci[i] ;
1887 xco[i] = som[i] ;
1888 } // Fin de la boucle sur r
1889 xci +=nr ;
1890 for (int i=0 ; i<nr ; i++ ) {
1891 som[i] = -0.5*xci[i] ;
1892 }
1893 } // Fin de la boucle sur theta
1894 xci -= nr ;
1895 for (int i=0; i<nr; i++) {
1896 xco[i] += 0.5*xci[i] ;
1897 }
1898 xco -= nr ;
1899 for (int i=0; i<nr; i++) {
1900 xco[i] = 0. ;
1901 }
1902
1903 // Positionnement phi suivant
1904 xci += nr*nt ;
1905 xco += nr*nt ;
1906 break ;
1907
1908 case 0: // m pair: sin(pair)
1909 // Dernier j: j = nt-1
1910 // Positionnement
1911 xci += nr * (nt-1) ;
1912 xco += nr * (nt-1) ;
1913
1914 // Initialisation de som
1915 for (int i=0 ; i<nr ; i++) {
1916 som[i] = 0.5*xci[i] ;
1917 xco[i] = 0. ;
1918 }
1919
1920 // j suivants
1921 for (int j=nt-2 ; j > 0 ; j--) {
1922 // Positionnement
1923 xci -= 2*nr ;
1924 xco -= nr ;
1925 // Calcul
1926 for (int i=0 ; i<nr ; i++ ) {
1927 som[i] -= 0.5*xci[i] ;
1928 xco[i] = som[i] ;
1929 } // Fin de la boucle sur r
1930 xci += nr ;
1931 for (int i=0 ; i<nr ; i++ ) {
1932 som[i] = 0.5*xci[i] ;
1933 }
1934 } // Fin de la boucle sur theta
1935 // j = 0
1936 xci -= nr ;
1937 xco -= nr ;
1938 for (int i=0; i<nr; i++) {
1939 xco[i] = som[i] ;
1940 }
1941 // Positionnement phi suivant
1942 xci += nr*nt ;
1943 xco += nr*nt ;
1944 break ;
1945 } // Fin du switch
1946 } // Fin de la boucle en phi
1947
1948 // On remet les choses la ou il faut
1949 delete [] tb->t ;
1950 tb->t = xo ;
1951
1952 //Menage
1953 delete [] som ;
1954
1955 // base de developpement
1956 int base_r = b & MSQ_R ;
1957 int base_p = b & MSQ_P ;
1958 switch(base_r){
1959 case(R_CHEBPI_P):
1960 b = R_CHEBPI_I | base_p | T_COSSIN_C ;
1961 break ;
1962 case(R_CHEBPI_I):
1963 b = R_CHEBPI_P | base_p | T_COSSIN_C ;
1964 break ;
1965 default:
1966 b = base_r | base_p | T_COSSIN_C ;
1967 break;
1968 }
1969}
1970}
#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