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 char 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 
79 namespace Lorene {
80 
81 void _sx_1d_r_legp(int, double* , double *) ;
82 void _sx_1d_r_legi(int, double* , double *) ;
83 
84 
85  //-----------------------------------
86  // Routine pour les cas non prevus --
87  //-----------------------------------
88 
89 void _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 
100 void _sx2_identite(Tbl* , int& ) {
101  return ;
102 }
103 
104  //---------------
105  // cas R_CHEBP --
106  //--------------
107 void _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 
177 void _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 
246 void _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 
363 void _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 
482 void _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  //------------------
515 void _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 
603 void _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  //--------------
690 void _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 
748 void _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
Lorene prototypes.
Definition: app_hor.h:64