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 
91 namespace Lorene {
92 void _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 
103 void _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 
261 void _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 
414 void _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 
547 void _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 
678 void _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  //---------------------
799 void _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  //----------------------
920 void _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  //----------------
1087 void _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  //----------------
1248 void _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  //----------------
1410 void _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  //----------------------
1577 void _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  //----------------
1780 void _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 }
R_CHEBPI_I
#define R_CHEBPI_I
Cheb. pair-impair suivant l impair pour l=0.
Definition: type_parite.h:174
Lorene
Lorene prototypes.
Definition: app_hor.h:67
MSQ_P
#define MSQ_P
Extraction de l&#39;info sur Phi.
Definition: type_parite.h:156
T_COS
#define T_COS
dev. cos seulement
Definition: type_parite.h:196
T_COSSIN_SP
#define T_COSSIN_SP
sin pair-cos impair alternes, sin pour m=0
Definition: type_parite.h:210
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_COS_P
#define T_COS_P
dev. cos seulement, harmoniques paires
Definition: type_parite.h:200
T_COSSIN_C
#define T_COSSIN_C
dev. cos-sin alternes, cos pour m=0
Definition: type_parite.h:192
MSQ_R
#define MSQ_R
Extraction de l&#39;info sur R.
Definition: type_parite.h:152
T_SIN_P
#define T_SIN_P
dev. sin seulement, harmoniques paires
Definition: type_parite.h:202
T_COSSIN_SI
#define T_COSSIN_SI
sin impair-cos pair alternes, sin pour m=0
Definition: type_parite.h:214
R_CHEBPI_P
#define R_CHEBPI_P
Cheb. pair-impair suivant l pair pour l=0.
Definition: type_parite.h:172
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_COSSIN_S
#define T_COSSIN_S
dev. cos-sin alternes, sin pour m=0
Definition: type_parite.h:194