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