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