LORENE
op_sx.C
1 /*
2  * Copyright (c) 1999-2000 Jean-Alain Marck
3  * Copyright (c) 1999-2001 Eric Gourgoulhon
4  * Copyright (c) 1999-2001 Philippe Grandclement
5  *
6  * This file is part of LORENE.
7  *
8  * LORENE is free software; you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation; either version 2 of the License, or
11  * (at your option) any later version.
12  *
13  * LORENE is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with LORENE; if not, write to the Free Software
20  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21  *
22  */
23 
24 
25 
26 
27 /*
28  * Ensemble des routines de base de division par rapport a x
29  * (Utilisation interne)
30  *
31  * void _sx_XXXX(Tbl * t, int & b)
32  * t pointeur sur le Tbl d'entree-sortie
33  * b base spectrale
34  *
35  */
36 
37  /*
38  * $Id: op_sx.C,v 1.7 2024/01/26 17:44:25 g_servignat Exp $
39  * $Log: op_sx.C,v $
40  * Revision 1.7 2024/01/26 17:44:25 g_servignat
41  * Updated the Pseudopolytrope_1D class to be consistent with the paper (i.e. with a GPP in the middle)
42  *
43  * Revision 1.5 2016/12/05 16:18:08 j_novak
44  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
45  *
46  * Revision 1.4 2015/03/05 08:49:32 j_novak
47  * Implemented operators with Legendre bases.
48  *
49  * Revision 1.3 2014/10/13 08:53:26 j_novak
50  * Lorene classes and functions now belong to the namespace Lorene.
51  *
52  * Revision 1.2 2004/11/23 15:16:01 m_forot
53  *
54  * Added the bases for the cases without any equatorial symmetry
55  * (T_COSSIN_C, T_COSSIN_S, T_LEG, R_CHEBPI_P, R_CHEBPI_I).
56  *
57  * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
58  * LORENE
59  *
60  * Revision 2.5 2000/09/07 12:50:48 phil
61  * *** empty log message ***
62  *
63  * Revision 2.4 2000/02/24 16:42:59 eric
64  * Initialisation a zero du tableau xo avant le calcul.
65  *
66  * Revision 2.3 1999/11/15 16:39:17 eric
67  * Tbl::dim est desormais un Dim_tbl et non plus un Dim_tbl*.
68  *
69  * Revision 2.2 1999/10/18 13:40:11 eric
70  * Suppression de la routine _sx_r_chebu car doublon avec _sxm1_cheb.
71  *
72  * Revision 2.1 1999/09/13 11:35:42 phil
73  * gestion des cas k==1 et k==np
74  *
75  * Revision 2.0 1999/04/26 14:59:42 phil
76  * *** empty log message ***
77  *
78  *
79  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_sx.C,v 1.7 2024/01/26 17:44:25 g_servignat Exp $
80  *
81  */
82 
83  // Fichier includes
84 #include "tbl.h"
85 
86 
87  //-----------------------------------
88  // Routine pour les cas non prevus --
89  //-----------------------------------
90 
91 namespace Lorene {
92 void _sx_pas_prevu(Tbl * tb, int& base) {
93  cout << "sx pas prevu..." << endl ;
94  cout << "Tbl: " << tb << " base: " << base << endl ;
95  abort () ;
96  exit(-1) ;
97 }
98 
99  //-------------
100  // Identite ---
101  //-------------
102 
103 void _sx_identite(Tbl* , int& ) {
104  return ;
105 }
106 
107  //---------------
108  // cas R_CHEBP --
109  //--------------
110 
111 void _sx_r_chebp(Tbl* tb, int& base)
112  {
113  // Peut-etre rien a faire ?
114  if (tb->get_etat() == ETATZERO) {
115  int base_t = base & MSQ_T ;
116  int base_p = base & MSQ_P ;
117  base = base_p | base_t | R_CHEBI ;
118  return ;
119  }
120 
121  // Pour le confort
122  int nr = (tb->dim).dim[0] ; // Nombre
123  int nt = (tb->dim).dim[1] ; // de points
124  int np = (tb->dim).dim[2] ; // physiques REELS
125  np = np - 2 ; // Nombre de points physiques
126 
127  // pt. sur le tableau de double resultat
128  double* xo = new double [tb->get_taille()];
129 
130  // Initialisation a zero :
131  for (int i=0; i<tb->get_taille(); i++) {
132  xo[i] = 0 ;
133  }
134 
135  // On y va...
136  double* xi = tb->t ;
137  double* xci = xi ; // Pointeurs
138  double* xco = xo ; // courants
139 
140  int borne_phi = np + 1 ;
141  if (np == 1) {
142  borne_phi = 1 ;
143  }
144 
145  for (int k=0 ; k< borne_phi ; k++)
146  if (k==1) {
147  xci += nr*nt ;
148  xco += nr*nt ;
149  }
150  else {
151  for (int j=0 ; j<nt ; j++) {
152 
153  double som ;
154  int sgn = 1 ;
155 
156  // Regularisation (cf Jordan Nicoules thesis)
157  double val_ori = 0 ;
158  for (int i=0 ; i<nr-1 ; i++) {
159  val_ori += sgn*xci[i] ;
160  sgn = -sgn ;
161  }
162  sgn = -sgn ;
163  xci[nr-1] = sgn*val_ori ;
164  sgn = 1 ;
165 
166  xco[nr-1] = 0 ;
167  som = 2 * sgn * xci[nr-1] ;
168  xco[nr-2] = som ;
169  for (int i = nr-3 ; i >= 0 ; i-- ) {
170  sgn = - sgn ;
171  som += 2 * sgn * xci[i+1] ;
172  xco[i] = som ;
173  } // Fin de la premiere boucle sur r
174  for (int i=0 ; i<nr ; i+=2) {
175  xco[i] = -xco[i] ;
176  } // Fin de la deuxieme boucle sur r
177 
178  xci += nr ;
179  xco += nr ;
180  } // Fin de la boucle sur theta
181  } // Fin de la boucle sur phi
182 
183  // On remet les choses la ou il faut
184  delete [] tb->t ;
185  tb->t = xo ;
186 
187  // base de developpement
188  // pair -> impair
189  int base_t = base & MSQ_T ;
190  int base_p = base & MSQ_P ;
191  base = base_p | base_t | R_CHEBI ;
192 
193 }
194 
195  //----------------
196  // cas R_CHEBI ---
197  //----------------
198 
199 void _sx_r_chebi(Tbl* tb, int& base)
200 {
201 
202  // Peut-etre rien a faire ?
203  if (tb->get_etat() == ETATZERO) {
204  int base_t = base & MSQ_T ;
205  int base_p = base & MSQ_P ;
206  base = base_p | base_t | R_CHEBP ;
207  return ;
208  }
209 
210  // Pour le confort
211  int nr = (tb->dim).dim[0] ; // Nombre
212  int nt = (tb->dim).dim[1] ; // de points
213  int np = (tb->dim).dim[2] ; // physiques REELS
214  np = np - 2 ; // Nombre de points physiques
215 
216  // pt. sur le tableau de double resultat
217  double* xo = new double [tb->get_taille()];
218 
219  // Initialisation a zero :
220  for (int i=0; i<tb->get_taille(); i++) {
221  xo[i] = 0 ;
222  }
223 
224  // On y va...
225  double* xi = tb->t ;
226  double* xci = xi ; // Pointeurs
227  double* xco = xo ; // courants
228 
229  int borne_phi = np + 1 ;
230  if (np == 1) {
231  borne_phi = 1 ;
232  }
233 
234  for (int k=0 ; k< borne_phi ; k++)
235  if (k == 1) {
236  xci += nr*nt ;
237  xco += nr*nt ;
238  }
239  else {
240  for (int j=0 ; j<nt ; j++) {
241  double som ;
242  int sgn = 1 ;
243 
244  xco[nr-1] = 0 ;
245  som = 2 * sgn * xci[nr-2] ;
246  xco[nr-2] = som ;
247  for (int i = nr-3 ; i >= 0 ; i-- ) {
248  sgn = - sgn ;
249  som += 2 * sgn * xci[i] ;
250  xco[i] = som ;
251  } // Fin de la premiere boucle sur r
252  for (int i=0 ; i<nr ; i+=2) {
253  xco[i] = -xco[i] ;
254  } // Fin de la deuxieme boucle sur r
255  xco[0] *= .5 ;
256 
257  xci += nr ;
258  xco += nr ;
259  } // Fin de la boucle sur theta
260  } // Fin de la boucle sur phi
261 
262  // On remet les choses la ou il faut
263  delete [] tb->t ;
264  tb->t = xo ;
265 
266  // base de developpement
267  // impair -> pair
268  int base_t = base & MSQ_T ;
269  int base_p = base & MSQ_P ;
270  base = base_p | base_t | R_CHEBP ;
271 }
272 
273  //--------------------
274  // cas R_CHEBPIM_P --
275  //------------------
276 
277 void _sx_r_chebpim_p(Tbl* tb, int& base)
278 {
279 
280  // Peut-etre rien a faire ?
281  if (tb->get_etat() == ETATZERO) {
282  int base_t = base & MSQ_T ;
283  int base_p = base & MSQ_P ;
284  base = base_p | base_t | R_CHEBPIM_I ;
285  return ;
286  }
287 
288  // Pour le confort
289  int nr = (tb->dim).dim[0] ; // Nombre
290  int nt = (tb->dim).dim[1] ; // de points
291  int np = (tb->dim).dim[2] ; // physiques REELS
292  np = np - 2 ;
293 
294  // pt. sur le tableau de double resultat
295  double* xo = new double [tb->get_taille()];
296 
297  // Initialisation a zero :
298  for (int i=0; i<tb->get_taille(); i++) {
299  xo[i] = 0 ;
300  }
301 
302  // On y va...
303  double* xi = tb->t ;
304  double* xci ; // Pointeurs
305  double* xco ; // courants
306 
307 
308  // Partie paire
309  xci = xi ;
310  xco = xo ;
311 
312  int auxiliaire ;
313 
314  for (int k=0 ; k<np+1 ; k += 4) {
315 
316  auxiliaire = (k==np) ? 1 : 2 ; // evite de prendre le dernier coef
317  for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
318 
319 
320  // On evite le sin (0*phi)
321 
322  if ((k==0) && (kmod == 1)) {
323  xco += nr*nt ;
324  xci += nr*nt ;
325  }
326 
327  else
328  for (int j=0 ; j<nt ; j++) {
329 
330  double som ;
331  int sgn = 1 ;
332 
333  xco[nr-1] = 0 ;
334  som = 2 * sgn * xci[nr-1] ;
335  xco[nr-2] = som ;
336  for (int i = nr-3 ; i >= 0 ; i-- ) {
337  sgn = - sgn ;
338  som += 2 * sgn * xci[i+1] ;
339  xco[i] = som ;
340  } // Fin de la premiere boucle sur r
341  for (int i=0 ; i<nr ; i+=2) {
342  xco[i] = -xco[i] ;
343  } // Fin de la deuxieme boucle sur r
344 
345  xci += nr ; // au
346  xco += nr ; // suivant
347  } // Fin de la boucle sur theta
348  } // Fin de la boucle sur kmod
349  xci += 2*nr*nt ; // saute
350  xco += 2*nr*nt ; // 2 phis
351  } // Fin de la boucle sur phi
352 
353  // Partie impaire
354  xci = xi + 2*nr*nt ;
355  xco = xo + 2*nr*nt ;
356  for (int k=2 ; k<np+1 ; k += 4) {
357 
358  auxiliaire = (k==np) ? 1 : 2 ; // evite de prendre le dernier coef
359  for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
360  for (int j=0 ; j<nt ; j++) {
361 
362  double som ;
363  int sgn = 1 ;
364 
365  xco[nr-1] = 0 ;
366  som = 2 * sgn * xci[nr-2] ;
367  xco[nr-2] = som ;
368  for (int i = nr-3 ; i >= 0 ; i-- ) {
369  sgn = - sgn ;
370  som += 2 * sgn * xci[i] ;
371  xco[i] = som ;
372  } // Fin de la premiere boucle sur r
373  for (int i=0 ; i<nr ; i+=2) {
374  xco[i] = -xco[i] ;
375  } // Fin de la deuxieme boucle sur r
376  xco[0] *= .5 ;
377 
378  xci += nr ;
379  xco += nr ;
380  } // Fin de la boucle sur theta
381  } // Fin de la boucle sur kmod
382  xci += 2*nr*nt ;
383  xco += 2*nr*nt ;
384  } // Fin de la boucle sur phi
385 
386  // On remet les choses la ou il faut
387  delete [] tb->t ;
388  tb->t = xo ;
389 
390  // base de developpement
391  // (pair,impair) -> (impair,pair)
392  int base_t = base & MSQ_T ;
393  int base_p = base & MSQ_P ;
394  base = base_p | base_t | R_CHEBPIM_I ;
395 }
396 
397  //-------------------
398  // cas R_CHEBPIM_I --
399  //-------------------
400 
401 void _sx_r_chebpim_i(Tbl* tb, int& base)
402 {
403 
404  // Peut-etre rien a faire ?
405  if (tb->get_etat() == ETATZERO) {
406  int base_t = base & MSQ_T ;
407  int base_p = base & MSQ_P ;
408  base = base_p | base_t | R_CHEBPIM_P ;
409  return ;
410  }
411 
412  // Pour le confort
413  int nr = (tb->dim).dim[0] ; // Nombre
414  int nt = (tb->dim).dim[1] ; // de points
415  int np = (tb->dim).dim[2] ; // physiques REELS
416  np = np - 2 ;
417 
418  // pt. sur le tableau de double resultat
419  double* xo = new double [tb->get_taille()];
420 
421  // Initialisation a zero :
422  for (int i=0; i<tb->get_taille(); i++) {
423  xo[i] = 0 ;
424  }
425 
426  // On y va...
427  double* xi = tb->t ;
428  double* xci ; // Pointeurs
429  double* xco ; // courants
430 
431  // Partie impaire
432  xci = xi ;
433  xco = xo ;
434 
435 
436  int auxiliaire ;
437 
438  for (int k=0 ; k<np+1 ; k += 4) {
439 
440  auxiliaire = (k==np) ? 1 : 2 ; // evite de prendre le dernier coef
441  for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
442 
443 
444  // On evite le sin (0*phi)
445 
446  if ((k==0) && (kmod == 1)) {
447  xco += nr*nt ;
448  xci += nr*nt ;
449  }
450 
451  else
452  for (int j=0 ; j<nt ; j++) {
453 
454  double som ;
455  int sgn = 1 ;
456 
457  xco[nr-1] = 0 ;
458  som = 2 * sgn * xci[nr-2] ;
459  xco[nr-2] = som ;
460  for (int i = nr-3 ; i >= 0 ; i-- ) {
461  sgn = - sgn ;
462  som += 2 * sgn * xci[i] ;
463  xco[i] = som ;
464  } // Fin de la premiere boucle sur r
465  for (int i=0 ; i<nr ; i+=2) {
466  xco[i] = -xco[i] ;
467  } // Fin de la deuxieme boucle sur r
468  xco[0] *= .5 ;
469 
470  xci += nr ;
471  xco += nr ;
472  } // Fin de la boucle sur theta
473  } // Fin de la boucle sur kmod
474  xci += 2*nr*nt ;
475  xco += 2*nr*nt ;
476  } // Fin de la boucle sur phi
477 
478  // Partie paire
479  xci = xi + 2*nr*nt ;
480  xco = xo + 2*nr*nt ;
481  for (int k=2 ; k<np+1 ; k += 4) {
482 
483  auxiliaire = (k==np) ? 1 : 2 ; // evite de prendre le dernier coef
484  for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
485  for (int j=0 ; j<nt ; j++) {
486 
487  double som ;
488  int i ;
489  int sgn = 1 ;
490 
491  xco[nr-1] = 0 ;
492  som = 2 * sgn * xci[nr-1] ;
493  xco[nr-2] = som ;
494  for (i = nr-3 ; i >= 0 ; i-- ) {
495  sgn = - sgn ;
496  som += 2 * sgn * xci[i+1] ;
497  xco[i] = som ;
498  } // Fin de la premiere boucle sur r
499  for (i=0 ; i<nr ; i+=2) {
500  xco[i] = -xco[i] ;
501  } // Fin de la deuxieme boucle sur r
502 
503  xci += nr ;
504  xco += nr ;
505  } // Fin de la boucle sur theta
506  } // Fin de la boucle sur kmod
507  xci += 2*nr*nt ;
508  xco += 2*nr*nt ;
509  } // Fin de la boucle sur phi
510 
511  // On remet les choses la ou il faut
512  delete [] tb->t ;
513  tb->t = xo ;
514 
515  // base de developpement
516  // (impair,pair) -> (pair,impair)
517  int base_t = base & MSQ_T ;
518  int base_p = base & MSQ_P ;
519  base = base_p | base_t | R_CHEBPIM_P ;
520 }
521 
522  //------------------
523  // cas R_CHEBPI_P --
524  //------------------
525 
526 void _sx_r_chebpi_p(Tbl* tb, int& base)
527  {
528  // Peut-etre rien a faire ?
529  if (tb->get_etat() == ETATZERO) {
530  int base_t = base & MSQ_T ;
531  int base_p = base & MSQ_P ;
532  base = base_p | base_t | R_CHEBPI_I ;
533  return ;
534  }
535 
536  // Pour le confort
537  int nr = (tb->dim).dim[0] ; // Nombre
538  int nt = (tb->dim).dim[1] ; // de points
539  int np = (tb->dim).dim[2] ; // physiques REELS
540  np = np - 2 ; // Nombre de points physiques
541 
542  // pt. sur le tableau de double resultat
543  double* xo = new double [tb->get_taille()];
544 
545  // Initialisation a zero :
546  for (int i=0; i<tb->get_taille(); i++) {
547  xo[i] = 0 ;
548  }
549 
550  // On y va...
551  double* xi = tb->t ;
552  double* xci = xi ; // Pointeurs
553  double* xco = xo ; // courants
554 
555  int borne_phi = np + 1 ;
556  if (np == 1) {
557  borne_phi = 1 ;
558  }
559 
560  for (int k=0 ; k< borne_phi ; k++)
561  if (k==1) {
562  xci += nr*nt ;
563  xco += nr*nt ;
564  }
565  else {
566  for (int j=0 ; j<nt ; j++) {
567  int l = j%2;
568  if(l==0){
569  double som ;
570  int sgn = 1 ;
571 
572  xco[nr-1] = 0 ;
573  som = 2 * sgn * xci[nr-1] ;
574  xco[nr-2] = som ;
575  for (int i = nr-3 ; i >= 0 ; i-- ) {
576  sgn = - sgn ;
577  som += 2 * sgn * xci[i+1] ;
578  xco[i] = som ;
579  } // Fin de la premiere boucle sur r
580  for (int i=0 ; i<nr ; i+=2) {
581  xco[i] = -xco[i] ;
582  } // Fin de la deuxieme boucle sur r
583  } else {
584  double som ;
585  int sgn = 1 ;
586 
587  xco[nr-1] = 0 ;
588  som = 2 * sgn * xci[nr-2] ;
589  xco[nr-2] = som ;
590  for (int i = nr-3 ; i >= 0 ; i-- ) {
591  sgn = - sgn ;
592  som += 2 * sgn * xci[i] ;
593  xco[i] = som ;
594  } // Fin de la premiere boucle sur r
595  for (int i=0 ; i<nr ; i+=2) {
596  xco[i] = -xco[i] ;
597  } // Fin de la deuxieme boucle sur r
598  xco[0] *= .5 ;
599  }
600  xci += nr ;
601  xco += nr ;
602  } // Fin de la boucle sur theta
603  } // Fin de la boucle sur phi
604 
605  // On remet les choses la ou il faut
606  delete [] tb->t ;
607  tb->t = xo ;
608 
609  // base de developpement
610  // pair -> impair
611  int base_t = base & MSQ_T ;
612  int base_p = base & MSQ_P ;
613  base = base_p | base_t | R_CHEBPI_I ;
614 
615 }
616 
617  //-------------------
618  // cas R_CHEBPI_I ---
619  //-------------------
620 
621 void _sx_r_chebpi_i(Tbl* tb, int& base)
622 {
623 
624  // Peut-etre rien a faire ?
625  if (tb->get_etat() == ETATZERO) {
626  int base_t = base & MSQ_T ;
627  int base_p = base & MSQ_P ;
628  base = base_p | base_t | R_CHEBPI_P ;
629  return ;
630  }
631 
632  // Pour le confort
633  int nr = (tb->dim).dim[0] ; // Nombre
634  int nt = (tb->dim).dim[1] ; // de points
635  int np = (tb->dim).dim[2] ; // physiques REELS
636  np = np - 2 ; // Nombre de points physiques
637 
638  // pt. sur le tableau de double resultat
639  double* xo = new double [tb->get_taille()];
640 
641  // Initialisation a zero :
642  for (int i=0; i<tb->get_taille(); i++) {
643  xo[i] = 0 ;
644  }
645 
646  // On y va...
647  double* xi = tb->t ;
648  double* xci = xi ; // Pointeurs
649  double* xco = xo ; // courants
650 
651  int borne_phi = np + 1 ;
652  if (np == 1) {
653  borne_phi = 1 ;
654  }
655 
656  for (int k=0 ; k< borne_phi ; k++)
657  if (k == 1) {
658  xci += nr*nt ;
659  xco += nr*nt ;
660  }
661  else {
662  for (int j=0 ; j<nt ; j++) {
663  int l = j%2;
664  if(l==1){
665  double som ;
666  int sgn = 1 ;
667 
668  xco[nr-1] = 0 ;
669  som = 2 * sgn * xci[nr-1] ;
670  xco[nr-2] = som ;
671  for (int i = nr-3 ; i >= 0 ; i-- ) {
672  sgn = - sgn ;
673  som += 2 * sgn * xci[i+1] ;
674  xco[i] = som ;
675  } // Fin de la premiere boucle sur r
676  for (int i=0 ; i<nr ; i+=2) {
677  xco[i] = -xco[i] ;
678  } // Fin de la deuxieme boucle sur r
679  } else {
680  double som ;
681  int sgn = 1 ;
682 
683  xco[nr-1] = 0 ;
684  som = 2 * sgn * xci[nr-2] ;
685  xco[nr-2] = som ;
686  for (int i = nr-3 ; i >= 0 ; i-- ) {
687  sgn = - sgn ;
688  som += 2 * sgn * xci[i] ;
689  xco[i] = som ;
690  } // Fin de la premiere boucle sur r
691  for (int i=0 ; i<nr ; i+=2) {
692  xco[i] = -xco[i] ;
693  } // Fin de la deuxieme boucle sur r
694  xco[0] *= .5 ;
695  }
696  xci += nr ;
697  xco += nr ;
698  } // Fin de la boucle sur theta
699  } // Fin de la boucle sur phi
700 
701  // On remet les choses la ou il faut
702  delete [] tb->t ;
703  tb->t = xo ;
704 
705  // base de developpement
706  // impair -> pair
707  int base_t = base & MSQ_T ;
708  int base_p = base & MSQ_P ;
709  base = base_p | base_t | R_CHEBPI_P ;
710 }
711 
712  //--------------
713  // cas R_LEGP --
714  //--------------
715 
716 void _sx_r_legp(Tbl* tb, int& base)
717 
718  {
719  // Peut-etre rien a faire ?
720  if (tb->get_etat() == ETATZERO) {
721  int base_t = base & MSQ_T ;
722  int base_p = base & MSQ_P ;
723  base = base_p | base_t | R_LEGI ;
724  return ;
725  }
726 
727  // Pour le confort
728  int nr = (tb->dim).dim[0] ; // Nombre
729  int nt = (tb->dim).dim[1] ; // de points
730  int np = (tb->dim).dim[2] ; // physiques REELS
731  np = np - 2 ; // Nombre de points physiques
732 
733  // pt. sur le tableau de double resultat
734  double* xo = new double [tb->get_taille()];
735 
736  // Initialisation a zero :
737  for (int i=0; i<tb->get_taille(); i++) {
738  xo[i] = 0 ;
739  }
740 
741  // On y va...
742  double* xi = tb->t ;
743  double* xci = xi ; // Pointeurs
744  double* xco = xo ; // courants
745 
746  int borne_phi = np + 1 ;
747  if (np == 1) {
748  borne_phi = 1 ;
749  }
750 
751  for (int k=0 ; k< borne_phi ; k++)
752  if (k==1) {
753  xci += nr*nt ;
754  xco += nr*nt ;
755  }
756  else {
757  for (int j=0 ; j<nt ; j++) {
758 
759  double som = 0 ;
760 
761  xco[nr-1] = 0 ;
762  for (int i=nr - 2; i>=0; i--) {
763  som += xci[i+1] ;
764  xco[i] = double(4*i+3)/double(2*i+2)*som ;
765  som *= -double(2*i+1)/double(2*i+2) ;
766  } //Fin de la boucle sur r
767 
768  xci += nr ;
769  xco += nr ;
770  } // Fin de la boucle sur theta
771  } // Fin de la boucle sur phi
772 
773  // On remet les choses la ou il faut
774  delete [] tb->t ;
775  tb->t = xo ;
776 
777  // base de developpement
778  // pair -> impair
779  int base_t = base & MSQ_T ;
780  int base_p = base & MSQ_P ;
781  base = base_p | base_t | R_LEGI ;
782 
783 }
784 
785  //---------------
786  // cas R_LEGI ---
787  //---------------
788 
789 void _sx_r_legi(Tbl* tb, int& base)
790 {
791 
792  // Peut-etre rien a faire ?
793  if (tb->get_etat() == ETATZERO) {
794  int base_t = base & MSQ_T ;
795  int base_p = base & MSQ_P ;
796  base = base_p | base_t | R_LEGP ;
797  return ;
798  }
799 
800  // Pour le confort
801  int nr = (tb->dim).dim[0] ; // Nombre
802  int nt = (tb->dim).dim[1] ; // de points
803  int np = (tb->dim).dim[2] ; // physiques REELS
804  np = np - 2 ; // Nombre de points physiques
805 
806  // pt. sur le tableau de double resultat
807  double* xo = new double [tb->get_taille()];
808 
809  // Initialisation a zero :
810  for (int i=0; i<tb->get_taille(); i++) {
811  xo[i] = 0 ;
812  }
813 
814  // On y va...
815  double* xi = tb->t ;
816  double* xci = xi ; // Pointeurs
817  double* xco = xo ; // courants
818 
819  int borne_phi = np + 1 ;
820  if (np == 1) {
821  borne_phi = 1 ;
822  }
823 
824  for (int k=0 ; k< borne_phi ; k++)
825  if (k == 1) {
826  xci += nr*nt ;
827  xco += nr*nt ;
828  }
829  else {
830  for (int j=0 ; j<nt ; j++) {
831  double som = 0 ;
832 
833  xco[nr-1] = 0 ;
834  for (int i = nr-2 ; i >= 0 ; i-- ) {
835  som += xci[i] ;
836  xco[i] = double(4*i+1)/double(2*i+1)*som ;
837  som *= -double(2*i)/double(2*i+1) ;
838  } // Fin de la boucle sur r
839 
840  xci += nr ;
841  xco += nr ;
842  } // Fin de la boucle sur theta
843  } // Fin de la boucle sur phi
844 
845  // On remet les choses la ou il faut
846  delete [] tb->t ;
847  tb->t = xo ;
848 
849  // base de developpement
850  // impair -> pair
851  int base_t = base & MSQ_T ;
852  int base_p = base & MSQ_P ;
853  base = base_p | base_t | R_LEGP ;
854 }
855 
856 
857 }
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
R_LEGP
#define R_LEGP
base de Legendre paire (rare) seulement
Definition: type_parite.h:184
R_LEGI
#define R_LEGI
base de Legendre impaire (rare) seulement
Definition: type_parite.h:186
R_CHEBI
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
MSQ_T
#define MSQ_T
Extraction de l&#39;info sur Theta.
Definition: type_parite.h:154
R_CHEBPIM_I
#define R_CHEBPIM_I
Cheb. pair-impair suivant m, impair pour m=0.
Definition: type_parite.h:178
R_CHEBPIM_P
#define R_CHEBPIM_P
Cheb. pair-impair suivant m, pair pour m=0.
Definition: type_parite.h:176
R_CHEBPI_P
#define R_CHEBPI_P
Cheb. pair-impair suivant l pair pour l=0.
Definition: type_parite.h:172