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
23char 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
80namespace Lorene {
81void _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
92void _mult_x_identite(Tbl* , int& ) {
93 return ;
94}
95
96 //---------------
97 // cas R_CHEBP --
98 //--------------
99
100void _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
169void _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
238void _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
344void _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
449void _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
525void _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
601void _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
664void _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
734void _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
#define R_LEGI
base de Legendre impaire (rare) seulement
#define R_CHEBI
base de Cheb. impaire (rare) seulement
#define R_CHEBPIM_I
Cheb. pair-impair suivant m, impair pour m=0.
#define R_CHEBPI_I
Cheb. pair-impair suivant l impair pour l=0.
#define R_CHEBPIM_P
Cheb. pair-impair suivant m, pair pour m=0.
#define MSQ_T
Extraction de l'info sur Theta.
#define R_CHEBP
base de Cheb. paire (rare) seulement
#define MSQ_P
Extraction de l'info sur Phi.
#define R_CHEBPI_P
Cheb. pair-impair suivant l pair pour l=0.
Lorene prototypes.
Definition app_hor.h:64