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