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