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//--------------------------------
97namespace Lorene {
98void _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//-----------
106void _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//------------
180void _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//------------
256void _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//------------
333void _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//----------------
408void _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//----------------
516void _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//------------
628void _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//------------
715void _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//-----------
800void _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//------------
873void _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//------------
950void _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//-----------
1024void _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}
Lorene::Tbl
Basic array class.
Definition tbl.h:161
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_CHEBPIM_I
#define R_CHEBPIM_I
Cheb. pair-impair suivant m, impair pour m=0.
Definition type_parite.h:178
R_CHEBPI_I
#define R_CHEBPI_I
Cheb. pair-impair suivant l impair pour l=0.
Definition type_parite.h:174
R_CHEBPIM_P
#define R_CHEBPIM_P
Cheb. pair-impair suivant m, pair pour m=0.
Definition type_parite.h:176
MSQ_T
#define MSQ_T
Extraction de l'info sur Theta.
Definition type_parite.h:154
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition type_parite.h:168
MSQ_P
#define MSQ_P
Extraction de l'info sur Phi.
Definition type_parite.h:156
R_CHEBPI_P
#define R_CHEBPI_P
Cheb. pair-impair suivant l pair pour l=0.
Definition type_parite.h:172
Lorene
Lorene prototypes.
Definition app_hor.h:67