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