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