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