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.6 2023/05/22 15:44:21 g_servignat Exp $
39  * $Log: op_sx.C,v $
40  * Revision 1.6 2023/05/22 15:44:21 g_servignat
41  * Added regularisation in division by \xi in a nucleus according to J. Nicoules' thesis
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.6 2023/05/22 15:44:21 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 ; ++i) {
159  val_ori += sgn*xci[i] ;
160  sgn = -sgn ;
161  }
162  xci[nr-1] += sgn*val_ori ;
163  sgn = 1 ;
164 
165  xco[nr-1] = 0 ;
166  som = 2 * sgn * xci[nr-1] ;
167  xco[nr-2] = som ;
168  for (int i = nr-3 ; i >= 0 ; i-- ) {
169  sgn = - sgn ;
170  som += 2 * sgn * xci[i+1] ;
171  xco[i] = som ;
172  } // Fin de la premiere boucle sur r
173  for (int i=0 ; i<nr ; i+=2) {
174  xco[i] = -xco[i] ;
175  } // Fin de la deuxieme boucle sur r
176 
177  xci += nr ;
178  xco += nr ;
179  } // Fin de la boucle sur theta
180  } // Fin de la boucle sur phi
181 
182  // On remet les choses la ou il faut
183  delete [] tb->t ;
184  tb->t = xo ;
185 
186  // base de developpement
187  // pair -> impair
188  int base_t = base & MSQ_T ;
189  int base_p = base & MSQ_P ;
190  base = base_p | base_t | R_CHEBI ;
191 
192 }
193 
194  //----------------
195  // cas R_CHEBI ---
196  //----------------
197 
198 void _sx_r_chebi(Tbl* tb, int& base)
199 {
200 
201  // Peut-etre rien a faire ?
202  if (tb->get_etat() == ETATZERO) {
203  int base_t = base & MSQ_T ;
204  int base_p = base & MSQ_P ;
205  base = base_p | base_t | R_CHEBP ;
206  return ;
207  }
208 
209  // Pour le confort
210  int nr = (tb->dim).dim[0] ; // Nombre
211  int nt = (tb->dim).dim[1] ; // de points
212  int np = (tb->dim).dim[2] ; // physiques REELS
213  np = np - 2 ; // Nombre de points physiques
214 
215  // pt. sur le tableau de double resultat
216  double* xo = new double [tb->get_taille()];
217 
218  // Initialisation a zero :
219  for (int i=0; i<tb->get_taille(); i++) {
220  xo[i] = 0 ;
221  }
222 
223  // On y va...
224  double* xi = tb->t ;
225  double* xci = xi ; // Pointeurs
226  double* xco = xo ; // courants
227 
228  int borne_phi = np + 1 ;
229  if (np == 1) {
230  borne_phi = 1 ;
231  }
232 
233  for (int k=0 ; k< borne_phi ; k++)
234  if (k == 1) {
235  xci += nr*nt ;
236  xco += nr*nt ;
237  }
238  else {
239  for (int j=0 ; j<nt ; j++) {
240  double som ;
241  int sgn = 1 ;
242 
243  xco[nr-1] = 0 ;
244  som = 2 * sgn * xci[nr-2] ;
245  xco[nr-2] = som ;
246  for (int i = nr-3 ; i >= 0 ; i-- ) {
247  sgn = - sgn ;
248  som += 2 * sgn * xci[i] ;
249  xco[i] = som ;
250  } // Fin de la premiere boucle sur r
251  for (int i=0 ; i<nr ; i+=2) {
252  xco[i] = -xco[i] ;
253  } // Fin de la deuxieme boucle sur r
254  xco[0] *= .5 ;
255 
256  xci += nr ;
257  xco += nr ;
258  } // Fin de la boucle sur theta
259  } // Fin de la boucle sur phi
260 
261  // On remet les choses la ou il faut
262  delete [] tb->t ;
263  tb->t = xo ;
264 
265  // base de developpement
266  // impair -> pair
267  int base_t = base & MSQ_T ;
268  int base_p = base & MSQ_P ;
269  base = base_p | base_t | R_CHEBP ;
270 }
271 
272  //--------------------
273  // cas R_CHEBPIM_P --
274  //------------------
275 
276 void _sx_r_chebpim_p(Tbl* tb, int& base)
277 {
278 
279  // Peut-etre rien a faire ?
280  if (tb->get_etat() == ETATZERO) {
281  int base_t = base & MSQ_T ;
282  int base_p = base & MSQ_P ;
283  base = base_p | base_t | R_CHEBPIM_I ;
284  return ;
285  }
286 
287  // Pour le confort
288  int nr = (tb->dim).dim[0] ; // Nombre
289  int nt = (tb->dim).dim[1] ; // de points
290  int np = (tb->dim).dim[2] ; // physiques REELS
291  np = np - 2 ;
292 
293  // pt. sur le tableau de double resultat
294  double* xo = new double [tb->get_taille()];
295 
296  // Initialisation a zero :
297  for (int i=0; i<tb->get_taille(); i++) {
298  xo[i] = 0 ;
299  }
300 
301  // On y va...
302  double* xi = tb->t ;
303  double* xci ; // Pointeurs
304  double* xco ; // courants
305 
306 
307  // Partie paire
308  xci = xi ;
309  xco = xo ;
310 
311  int auxiliaire ;
312 
313  for (int k=0 ; k<np+1 ; k += 4) {
314 
315  auxiliaire = (k==np) ? 1 : 2 ; // evite de prendre le dernier coef
316  for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
317 
318 
319  // On evite le sin (0*phi)
320 
321  if ((k==0) && (kmod == 1)) {
322  xco += nr*nt ;
323  xci += nr*nt ;
324  }
325 
326  else
327  for (int j=0 ; j<nt ; j++) {
328 
329  double som ;
330  int sgn = 1 ;
331 
332  xco[nr-1] = 0 ;
333  som = 2 * sgn * xci[nr-1] ;
334  xco[nr-2] = som ;
335  for (int i = nr-3 ; i >= 0 ; i-- ) {
336  sgn = - sgn ;
337  som += 2 * sgn * xci[i+1] ;
338  xco[i] = som ;
339  } // Fin de la premiere boucle sur r
340  for (int i=0 ; i<nr ; i+=2) {
341  xco[i] = -xco[i] ;
342  } // Fin de la deuxieme boucle sur r
343 
344  xci += nr ; // au
345  xco += nr ; // suivant
346  } // Fin de la boucle sur theta
347  } // Fin de la boucle sur kmod
348  xci += 2*nr*nt ; // saute
349  xco += 2*nr*nt ; // 2 phis
350  } // Fin de la boucle sur phi
351 
352  // Partie impaire
353  xci = xi + 2*nr*nt ;
354  xco = xo + 2*nr*nt ;
355  for (int k=2 ; k<np+1 ; k += 4) {
356 
357  auxiliaire = (k==np) ? 1 : 2 ; // evite de prendre le dernier coef
358  for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
359  for (int j=0 ; j<nt ; j++) {
360 
361  double som ;
362  int sgn = 1 ;
363 
364  xco[nr-1] = 0 ;
365  som = 2 * sgn * xci[nr-2] ;
366  xco[nr-2] = som ;
367  for (int i = nr-3 ; i >= 0 ; i-- ) {
368  sgn = - sgn ;
369  som += 2 * sgn * xci[i] ;
370  xco[i] = som ;
371  } // Fin de la premiere boucle sur r
372  for (int i=0 ; i<nr ; i+=2) {
373  xco[i] = -xco[i] ;
374  } // Fin de la deuxieme boucle sur r
375  xco[0] *= .5 ;
376 
377  xci += nr ;
378  xco += nr ;
379  } // Fin de la boucle sur theta
380  } // Fin de la boucle sur kmod
381  xci += 2*nr*nt ;
382  xco += 2*nr*nt ;
383  } // Fin de la boucle sur phi
384 
385  // On remet les choses la ou il faut
386  delete [] tb->t ;
387  tb->t = xo ;
388 
389  // base de developpement
390  // (pair,impair) -> (impair,pair)
391  int base_t = base & MSQ_T ;
392  int base_p = base & MSQ_P ;
393  base = base_p | base_t | R_CHEBPIM_I ;
394 }
395 
396  //-------------------
397  // cas R_CHEBPIM_I --
398  //-------------------
399 
400 void _sx_r_chebpim_i(Tbl* tb, int& base)
401 {
402 
403  // Peut-etre rien a faire ?
404  if (tb->get_etat() == ETATZERO) {
405  int base_t = base & MSQ_T ;
406  int base_p = base & MSQ_P ;
407  base = base_p | base_t | R_CHEBPIM_P ;
408  return ;
409  }
410 
411  // Pour le confort
412  int nr = (tb->dim).dim[0] ; // Nombre
413  int nt = (tb->dim).dim[1] ; // de points
414  int np = (tb->dim).dim[2] ; // physiques REELS
415  np = np - 2 ;
416 
417  // pt. sur le tableau de double resultat
418  double* xo = new double [tb->get_taille()];
419 
420  // Initialisation a zero :
421  for (int i=0; i<tb->get_taille(); i++) {
422  xo[i] = 0 ;
423  }
424 
425  // On y va...
426  double* xi = tb->t ;
427  double* xci ; // Pointeurs
428  double* xco ; // courants
429 
430  // Partie impaire
431  xci = xi ;
432  xco = xo ;
433 
434 
435  int auxiliaire ;
436 
437  for (int k=0 ; k<np+1 ; k += 4) {
438 
439  auxiliaire = (k==np) ? 1 : 2 ; // evite de prendre le dernier coef
440  for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
441 
442 
443  // On evite le sin (0*phi)
444 
445  if ((k==0) && (kmod == 1)) {
446  xco += nr*nt ;
447  xci += nr*nt ;
448  }
449 
450  else
451  for (int j=0 ; j<nt ; j++) {
452 
453  double som ;
454  int sgn = 1 ;
455 
456  xco[nr-1] = 0 ;
457  som = 2 * sgn * xci[nr-2] ;
458  xco[nr-2] = som ;
459  for (int i = nr-3 ; i >= 0 ; i-- ) {
460  sgn = - sgn ;
461  som += 2 * sgn * xci[i] ;
462  xco[i] = som ;
463  } // Fin de la premiere boucle sur r
464  for (int i=0 ; i<nr ; i+=2) {
465  xco[i] = -xco[i] ;
466  } // Fin de la deuxieme boucle sur r
467  xco[0] *= .5 ;
468 
469  xci += nr ;
470  xco += nr ;
471  } // Fin de la boucle sur theta
472  } // Fin de la boucle sur kmod
473  xci += 2*nr*nt ;
474  xco += 2*nr*nt ;
475  } // Fin de la boucle sur phi
476 
477  // Partie paire
478  xci = xi + 2*nr*nt ;
479  xco = xo + 2*nr*nt ;
480  for (int k=2 ; k<np+1 ; k += 4) {
481 
482  auxiliaire = (k==np) ? 1 : 2 ; // evite de prendre le dernier coef
483  for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
484  for (int j=0 ; j<nt ; j++) {
485 
486  double som ;
487  int i ;
488  int sgn = 1 ;
489 
490  xco[nr-1] = 0 ;
491  som = 2 * sgn * xci[nr-1] ;
492  xco[nr-2] = som ;
493  for (i = nr-3 ; i >= 0 ; i-- ) {
494  sgn = - sgn ;
495  som += 2 * sgn * xci[i+1] ;
496  xco[i] = som ;
497  } // Fin de la premiere boucle sur r
498  for (i=0 ; i<nr ; i+=2) {
499  xco[i] = -xco[i] ;
500  } // Fin de la deuxieme boucle sur r
501 
502  xci += nr ;
503  xco += nr ;
504  } // Fin de la boucle sur theta
505  } // Fin de la boucle sur kmod
506  xci += 2*nr*nt ;
507  xco += 2*nr*nt ;
508  } // Fin de la boucle sur phi
509 
510  // On remet les choses la ou il faut
511  delete [] tb->t ;
512  tb->t = xo ;
513 
514  // base de developpement
515  // (impair,pair) -> (pair,impair)
516  int base_t = base & MSQ_T ;
517  int base_p = base & MSQ_P ;
518  base = base_p | base_t | R_CHEBPIM_P ;
519 }
520 
521  //------------------
522  // cas R_CHEBPI_P --
523  //------------------
524 
525 void _sx_r_chebpi_p(Tbl* tb, int& base)
526  {
527  // Peut-etre rien a faire ?
528  if (tb->get_etat() == ETATZERO) {
529  int base_t = base & MSQ_T ;
530  int base_p = base & MSQ_P ;
531  base = base_p | base_t | R_CHEBPI_I ;
532  return ;
533  }
534 
535  // Pour le confort
536  int nr = (tb->dim).dim[0] ; // Nombre
537  int nt = (tb->dim).dim[1] ; // de points
538  int np = (tb->dim).dim[2] ; // physiques REELS
539  np = np - 2 ; // Nombre de points physiques
540 
541  // pt. sur le tableau de double resultat
542  double* xo = new double [tb->get_taille()];
543 
544  // Initialisation a zero :
545  for (int i=0; i<tb->get_taille(); i++) {
546  xo[i] = 0 ;
547  }
548 
549  // On y va...
550  double* xi = tb->t ;
551  double* xci = xi ; // Pointeurs
552  double* xco = xo ; // courants
553 
554  int borne_phi = np + 1 ;
555  if (np == 1) {
556  borne_phi = 1 ;
557  }
558 
559  for (int k=0 ; k< borne_phi ; k++)
560  if (k==1) {
561  xci += nr*nt ;
562  xco += nr*nt ;
563  }
564  else {
565  for (int j=0 ; j<nt ; j++) {
566  int l = j%2;
567  if(l==0){
568  double som ;
569  int sgn = 1 ;
570 
571  xco[nr-1] = 0 ;
572  som = 2 * sgn * xci[nr-1] ;
573  xco[nr-2] = som ;
574  for (int i = nr-3 ; i >= 0 ; i-- ) {
575  sgn = - sgn ;
576  som += 2 * sgn * xci[i+1] ;
577  xco[i] = som ;
578  } // Fin de la premiere boucle sur r
579  for (int i=0 ; i<nr ; i+=2) {
580  xco[i] = -xco[i] ;
581  } // Fin de la deuxieme boucle sur r
582  } else {
583  double som ;
584  int sgn = 1 ;
585 
586  xco[nr-1] = 0 ;
587  som = 2 * sgn * xci[nr-2] ;
588  xco[nr-2] = som ;
589  for (int i = nr-3 ; i >= 0 ; i-- ) {
590  sgn = - sgn ;
591  som += 2 * sgn * xci[i] ;
592  xco[i] = som ;
593  } // Fin de la premiere boucle sur r
594  for (int i=0 ; i<nr ; i+=2) {
595  xco[i] = -xco[i] ;
596  } // Fin de la deuxieme boucle sur r
597  xco[0] *= .5 ;
598  }
599  xci += nr ;
600  xco += nr ;
601  } // Fin de la boucle sur theta
602  } // Fin de la boucle sur phi
603 
604  // On remet les choses la ou il faut
605  delete [] tb->t ;
606  tb->t = xo ;
607 
608  // base de developpement
609  // pair -> impair
610  int base_t = base & MSQ_T ;
611  int base_p = base & MSQ_P ;
612  base = base_p | base_t | R_CHEBPI_I ;
613 
614 }
615 
616  //-------------------
617  // cas R_CHEBPI_I ---
618  //-------------------
619 
620 void _sx_r_chebpi_i(Tbl* tb, int& base)
621 {
622 
623  // Peut-etre rien a faire ?
624  if (tb->get_etat() == ETATZERO) {
625  int base_t = base & MSQ_T ;
626  int base_p = base & MSQ_P ;
627  base = base_p | base_t | R_CHEBPI_P ;
628  return ;
629  }
630 
631  // Pour le confort
632  int nr = (tb->dim).dim[0] ; // Nombre
633  int nt = (tb->dim).dim[1] ; // de points
634  int np = (tb->dim).dim[2] ; // physiques REELS
635  np = np - 2 ; // Nombre de points physiques
636 
637  // pt. sur le tableau de double resultat
638  double* xo = new double [tb->get_taille()];
639 
640  // Initialisation a zero :
641  for (int i=0; i<tb->get_taille(); i++) {
642  xo[i] = 0 ;
643  }
644 
645  // On y va...
646  double* xi = tb->t ;
647  double* xci = xi ; // Pointeurs
648  double* xco = xo ; // courants
649 
650  int borne_phi = np + 1 ;
651  if (np == 1) {
652  borne_phi = 1 ;
653  }
654 
655  for (int k=0 ; k< borne_phi ; k++)
656  if (k == 1) {
657  xci += nr*nt ;
658  xco += nr*nt ;
659  }
660  else {
661  for (int j=0 ; j<nt ; j++) {
662  int l = j%2;
663  if(l==1){
664  double som ;
665  int sgn = 1 ;
666 
667  xco[nr-1] = 0 ;
668  som = 2 * sgn * xci[nr-1] ;
669  xco[nr-2] = som ;
670  for (int i = nr-3 ; i >= 0 ; i-- ) {
671  sgn = - sgn ;
672  som += 2 * sgn * xci[i+1] ;
673  xco[i] = som ;
674  } // Fin de la premiere boucle sur r
675  for (int i=0 ; i<nr ; i+=2) {
676  xco[i] = -xco[i] ;
677  } // Fin de la deuxieme boucle sur r
678  } else {
679  double som ;
680  int sgn = 1 ;
681 
682  xco[nr-1] = 0 ;
683  som = 2 * sgn * xci[nr-2] ;
684  xco[nr-2] = som ;
685  for (int i = nr-3 ; i >= 0 ; i-- ) {
686  sgn = - sgn ;
687  som += 2 * sgn * xci[i] ;
688  xco[i] = som ;
689  } // Fin de la premiere boucle sur r
690  for (int i=0 ; i<nr ; i+=2) {
691  xco[i] = -xco[i] ;
692  } // Fin de la deuxieme boucle sur r
693  xco[0] *= .5 ;
694  }
695  xci += nr ;
696  xco += nr ;
697  } // Fin de la boucle sur theta
698  } // Fin de la boucle sur phi
699 
700  // On remet les choses la ou il faut
701  delete [] tb->t ;
702  tb->t = xo ;
703 
704  // base de developpement
705  // impair -> pair
706  int base_t = base & MSQ_T ;
707  int base_p = base & MSQ_P ;
708  base = base_p | base_t | R_CHEBPI_P ;
709 }
710 
711  //--------------
712  // cas R_LEGP --
713  //--------------
714 
715 void _sx_r_legp(Tbl* tb, int& base)
716 
717  {
718  // Peut-etre rien a faire ?
719  if (tb->get_etat() == ETATZERO) {
720  int base_t = base & MSQ_T ;
721  int base_p = base & MSQ_P ;
722  base = base_p | base_t | R_LEGI ;
723  return ;
724  }
725 
726  // Pour le confort
727  int nr = (tb->dim).dim[0] ; // Nombre
728  int nt = (tb->dim).dim[1] ; // de points
729  int np = (tb->dim).dim[2] ; // physiques REELS
730  np = np - 2 ; // Nombre de points physiques
731 
732  // pt. sur le tableau de double resultat
733  double* xo = new double [tb->get_taille()];
734 
735  // Initialisation a zero :
736  for (int i=0; i<tb->get_taille(); i++) {
737  xo[i] = 0 ;
738  }
739 
740  // On y va...
741  double* xi = tb->t ;
742  double* xci = xi ; // Pointeurs
743  double* xco = xo ; // courants
744 
745  int borne_phi = np + 1 ;
746  if (np == 1) {
747  borne_phi = 1 ;
748  }
749 
750  for (int k=0 ; k< borne_phi ; k++)
751  if (k==1) {
752  xci += nr*nt ;
753  xco += nr*nt ;
754  }
755  else {
756  for (int j=0 ; j<nt ; j++) {
757 
758  double som = 0 ;
759 
760  xco[nr-1] = 0 ;
761  for (int i=nr - 2; i>=0; i--) {
762  som += xci[i+1] ;
763  xco[i] = double(4*i+3)/double(2*i+2)*som ;
764  som *= -double(2*i+1)/double(2*i+2) ;
765  } //Fin de la boucle sur r
766 
767  xci += nr ;
768  xco += nr ;
769  } // Fin de la boucle sur theta
770  } // Fin de la boucle sur phi
771 
772  // On remet les choses la ou il faut
773  delete [] tb->t ;
774  tb->t = xo ;
775 
776  // base de developpement
777  // pair -> impair
778  int base_t = base & MSQ_T ;
779  int base_p = base & MSQ_P ;
780  base = base_p | base_t | R_LEGI ;
781 
782 }
783 
784  //---------------
785  // cas R_LEGI ---
786  //---------------
787 
788 void _sx_r_legi(Tbl* tb, int& base)
789 {
790 
791  // Peut-etre rien a faire ?
792  if (tb->get_etat() == ETATZERO) {
793  int base_t = base & MSQ_T ;
794  int base_p = base & MSQ_P ;
795  base = base_p | base_t | R_LEGP ;
796  return ;
797  }
798 
799  // Pour le confort
800  int nr = (tb->dim).dim[0] ; // Nombre
801  int nt = (tb->dim).dim[1] ; // de points
802  int np = (tb->dim).dim[2] ; // physiques REELS
803  np = np - 2 ; // Nombre de points physiques
804 
805  // pt. sur le tableau de double resultat
806  double* xo = new double [tb->get_taille()];
807 
808  // Initialisation a zero :
809  for (int i=0; i<tb->get_taille(); i++) {
810  xo[i] = 0 ;
811  }
812 
813  // On y va...
814  double* xi = tb->t ;
815  double* xci = xi ; // Pointeurs
816  double* xco = xo ; // courants
817 
818  int borne_phi = np + 1 ;
819  if (np == 1) {
820  borne_phi = 1 ;
821  }
822 
823  for (int k=0 ; k< borne_phi ; k++)
824  if (k == 1) {
825  xci += nr*nt ;
826  xco += nr*nt ;
827  }
828  else {
829  for (int j=0 ; j<nt ; j++) {
830  double som = 0 ;
831 
832  xco[nr-1] = 0 ;
833  for (int i = nr-2 ; i >= 0 ; i-- ) {
834  som += xci[i] ;
835  xco[i] = double(4*i+1)/double(2*i+1)*som ;
836  som *= -double(2*i)/double(2*i+1) ;
837  } // Fin de la boucle sur r
838 
839  xci += nr ;
840  xco += nr ;
841  } // Fin de la boucle sur theta
842  } // Fin de la boucle sur phi
843 
844  // On remet les choses la ou il faut
845  delete [] tb->t ;
846  tb->t = xo ;
847 
848  // base de developpement
849  // impair -> pair
850  int base_t = base & MSQ_T ;
851  int base_p = base & MSQ_P ;
852  base = base_p | base_t | R_LEGP ;
853 }
854 
855 
856 }
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