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
25char 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//--------------------------------
94namespace Lorene {
95void _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//-----------
103void _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//------------
177void _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//------------
253void _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//------------
330void _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//----------------
405void _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//----------------
513void _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//------------
625void _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//------------
712void _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//-----------
797void _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//------------
870void _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//------------
947void _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//-----------
1021void _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}
#define R_LEGP
base de Legendre paire (rare) seulement
#define R_LEGI
base de Legendre impaire (rare) seulement
#define R_CHEBI
base de Cheb. impaire (rare) seulement
#define R_CHEBPIM_I
Cheb. pair-impair suivant m, impair pour m=0.
#define R_CHEBPI_I
Cheb. pair-impair suivant l impair pour l=0.
#define R_CHEBPIM_P
Cheb. pair-impair suivant m, pair pour m=0.
#define MSQ_T
Extraction de l'info sur Theta.
#define R_CHEBP
base de Cheb. paire (rare) seulement
#define MSQ_P
Extraction de l'info sur Phi.
#define R_CHEBPI_P
Cheb. pair-impair suivant l pair pour l=0.
Lorene prototypes.
Definition app_hor.h:64