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 
82 namespace Lorene {
83 
84 void _sx_1d_r_legp(int, double* , double *) ;
85 void _sx_1d_r_legi(int, double* , double *) ;
86 
87 
88  //-----------------------------------
89  // Routine pour les cas non prevus --
90  //-----------------------------------
91 
92 void _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 
103 void _sx2_identite(Tbl* , int& ) {
104  return ;
105 }
106 
107  //---------------
108  // cas R_CHEBP --
109  //--------------
110 void _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 
180 void _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 
249 void _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 
366 void _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 
485 void _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  //------------------
518 void _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 
606 void _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  //--------------
693 void _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 
751 void _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
Lorene prototypes.
Definition: app_hor.h:67