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