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