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