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