LORENE
op_dsdx.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 derivation par rapport a r
29  * (Utilisation interne)
30  *
31  * void _dsdx_XXXX(Tbl * t, int & b)
32  * t pointeur sur le Tbl d'entree-sortie
33  * b base spectrale
34  *
35  */
36 
37 /*
38  * $Id: op_dsdx.C,v 1.9 2016/12/05 16:18:07 j_novak Exp $
39  * $Log: op_dsdx.C,v $
40  * Revision 1.9 2016/12/05 16:18:07 j_novak
41  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
42  *
43  * Revision 1.8 2014/10/13 08:53:25 j_novak
44  * Lorene classes and functions now belong to the namespace Lorene.
45  *
46  * Revision 1.7 2013/06/14 15:54:06 j_novak
47  * Inclusion of Legendre bases.
48  *
49  * Revision 1.6 2007/12/21 13:59:02 j_novak
50  * Suppression of call to pow(-1, something).
51  *
52  * Revision 1.5 2007/12/20 09:11:08 jl_cornou
53  * Correction of an error in op_sxpun about Jacobi(0,2) polynomials
54  *
55  * Revision 1.4 2007/12/11 15:28:18 jl_cornou
56  * Jacobi(0,2) polynomials partially implemented
57  *
58  * Revision 1.3 2006/05/17 13:15:18 j_novak
59  * Added a test for the pure angular grid case (nr=1), in the shell.
60  *
61  * Revision 1.2 2004/11/23 15:16:01 m_forot
62  *
63  * Added the bases for the cases without any equatorial symmetry
64  * (T_COSSIN_C, T_COSSIN_S, T_LEG, R_CHEBPI_P, R_CHEBPI_I).
65  *
66  * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
67  * LORENE
68  *
69  * Revision 2.6 2000/09/07 12:49:04 phil
70  * *** empty log message ***
71  *
72  * Revision 2.5 2000/02/24 16:41:25 eric
73  * Initialisation a zero du tableau xo avant le calcul.
74  *
75  * Revision 2.4 1999/11/15 16:38:26 eric
76  * Tbl::dim est desormais un Dim_tbl et non plus un Dim_tbl*.
77  *
78  * Revision 2.3 1999/09/13 11:31:49 phil
79  * gestion des cas k==1 et k==np+1\
80  *
81  * Revision 2.2 1999/02/22 17:25:15 hyc
82  * *** empty log message ***
83  *
84  * Revision 2.1 1999/02/22 16:17:36 hyc
85  * *** empty log message ***
86  *
87  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_dsdx.C,v 1.9 2016/12/05 16:18:07 j_novak Exp $
88  *
89  */
90 
91 // Fichier includes
92 #include "tbl.h"
93 #include <cmath>
94 
95 // Routine pour les cas non prevus
96 //--------------------------------
97 namespace Lorene {
98 void _dsdx_pas_prevu(Tbl* , int & b) {
99  cout << "dsdx pas prevu..." << endl ;
100  cout << " base: " << b << endl ;
101  abort () ;
102 }
103 
104 // cas R_CHEB
105 //-----------
106 void _dsdx_r_cheb(Tbl *tb, int & )
107 {
108 
109  // Peut-etre rien a faire ?
110  if (tb->get_etat() == ETATZERO) {
111  return ;
112  }
113 
114  // Protection
115  assert(tb->get_etat() == ETATQCQ) ;
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->dim).taille] ;
125 
126  // Initialisation a zero :
127  for (int i=0; i<(tb->dim).taille; i++) {
128  xo[i] = 0 ;
129  }
130  if (nr > 2) { // If not an angular grid...
131  // On y va...
132  double* xi = tb->t ;
133  double* xci = xi ; // Pointeurs
134  double* xco = xo ; // courants
135 
136  int borne_phi = np + 1 ;
137  if (np == 1) borne_phi = 1 ;
138 
139  for (int k=0 ; k< borne_phi ; k++)
140  // On evite le coefficient de sin(0*phi)
141  if (k==1) {
142  xci += nr*nt ;
143  xco += nr*nt ;
144  }
145  else {
146  for (int j=0 ; j<nt ; j++) {
147 
148  double som ;
149 
150  xco[nr-1] = 0 ;
151  som = 2*(nr-1) * xci[nr-1] ;
152  xco[nr-2] = som ;
153  for (int i = nr-4 ; i >= 0 ; i -= 2 ) {
154  som += 2*(i+1) * xci[i+1] ;
155  xco[i] = som ;
156  } // Fin de la premiere boucle sur r
157  som = 2*(nr-2) * xci[nr-2] ;
158  xco[nr-3] = som ;
159  for (int i = nr-5 ; i >= 0 ; i -= 2 ) {
160  som += 2*(i+1) * xci[i+1] ;
161  xco[i] = som ;
162  } // Fin de la deuxieme boucle sur r
163  xco[0] *= .5 ;
164 
165  xci += nr ;
166  xco += nr ;
167  } // Fin de la boucle sur theta
168  } // Fin de la boucle sur phi
169  }
170  // On remet les choses la ou il faut
171  delete [] tb->t ;
172  tb->t = xo ;
173 
174  // base de developpement
175  // inchangee
176 }
177 
178 // cas R_CHEBU
179 //------------
180 void _dsdx_r_chebu(Tbl *tb, int & )
181 {
182 
183  // Peut-etre rien a faire ?
184  if (tb->get_etat() == ETATZERO) {
185  return ;
186  }
187 
188  // Protection
189  assert(tb->get_etat() == ETATQCQ) ;
190 
191  // Pour le confort
192  int nr = (tb->dim).dim[0] ; // Nombre
193  int nt = (tb->dim).dim[1] ; // de points
194  int np = (tb->dim).dim[2] ; // physiques REELS
195  np = np - 2 ; // Nombre de points physiques
196 
197  // pt. sur le tableau de double resultat
198  double* xo = new double[(tb->dim).taille] ;
199 
200  // Initialisation a zero :
201  for (int i=0; i<(tb->dim).taille; i++) {
202  xo[i] = 0 ;
203  }
204 
205  // On y va...
206  double* xi = tb->t ;
207  double* xci = xi ; // Pointeurs
208  double* xco = xo ; // courants
209 
210  int borne_phi = np + 1 ;
211  if (np == 1) borne_phi = 1 ;
212 
213  for (int k=0 ; k< borne_phi ; k++)
214 
215  // On evite le coefficient de sin(0*phi)
216  if (k==1) {
217  xci += nr*nt ;
218  xco += nr*nt ;
219  }
220 
221  else {
222  for (int j=0 ; j<nt ; j++) {
223 
224  double som ;
225 
226  xco[nr-1] = 0 ;
227  som = 2*(nr-1) * xci[nr-1] ;
228  xco[nr-2] = som ;
229  for (int i = nr-4 ; i >= 0 ; i -= 2 ) {
230  som += 2*(i+1) * xci[i+1] ;
231  xco[i] = som ;
232  } // Fin de la premiere boucle sur r
233  som = 2*(nr-2) * xci[nr-2] ;
234  xco[nr-3] = som ;
235  for (int i = nr-5 ; i >= 0 ; i -= 2 ) {
236  som += 2*(i+1) * xci[i+1] ;
237  xco[i] = som ;
238  } // Fin de la deuxieme boucle sur r
239  xco[0] *= .5 ;
240 
241  xci += nr ;
242  xco += nr ;
243  } // Fin de la boucle sur theta
244  } // Fin de la boucle sur phi
245 
246  // On remet les choses la ou il faut
247  delete [] tb->t ;
248  tb->t = xo ;
249 
250  // base de developpement
251  // inchangee
252 }
253 
254 // cas R_CHEBP
255 //------------
256 void _dsdx_r_chebp(Tbl *tb, int & b)
257 {
258 
259  // Peut-etre rien a faire ?
260  if (tb->get_etat() == ETATZERO) {
261  int base_t = b & MSQ_T ;
262  int base_p = b & MSQ_P ;
263  b = base_p | base_t | R_CHEBI ;
264  return ;
265  }
266 
267  // Protection
268  assert(tb->get_etat() == ETATQCQ) ;
269 
270  // Pour le confort
271  int nr = (tb->dim).dim[0] ; // Nombre
272  int nt = (tb->dim).dim[1] ; // de points
273  int np = (tb->dim).dim[2] ; // physiques REELS
274  np = np - 2 ; // Nombre de points physiques
275 
276  // pt. sur le tableau de double resultat
277  double* xo = new double[(tb->dim).taille] ;
278 
279  // Initialisation a zero :
280  for (int i=0; i<(tb->dim).taille; i++) {
281  xo[i] = 0 ;
282  }
283 
284  // On y va...
285  double* xi = tb->t ;
286  double* xci = xi ; // Pointeurs
287  double* xco = xo ; // courants
288 
289  int borne_phi = np + 1 ;
290  if (np == 1) borne_phi = 1 ;
291 
292  for (int k=0 ; k< borne_phi ; k++)
293 
294  // On evite le coefficient de sin(0*phi)
295 
296  if (k==1) {
297  xco += nr*nt ;
298  xci += nr*nt ;
299  }
300 
301 
302  else {
303  for (int j=0 ; j<nt ; j++) {
304 
305  double som ;
306 
307  xco[nr-1] = 0 ;
308  som = 4*(nr-1) * xci[nr-1] ;
309  xco[nr-2] = som ;
310  for (int i = nr-3 ; i >= 0 ; i-- ) {
311  som += 4*(i+1) * xci[i+1] ;
312  xco[i] = som ;
313  } // Fin de la boucle sur r
314 
315  xci += nr ;
316  xco += nr ;
317  } // Fin de la boucle sur theta
318  } // Fin de la boucle sur phi
319 
320  // On remet les choses la ou il faut
321  delete [] tb->t ;
322  tb->t = xo ;
323 
324  // base de developpement
325  // pair -> impair
326  int base_t = b & MSQ_T ;
327  int base_p = b & MSQ_P ;
328  b = base_p | base_t | R_CHEBI ;
329 }
330 
331 // cas R_CHEBI
332 //------------
333 void _dsdx_r_chebi(Tbl *tb, int & b)
334 {
335 
336  // Peut-etre rien a faire ?
337  if (tb->get_etat() == ETATZERO) {
338  int base_t = b & MSQ_T ;
339  int base_p = b & MSQ_P ;
340  b = base_p | base_t | R_CHEBP ;
341  return ;
342  }
343 
344  // Protection
345  assert(tb->get_etat() == ETATQCQ) ;
346 
347  // Pour le confort
348  int nr = (tb->dim).dim[0] ; // Nombre
349  int nt = (tb->dim).dim[1] ; // de points
350  int np = (tb->dim).dim[2] ; // physiques REELS
351  np = np - 2 ; // Nombre de points physiques
352 
353  // pt. sur le tableau de double resultat
354  double* xo = new double[(tb->dim).taille] ;
355 
356  // Initialisation a zero :
357  for (int i=0; i<(tb->dim).taille; i++) {
358  xo[i] = 0 ;
359  }
360 
361  // On y va...
362  double* xi = tb->t ;
363  double* xci = xi ; // Pointeurs
364  double* xco = xo ; // courants
365 
366  int borne_phi = np + 1 ;
367  if (np == 1) borne_phi = 1 ;
368 
369  for (int k=0 ; k< borne_phi ; k++)
370  // On evite le coefficient de sin(0*phi)
371  if (k==1) {
372  xci += nr*nt ;
373  xco += nr*nt ;
374  }
375 
376  else {
377  for (int j=0 ; j<nt ; j++) {
378 
379  double som ;
380 
381  xco[nr-1] = 0 ;
382  som = 2*(2*nr-3) * xci[nr-2] ;
383  xco[nr-2] = som ;
384  for (int i = nr-3 ; i >= 0 ; i-- ) {
385  som += 2*(2*i+1) * xci[i] ;
386  xco[i] = som ;
387  } // Fin de la boucle sur r
388  xco[0] *= .5 ;
389 
390  xci += nr ;
391  xco += nr ;
392  } // Fin de la boucle sur theta
393  } // Fin de la boucle sur phi
394 
395  // On remet les choses la ou il faut
396  delete [] tb->t ;
397  tb->t = xo ;
398 
399  // base de developpement
400  // impair -> pair
401  int base_t = b & MSQ_T ;
402  int base_p = b & MSQ_P ;
403  b = base_p | base_t | R_CHEBP ;
404 }
405 
406 // cas R_CHEBPIM_P
407 //----------------
408 void _dsdx_r_chebpim_p(Tbl *tb, int & b)
409 {
410 
411  // Peut-etre rien a faire ?
412  if (tb->get_etat() == ETATZERO) {
413  int base_t = b & MSQ_T ;
414  int base_p = b & MSQ_P ;
415  b = base_p | base_t | R_CHEBPIM_I ;
416  return ;
417  }
418 
419  // Pour le confort
420  int nr = (tb->dim).dim[0] ; // Nombre
421  int nt = (tb->dim).dim[1] ; // de points
422  int np = (tb->dim).dim[2] ; // physiques REELS
423  np = np - 2 ; // Nombre de points physiques
424 
425  // pt. sur le tableau de double resultat
426  double* xo = new double[(tb->dim).taille] ;
427 
428  // Initialisation a zero :
429  for (int i=0; i<(tb->dim).taille; i++) {
430  xo[i] = 0 ;
431  }
432 
433  // On y va...
434  double* xi = tb->t ;
435  double* xci ; // Pointeurs
436  double* xco ; // courants
437 
438  int auxiliaire ;
439  // Partie paire
440  xci = xi ;
441  xco = xo ;
442  for (int k=0 ; k<np+1 ; k += 4) {
443  auxiliaire = (k==np) ? 1 : 2 ; // evite de prendre le dernier coef
444  for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
445 
446 
447  // On evite le sin (0*phi)
448 
449  if ((k==0) && (kmod == 1)) {
450  xco += nr*nt ;
451  xci += nr*nt ;
452  }
453 
454  else
455 
456  for (int j=0 ; j<nt ; j++) {
457 
458  double som ;
459 
460  xco[nr-1] = 0 ;
461  som = 4*(nr-1) * xci[nr-1] ;
462  xco[nr-2] = som ;
463  for (int i = nr-3 ; i >= 0 ; i-- ) {
464  som += 4*(i+1) * xci[i+1] ;
465  xco[i] = som ;
466  } // Fin de la boucle sur r
467 
468  xci += nr ; // au
469  xco += nr ; // suivant
470  } // Fin de la boucle sur theta
471  } // Fin de la boucle sur kmod
472  xci += 2*nr*nt ; // saute
473  xco += 2*nr*nt ; // 2 phis
474  } // Fin de la boucle sur phi
475 
476  // Partie impaire
477  xci = xi + 2*nr*nt ;
478  xco = xo + 2*nr*nt ;
479  for (int k=2 ; k<np+1 ; k += 4) {
480  auxiliaire = (k==np) ? 1 : 2 ; // evite de prendre le dernier coef
481  for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
482  for (int j=0 ; j<nt ; j++) {
483 
484  double som ;
485 
486  xco[nr-1] = 0 ;
487  som = 2*(2*nr-3) * xci[nr-2] ;
488  xco[nr-2] = som ;
489  for (int i = nr-3 ; i >= 0 ; i-- ) {
490  som += 2*(2*i+1) * xci[i] ;
491  xco[i] = som ;
492  } // Fin de la boucle sur r
493  xco[0] *= .5 ;
494 
495  xci += nr ;
496  xco += nr ;
497  } // Fin de la boucle sur theta
498  } // Fin de la boucle sur kmod
499  xci += 2*nr*nt ;
500  xco += 2*nr*nt ;
501  } // Fin de la boucle sur phi
502 
503  // On remet les choses la ou il faut
504  delete [] tb->t ;
505  tb->t = xo ;
506 
507  // base de developpement
508  // (pair,impair) -> (impair,pair)
509  int base_t = b & MSQ_T ;
510  int base_p = b & MSQ_P ;
511  b = base_p | base_t | R_CHEBPIM_I ;
512 }
513 
514 // cas R_CHEBPIM_I
515 //----------------
516 void _dsdx_r_chebpim_i(Tbl *tb, int & b)
517 {
518 
519  // Peut-etre rien a faire ?
520  if (tb->get_etat() == ETATZERO) {
521  int base_t = b & MSQ_T ;
522  int base_p = b & MSQ_P ;
523  b = base_p | base_t | R_CHEBPIM_P ;
524  return ;
525  }
526 
527  // Protection
528  assert(tb->get_etat() == ETATQCQ) ;
529 
530  // Pour le confort
531  // Pour le confort
532  int nr = (tb->dim).dim[0] ; // Nombre
533  int nt = (tb->dim).dim[1] ; // de points
534  int np = (tb->dim).dim[2] ; // physiques REELS
535  np = np - 2 ; // Nombre de points physiques
536 
537  // pt. sur le tableau de double resultat
538  double* xo = new double[(tb->dim).taille] ;
539 
540  // Initialisation a zero :
541  for (int i=0; i<(tb->dim).taille; i++) {
542  xo[i] = 0 ;
543  }
544 
545  // On y va...
546  double* xi = tb->t ;
547  double* xci ; // Pointeurs
548  double* xco ; // courants
549  int auxiliaire ;
550 
551  // Partie impaire
552  xci = xi ;
553  xco = xo ;
554  for (int k=0 ; k<np+1 ; k += 4) {
555  auxiliaire = (k==np) ? 1 : 2 ; // evite de prendre le dernier coef
556  for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
557 
558 
559  // On evite le sin (0*phi)
560 
561  if ((k==0) && (kmod == 1)) {
562  xco += nr*nt ;
563  xci += nr*nt ;
564  }
565 
566  else
567 
568  for (int j=0 ; j<nt ; j++) {
569 
570  double som ;
571 
572  xco[nr-1] = 0 ;
573  som = 2*(2*nr-3) * xci[nr-2] ;
574  xco[nr-2] = som ;
575  for (int i = nr-3 ; i >= 0 ; i-- ) {
576  som += 2*(2*i+1) * xci[i] ;
577  xco[i] = som ;
578  } // Fin de la boucle sur r
579  xco[0] *= .5 ;
580 
581  xci += nr ;
582  xco += nr ;
583  } // Fin de la boucle sur theta
584  } // Fin de la boucle sur kmod
585  xci += 2*nr*nt ;
586  xco += 2*nr*nt ;
587  } // Fin de la boucle sur phi
588 
589  // Partie paire
590  xci = xi + 2*nr*nt ;
591  xco = xo + 2*nr*nt ;
592  for (int k=2 ; k<np+1 ; k += 4) {
593  auxiliaire = (k==np) ? 1 : 2 ; // evite de prendre le dernier coef
594  for (int kmod=0 ; kmod<auxiliaire ; kmod++) {
595  for (int j=0 ; j<nt ; j++) {
596 
597  double som ;
598 
599  xco[nr-1] = 0 ;
600  som = 4*(nr-1) * xci[nr-1] ;
601  xco[nr-2] = som ;
602  for (int i = nr-3 ; i >= 0 ; i-- ) {
603  som += 4*(i+1) * xci[i+1] ;
604  xco[i] = som ;
605  } // Fin de la boucle sur r
606 
607  xci += nr ;
608  xco += nr ;
609  } // Fin de la boucle sur theta
610  } // Fin de la boucle sur kmod
611  xci += 2*nr*nt ;
612  xco += 2*nr*nt ;
613  } // Fin de la boucle sur phi
614 
615  // On remet les choses la ou il faut
616  delete [] tb->t ;
617  tb->t = xo ;
618 
619  // base de developpement
620  // (impair,pair) -> (pair,impair)
621  int base_t = b & MSQ_T ;
622  int base_p = b & MSQ_P ;
623  b = base_p | base_t | R_CHEBPIM_P ;
624 }
625 
626 // cas R_CHEBPI_P
627 //------------
628 void _dsdx_r_chebpi_p(Tbl *tb, int & b)
629 {
630 
631  // Peut-etre rien a faire ?
632  if (tb->get_etat() == ETATZERO) {
633  int base_t = b & MSQ_T ;
634  int base_p = b & MSQ_P ;
635  b = base_p | base_t | R_CHEBPI_I ;
636  return ;
637  }
638 
639  // Protection
640  assert(tb->get_etat() == ETATQCQ) ;
641 
642  // Pour le confort
643  int nr = (tb->dim).dim[0] ; // Nombre
644  int nt = (tb->dim).dim[1] ; // de points
645  int np = (tb->dim).dim[2] ; // physiques REELS
646  np = np - 2 ; // Nombre de points physiques
647 
648  // pt. sur le tableau de double resultat
649  double* xo = new double[(tb->dim).taille] ;
650 
651  // Initialisation a zero :
652  for (int i=0; i<(tb->dim).taille; i++) {
653  xo[i] = 0 ;
654  }
655 
656  // On y va...
657  double* xi = tb->t ;
658  double* xci = xi ; // Pointeurs
659  double* xco = xo ; // courants
660 
661  int borne_phi = np + 1 ;
662  if (np == 1) borne_phi = 1 ;
663 
664  for (int k=0 ; k< borne_phi ; k++)
665 
666  // On evite le coefficient de sin(0*phi)
667 
668  if (k==1) {
669  xco += nr*nt ;
670  xci += nr*nt ;
671  }
672 
673 
674  else {
675  for (int j=0 ; j<nt ; j++) {
676  int l = j%2 ;
677  double som ;
678 
679  if(l==0){
680  xco[nr-1] = 0 ;
681  som = 4*(nr-1) * xci[nr-1] ;
682  xco[nr-2] = som ;
683  for (int i = nr-3 ; i >= 0 ; i-- ) {
684  som += 4*(i+1) * xci[i+1] ;
685  xco[i] = som ;
686  } // Fin de la boucle sur r
687  } else {
688  xco[nr-1] = 0 ;
689  som = 2*(2*nr-3) * xci[nr-2] ;
690  xco[nr-2] = som ;
691  for (int i = nr-3 ; i >= 0 ; i-- ) {
692  som += 2*(2*i+1) * xci[i] ;
693  xco[i] = som ;
694  } // Fin de la boucle sur r
695  xco[0] *= .5 ;
696  }
697  xci += nr ;
698  xco += nr ;
699  } // Fin de la boucle sur theta
700  } // Fin de la boucle sur phi
701 
702  // On remet les choses la ou il faut
703  delete [] tb->t ;
704  tb->t = xo ;
705 
706  // base de developpement
707  // pair -> impair
708  int base_t = b & MSQ_T ;
709  int base_p = b & MSQ_P ;
710  b = base_p | base_t | R_CHEBPI_I ;
711 }
712 
713 // cas R_CHEBPI_I
714 //------------
715 void _dsdx_r_chebpi_i(Tbl *tb, int & b)
716 {
717 
718  // Peut-etre rien a faire ?
719  if (tb->get_etat() == ETATZERO) {
720  int base_t = b & MSQ_T ;
721  int base_p = b & MSQ_P ;
722  b = base_p | base_t | R_CHEBPI_P ;
723  return ;
724  }
725 
726  // Protection
727  assert(tb->get_etat() == ETATQCQ) ;
728 
729  // Pour le confort
730  int nr = (tb->dim).dim[0] ; // Nombre
731  int nt = (tb->dim).dim[1] ; // de points
732  int np = (tb->dim).dim[2] ; // physiques REELS
733  np = np - 2 ; // Nombre de points physiques
734 
735  // pt. sur le tableau de double resultat
736  double* xo = new double[(tb->dim).taille] ;
737 
738  // Initialisation a zero :
739  for (int i=0; i<(tb->dim).taille; i++) {
740  xo[i] = 0 ;
741  }
742 
743  // On y va...
744  double* xi = tb->t ;
745  double* xci = xi ; // Pointeurs
746  double* xco = xo ; // courants
747 
748  int borne_phi = np + 1 ;
749  if (np == 1) borne_phi = 1 ;
750 
751  for (int k=0 ; k< borne_phi ; k++)
752  // On evite le coefficient de sin(0*phi)
753  if (k==1) {
754  xci += nr*nt ;
755  xco += nr*nt ;
756  }
757 
758  else {
759  for (int j=0 ; j<nt ; j++) {
760  int l = j%2 ;
761  double som ;
762 
763  if(l==1){
764  xco[nr-1] = 0 ;
765  som = 4*(nr-1) * xci[nr-1] ;
766  xco[nr-2] = som ;
767  for (int i = nr-3 ; i >= 0 ; i-- ) {
768  som += 4*(i+1) * xci[i+1] ;
769  xco[i] = som ;
770  } // Fin de la boucle sur r
771  } else {
772  xco[nr-1] = 0 ;
773  som = 2*(2*nr-3) * xci[nr-2] ;
774  xco[nr-2] = som ;
775  for (int i = nr-3 ; i >= 0 ; i-- ) {
776  som += 2*(2*i+1) * xci[i] ;
777  xco[i] = som ;
778  } // Fin de la boucle sur r
779  xco[0] *= .5 ;
780  }
781  xci += nr ;
782  xco += nr ;
783  } // Fin de la boucle sur theta
784  } // Fin de la boucle sur phi
785 
786  // On remet les choses la ou il faut
787  delete [] tb->t ;
788  tb->t = xo ;
789 
790  // base de developpement
791  // impair -> pair
792  int base_t = b & MSQ_T ;
793  int base_p = b & MSQ_P ;
794  b = base_p | base_t | R_CHEBPI_P ;
795 }
796 
797 
798 // cas R_LEG
799 //-----------
800 void _dsdx_r_leg(Tbl *tb, int & )
801 {
802 
803  // Peut-etre rien a faire ?
804  if (tb->get_etat() == ETATZERO) {
805  return ;
806  }
807 
808  // Protection
809  assert(tb->get_etat() == ETATQCQ) ;
810 
811  // Pour le confort
812  int nr = (tb->dim).dim[0] ; // Nombre
813  int nt = (tb->dim).dim[1] ; // de points
814  int np = (tb->dim).dim[2] ; // physiques REELS
815  np = np - 2 ; // Nombre de points physiques
816 
817  // pt. sur le tableau de double resultat
818  double* xo = new double[(tb->dim).taille] ;
819 
820  // Initialisation a zero :
821  for (int i=0; i<(tb->dim).taille; i++) {
822  xo[i] = 0 ;
823  }
824  if (nr > 1) { // If not an angular grid...
825  // On y va...
826  double* xi = tb->t ;
827  double* xci = xi ; // Pointeurs
828  double* xco = xo ; // courants
829 
830  int borne_phi = np + 1 ;
831  if (np == 1) borne_phi = 1 ;
832 
833  for (int k=0 ; k< borne_phi ; k++)
834  // On evite le coefficient de sin(0*phi)
835  if (k==1) {
836  xci += nr*nt ;
837  xco += nr*nt ;
838  }
839  else {
840  for (int j=0 ; j<nt ; j++) {
841 
842  double som ;
843 
844  xco[nr-1] = 0 ;
845  som = xci[nr-1] ;
846  xco[nr-2] = double(2*nr-3)*som ;
847  for (int i = nr-4 ; i >= 0 ; i -= 2 ) {
848  som += xci[i+1] ;
849  xco[i] = double(2*i+1)*som ;
850  } // Fin de la premiere boucle sur r
851  som = xci[nr-2] ;
852  if (nr > 2) xco[nr-3] = double(2*nr-5)*som ;
853  for (int i = nr-5 ; i >= 0 ; i -= 2 ) {
854  som += xci[i+1] ;
855  xco[i] = double(2*i+1) * som ;
856  } // Fin de la deuxieme boucle sur r
857 
858  xci += nr ;
859  xco += nr ;
860  } // Fin de la boucle sur theta
861  } // Fin de la boucle sur phi
862  }
863  // On remet les choses la ou il faut
864  delete [] tb->t ;
865  tb->t = xo ;
866 
867  // base de developpement
868  // inchangee
869 }
870 
871 // cas R_LEGP
872 //------------
873 void _dsdx_r_legp(Tbl *tb, int & b)
874 {
875 
876  // Peut-etre rien a faire ?
877  if (tb->get_etat() == ETATZERO) {
878  int base_t = b & MSQ_T ;
879  int base_p = b & MSQ_P ;
880  b = base_p | base_t | R_LEGI ;
881  return ;
882  }
883 
884  // Protection
885  assert(tb->get_etat() == ETATQCQ) ;
886 
887  // Pour le confort
888  int nr = (tb->dim).dim[0] ; // Nombre
889  int nt = (tb->dim).dim[1] ; // de points
890  int np = (tb->dim).dim[2] ; // physiques REELS
891  np = np - 2 ; // Nombre de points physiques
892 
893  // pt. sur le tableau de double resultat
894  double* xo = new double[(tb->dim).taille] ;
895 
896  // Initialisation a zero :
897  for (int i=0; i<(tb->dim).taille; i++) {
898  xo[i] = 0 ;
899  }
900 
901  // On y va...
902  double* xi = tb->t ;
903  double* xci = xi ; // Pointeurs
904  double* xco = xo ; // courants
905 
906  int borne_phi = np + 1 ;
907  if (np == 1) borne_phi = 1 ;
908 
909  for (int k=0 ; k< borne_phi ; k++)
910 
911  // On evite le coefficient de sin(0*phi)
912 
913  if (k==1) {
914  xco += nr*nt ;
915  xci += nr*nt ;
916  }
917 
918 
919  else {
920  for (int j=0 ; j<nt ; j++) {
921 
922  double som ;
923 
924  xco[nr-1] = 0 ;
925  som = xci[nr-1] ;
926  if (nr > 1) xco[nr-2] = double(4*nr-5) * som ;
927  for (int i = nr-3 ; i >= 0 ; i-- ) {
928  som += xci[i+1] ;
929  xco[i] = double(4*i+3) * som ;
930  } // Fin de la boucle sur r
931 
932  xci += nr ;
933  xco += nr ;
934  } // Fin de la boucle sur theta
935  } // Fin de la boucle sur phi
936 
937  // On remet les choses la ou il faut
938  delete [] tb->t ;
939  tb->t = xo ;
940 
941  // base de developpement
942  // pair -> impair
943  int base_t = b & MSQ_T ;
944  int base_p = b & MSQ_P ;
945  b = base_p | base_t | R_LEGI ;
946 }
947 
948 // cas R_LEGI
949 //------------
950 void _dsdx_r_legi(Tbl *tb, int & b)
951 {
952 
953  // Peut-etre rien a faire ?
954  if (tb->get_etat() == ETATZERO) {
955  int base_t = b & MSQ_T ;
956  int base_p = b & MSQ_P ;
957  b = base_p | base_t | R_LEGP ;
958  return ;
959  }
960 
961  // Protection
962  assert(tb->get_etat() == ETATQCQ) ;
963 
964  // Pour le confort
965  int nr = (tb->dim).dim[0] ; // Nombre
966  int nt = (tb->dim).dim[1] ; // de points
967  int np = (tb->dim).dim[2] ; // physiques REELS
968  np = np - 2 ; // Nombre de points physiques
969 
970  // pt. sur le tableau de double resultat
971  double* xo = new double[(tb->dim).taille] ;
972 
973  // Initialisation a zero :
974  for (int i=0; i<(tb->dim).taille; i++) {
975  xo[i] = 0 ;
976  }
977 
978  // On y va...
979  double* xi = tb->t ;
980  double* xci = xi ; // Pointeurs
981  double* xco = xo ; // courants
982 
983  int borne_phi = np + 1 ;
984  if (np == 1) borne_phi = 1 ;
985 
986  for (int k=0 ; k< borne_phi ; k++)
987  // On evite le coefficient de sin(0*phi)
988  if (k==1) {
989  xci += nr*nt ;
990  xco += nr*nt ;
991  }
992 
993  else {
994  for (int j=0 ; j<nt ; j++) {
995 
996  double som ;
997 
998  xco[nr-1] = 0 ;
999  som = xci[nr-2] ;
1000  if (nr > 1) xco[nr-2] = double(4*nr - 7) * som ;
1001  for (int i = nr-3 ; i >= 0 ; i-- ) {
1002  som += xci[i] ;
1003  xco[i] = double(4*i+1) * som ;
1004  } // Fin de la boucle sur r
1005 
1006  xci += nr ;
1007  xco += nr ;
1008  } // Fin de la boucle sur theta
1009  } // Fin de la boucle sur phi
1010 
1011  // On remet les choses la ou il faut
1012  delete [] tb->t ;
1013  tb->t = xo ;
1014 
1015  // base de developpement
1016  // impair -> pair
1017  int base_t = b & MSQ_T ;
1018  int base_p = b & MSQ_P ;
1019  b = base_p | base_t | R_LEGP ;
1020 }
1021 
1022 // cas R_JACO02
1023 //-----------
1024 void _dsdx_r_jaco02(Tbl *tb, int & )
1025 {
1026 
1027  // Peut-etre rien a faire ?
1028  if (tb->get_etat() == ETATZERO) {
1029  return ;
1030  }
1031 
1032  // Protection
1033  assert(tb->get_etat() == ETATQCQ) ;
1034 
1035  // Pour le confort
1036  int nr = (tb->dim).dim[0] ; // Nombre
1037  int nt = (tb->dim).dim[1] ; // de points
1038  int np = (tb->dim).dim[2] ; // physiques REELS
1039  np = np - 2 ; // Nombre de points physiques
1040 
1041  // pt. sur le tableau de double resultat
1042  double* xo = new double[(tb->dim).taille] ;
1043 
1044  // Initialisation a zero :
1045  for (int i=0; i<(tb->dim).taille; i++) {
1046  xo[i] = 0 ;
1047  }
1048  if (nr > 2) { // If not an angular grid...
1049  // On y va...
1050  double* xi = tb->t ;
1051  double* xci = xi ; // Pointeurs
1052  double* xco = xo ; // courants
1053 
1054  int borne_phi = np + 1 ;
1055  if (np == 1) borne_phi = 1 ;
1056 
1057  for (int k=0 ; k< borne_phi ; k++)
1058  // On evite le coefficient de sin(0*phi)
1059  if (k==1) {
1060  xci += nr*nt ;
1061  xco += nr*nt ;
1062  }
1063  else {
1064  for (int j=0 ; j<nt ; j++) {
1065 
1066  double som ;
1067  xco[nr-1] = 0 ;
1068 
1069  for (int i = 0 ; i < nr-1 ; i++ ) {
1070 
1071  som = 0 ;
1072  for (int m = i+1 ; m < nr ; m++ ) {
1073  int signe = ((m-i)%2 == 0 ? 1 : -1) ;
1074  som += (1-signe*(i+1)*(i+2)/double((m+1)*(m+2)))* xci[m] ;
1075  } // Fin de la boucle annexe
1076  xco[i] = (i+1.5)*som ;
1077  }// Fin de la boucle sur r
1078 
1079  xci += nr ;
1080  xco += nr ;
1081  } // Fin de la boucle sur theta
1082  } // Fin de la boucle sur phi
1083  }
1084  // On remet les choses la ou il faut
1085  delete [] tb->t ;
1086  tb->t = xo ;
1087 
1088  // base de developpement
1089  // inchangee
1090 }
1091 }
R_CHEBPI_I
#define R_CHEBPI_I
Cheb. pair-impair suivant l impair pour l=0.
Definition: type_parite.h:174
Lorene
Lorene prototypes.
Definition: app_hor.h:67
MSQ_P
#define MSQ_P
Extraction de l&#39;info sur Phi.
Definition: type_parite.h:156
R_LEGP
#define R_LEGP
base de Legendre paire (rare) seulement
Definition: type_parite.h:184
R_LEGI
#define R_LEGI
base de Legendre impaire (rare) seulement
Definition: type_parite.h:186
R_CHEBI
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
MSQ_T
#define MSQ_T
Extraction de l&#39;info sur Theta.
Definition: type_parite.h:154
R_CHEBPIM_I
#define R_CHEBPIM_I
Cheb. pair-impair suivant m, impair pour m=0.
Definition: type_parite.h:178
R_CHEBPIM_P
#define R_CHEBPIM_P
Cheb. pair-impair suivant m, pair pour m=0.
Definition: type_parite.h:176
R_CHEBPI_P
#define R_CHEBPI_P
Cheb. pair-impair suivant l pair pour l=0.
Definition: type_parite.h:172