LORENE
op_dsdtet.C
1/*
2 * Copyright (c) 1999-2000 Jean-Alain Marck
3 * Copyright (c) 1999-2001 Eric Gourgoulhon
4 *
5 * This file is part of LORENE.
6 *
7 * LORENE is free software; you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation; either version 2 of the License, or
10 * (at your option) any later version.
11 *
12 * LORENE is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with LORENE; if not, write to the Free Software
19 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 *
21 */
22
23
24
25
26/*
27 * Ensemble des routines de base pour la derivation par rapport a theta
28 * (Utilisation interne)
29 *
30 * void _dsdtet_XXXX(Tbl * t, int & b)
31 * t pointeur sur le Tbl d'entree-sortie
32 * b base spectrale
33 *
34 */
35
36/*
37 * $Id: op_dsdtet.C,v 1.7 2016/12/05 16:18:07 j_novak Exp $
38 * $Log: op_dsdtet.C,v $
39 * Revision 1.7 2016/12/05 16:18:07 j_novak
40 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
41 *
42 * Revision 1.6 2014/10/13 08:53:25 j_novak
43 * Lorene classes and functions now belong to the namespace Lorene.
44 *
45 * Revision 1.5 2009/10/08 16:22:19 j_novak
46 * Addition of new bases T_COS and T_SIN.
47 *
48 * Revision 1.4 2006/03/10 12:45:38 j_novak
49 * Use of C++-style cast.
50 *
51 * Revision 1.3 2004/11/23 15:16:01 m_forot
52 *
53 * Added the bases for the cases without any equatorial symmetry
54 * (T_COSSIN_C, T_COSSIN_S, T_LEG, R_CHEBPI_P, R_CHEBPI_I).
55 *
56 * Revision 1.2 2003/12/19 16:21:48 j_novak
57 * Shadow hunt
58 *
59 * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
60 * LORENE
61 *
62 * Revision 2.6 2000/10/04 11:50:44 eric
63 * Ajout d' abort() dans le cas non prevu.
64 *
65 * Revision 2.5 2000/02/24 16:41:17 eric
66 * Initialisation a zero du tableau xo avant le calcul.
67 *
68 * Revision 2.4 1999/11/15 16:38:13 eric
69 * Tbl::dim est desormais un Dim_tbl et non plus un Dim_tbl*.
70 *
71 * Revision 2.3 1999/03/01 15:07:01 eric
72 * *** empty log message ***
73 *
74 * Revision 2.2 1999/02/23 11:04:52 hyc
75 * *** empty log message ***
76 *
77 * Revision 2.1 1999/02/22 17:12:06 hyc
78 * *** empty log message ***
79 *
80 * Revision 2.0 1999/02/22 15:51:02 hyc
81 * *** empty log message ***
82 *
83 *
84 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_dsdtet.C,v 1.7 2016/12/05 16:18:07 j_novak Exp $
85 *
86 */
87
88// Headers Lorene
89#include "tbl.h"
90
91
92// Routine pour les cas non prevus
93//--------------------------------
94namespace Lorene {
95void _dsdtet_pas_prevu(Tbl* , int & b) {
96 cout << "Unknown theta basis in Mtbl_cf::dsdt() !" << endl ;
97 cout << " basis: " << hex << b << endl ;
98 abort() ;
99}
100
101// cas T_COS
102//------------
103void _dsdtet_t_cos(Tbl* tb, int & b)
104{
105
106 // Peut-etre rien a faire ?
107 if (tb->get_etat() == ETATZERO) {
108 int base_r = b & MSQ_R ;
109 int base_p = b & MSQ_P ;
110 b = base_r | base_p | T_SIN ;
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 // Variables statiques
124 static double* cx = 0 ;
125 static int nt_pre =0 ;
126
127 // Test sur np pour initialisation eventuelle
128 {
129 if (nt > nt_pre) {
130 nt_pre = nt ;
131 cx = reinterpret_cast<double*>(realloc(cx, nt * sizeof(double))) ;
132 for (int i=0 ; i<nt ; i++) {
133 cx[i] = - double(i) ;
134 }
135 }
136 } // Fin de region critique
137
138 // pt. sur le tableau de double resultat
139 double* xo = new double[(tb->dim).taille] ;
140
141 // Initialisation a zero :
142 for (int i=0; i<(tb->dim).taille; i++) {
143 xo[i] = 0 ;
144 }
145
146 // On y va...
147 double* xi = tb->t ;
148 double* xci = xi ; // Pointeurs
149 double* xco = xo ; // courants
150
151 // k = 0
152 for (int j=0 ; j<nt ; j++) {
153 for (int i=0 ; i<nr ; i++ ) {
154 *xco = cx[j] * (*xci) ;
155 xci++ ;
156 xco++ ;
157 } // Fin de la boucle sur r
158 } // Fin de la boucle sur theta
159
160 // k = 1
161 xci += nr*nt ;
162 xco += nr*nt ;
163
164 // k >=2
165 for (int k=2 ; k<np+1 ; k++) {
166 for (int j=0 ; j<nt ; j++) {
167 for (int i=0 ; i<nr ; i++ ) {
168 *xco = cx[j] * (*xci) ;
169 xci++ ;
170 xco++ ;
171 } // Fin de la boucle sur r
172 } // Fin de la boucle sur theta
173 } // Fin de la boucle sur phi
174
175 // On remet les choses la ou il faut
176 delete [] tb->t ;
177 tb->t = xo ;
178
179 // base de developpement
180 int base_r = b & MSQ_R ;
181 int base_p = b & MSQ_P ;
182 b = base_r | base_p | T_SIN ;
183}
184
185// cas T_SIN
186//------------
187void _dsdtet_t_sin(Tbl* tb, int & b)
188{
189
190 // Peut-etre rien a faire ?
191 if (tb->get_etat() == ETATZERO) {
192 int base_r = b & MSQ_R ;
193 int base_p = b & MSQ_P ;
194 b = base_r | base_p | T_COS ;
195 return ;
196 }
197
198 // Protection
199 assert(tb->get_etat() == ETATQCQ) ;
200
201 // Pour le confort
202 int nr = (tb->dim).dim[0] ; // Nombre
203 int nt = (tb->dim).dim[1] ; // de points
204 int np = (tb->dim).dim[2] ; // physiques REELS
205 np = np - 2 ; // Nombre de points physiques
206
207 // Variables statiques
208 static double* cx = 0 ;
209 static int nt_pre = 0 ;
210
211 // Test sur np pour initialisation eventuelle
212 {
213 if (nt > nt_pre) {
214 nt_pre = nt ;
215 cx = reinterpret_cast<double*>(realloc(cx, nt * sizeof(double))) ;
216 for (int i=0 ; i<nt ; i++) {
217 cx[i] = double(i) ;
218 }
219 }
220 } // Fin de region critique
221
222 // pt. sur le tableau de double resultat
223 double* xo = new double[(tb->dim).taille] ;
224
225 // Initialisation a zero :
226 for (int i=0; i<(tb->dim).taille; i++) {
227 xo[i] = 0 ;
228 }
229
230 // On y va...
231 double* xi = tb->t ;
232 double* xci = xi ; // Pointeurs
233 double* xco = xo ; // courants
234
235 // k = 0
236 for (int j=0 ; j<nt ; j++) {
237 for (int i=0 ; i<nr ; i++ ) {
238 *xco = cx[j] * (*xci) ;
239 xci++ ;
240 xco++ ;
241 } // Fin de la boucle sur r
242 } // Fin de la boucle sur theta
243
244 // k = 1
245 xci += nr*nt ;
246 xco += nr*nt ;
247
248 // k >=2
249 for (int k=2 ; k<np+1 ; k++) {
250 for (int j=0 ; j<nt ; j++) {
251 for (int i=0 ; i<nr ; i++ ) {
252 *xco = cx[j] * (*xci) ;
253 xci++ ;
254 xco++ ;
255 } // Fin de la boucle sur r
256 } // Fin de la boucle sur theta
257 } // Fin de la boucle sur phi
258
259 // On remet les choses la ou il faut
260 delete [] tb->t ;
261 tb->t = xo ;
262
263 // base de developpement
264 int base_r = b & MSQ_R ;
265 int base_p = b & MSQ_P ;
266 b = base_r | base_p | T_COS ;
267}
268
269// cas T_COS_P
270//------------
271void _dsdtet_t_cos_p(Tbl* tb, int & b)
272{
273
274 // Peut-etre rien a faire ?
275 if (tb->get_etat() == ETATZERO) {
276 int base_r = b & MSQ_R ;
277 int base_p = b & MSQ_P ;
278 b = base_r | base_p | T_SIN_P ;
279 return ;
280 }
281
282 // Protection
283 assert(tb->get_etat() == ETATQCQ) ;
284
285 // Pour le confort
286 int nr = (tb->dim).dim[0] ; // Nombre
287 int nt = (tb->dim).dim[1] ; // de points
288 int np = (tb->dim).dim[2] ; // physiques REELS
289 np = np - 2 ; // Nombre de points physiques
290
291 // Variables statiques
292 static double* cx = 0 ;
293 static int nt_pre =0 ;
294
295 // Test sur np pour initialisation eventuelle
296 {
297 if (nt > nt_pre) {
298 nt_pre = nt ;
299 cx = reinterpret_cast<double*>(realloc(cx, nt * sizeof(double))) ;
300 for (int i=0 ; i<nt ; i++) {
301 cx[i] = - (2*i) ;
302 }
303 }
304 } // Fin de region critique
305
306 // pt. sur le tableau de double resultat
307 double* xo = new double[(tb->dim).taille] ;
308
309 // Initialisation a zero :
310 for (int i=0; i<(tb->dim).taille; i++) {
311 xo[i] = 0 ;
312 }
313
314 // On y va...
315 double* xi = tb->t ;
316 double* xci = xi ; // Pointeurs
317 double* xco = xo ; // courants
318
319 // k = 0
320 for (int j=0 ; j<nt ; j++) {
321 for (int i=0 ; i<nr ; i++ ) {
322 *xco = cx[j] * (*xci) ;
323 xci++ ;
324 xco++ ;
325 } // Fin de la boucle sur r
326 } // Fin de la boucle sur theta
327
328 // k = 1
329 xci += nr*nt ;
330 xco += nr*nt ;
331
332 // k >=2
333 for (int k=2 ; k<np+1 ; k++) {
334 for (int j=0 ; j<nt ; j++) {
335 for (int i=0 ; i<nr ; i++ ) {
336 *xco = cx[j] * (*xci) ;
337 xci++ ;
338 xco++ ;
339 } // Fin de la boucle sur r
340 } // Fin de la boucle sur theta
341 } // Fin de la boucle sur phi
342
343 // On remet les choses la ou il faut
344 delete [] tb->t ;
345 tb->t = xo ;
346
347 // base de developpement
348 int base_r = b & MSQ_R ;
349 int base_p = b & MSQ_P ;
350 b = base_r | base_p | T_SIN_P ;
351}
352
353// cas T_SIN_P
354//------------
355void _dsdtet_t_sin_p(Tbl* tb, int & b)
356{
357
358 // Peut-etre rien a faire ?
359 if (tb->get_etat() == ETATZERO) {
360 int base_r = b & MSQ_R ;
361 int base_p = b & MSQ_P ;
362 b = base_r | base_p | T_COS_P ;
363 return ;
364 }
365
366 // Protection
367 assert(tb->get_etat() == ETATQCQ) ;
368
369 // Pour le confort
370 int nr = (tb->dim).dim[0] ; // Nombre
371 int nt = (tb->dim).dim[1] ; // de points
372 int np = (tb->dim).dim[2] ; // physiques REELS
373 np = np - 2 ; // Nombre de points physiques
374
375 // Variables statiques
376 static double* cx = 0 ;
377 static int nt_pre = 0 ;
378
379 // Test sur np pour initialisation eventuelle
380 {
381 if (nt > nt_pre) {
382 nt_pre = nt ;
383 cx = reinterpret_cast<double*>(realloc(cx, nt * sizeof(double))) ;
384 for (int i=0 ; i<nt ; i++) {
385 cx[i] = (2*i) ;
386 }
387 }
388 } // Fin de region critique
389
390 // pt. sur le tableau de double resultat
391 double* xo = new double[(tb->dim).taille] ;
392
393 // Initialisation a zero :
394 for (int i=0; i<(tb->dim).taille; i++) {
395 xo[i] = 0 ;
396 }
397
398 // On y va...
399 double* xi = tb->t ;
400 double* xci = xi ; // Pointeurs
401 double* xco = xo ; // courants
402
403 // k = 0
404 for (int j=0 ; j<nt ; j++) {
405 for (int i=0 ; i<nr ; i++ ) {
406 *xco = cx[j] * (*xci) ;
407 xci++ ;
408 xco++ ;
409 } // Fin de la boucle sur r
410 } // Fin de la boucle sur theta
411
412 // k = 1
413 xci += nr*nt ;
414 xco += nr*nt ;
415
416 // k >=2
417 for (int k=2 ; k<np+1 ; k++) {
418 for (int j=0 ; j<nt ; j++) {
419 for (int i=0 ; i<nr ; i++ ) {
420 *xco = cx[j] * (*xci) ;
421 xci++ ;
422 xco++ ;
423 } // Fin de la boucle sur r
424 } // Fin de la boucle sur theta
425 } // Fin de la boucle sur phi
426
427 // On remet les choses la ou il faut
428 delete [] tb->t ;
429 tb->t = xo ;
430
431 // base de developpement
432 int base_r = b & MSQ_R ;
433 int base_p = b & MSQ_P ;
434 b = base_r | base_p | T_COS_P ;
435}
436
437// cas T_SIN_I
438//------------
439void _dsdtet_t_sin_i(Tbl* tb, int & b)
440{
441
442 // Peut-etre rien a faire ?
443 if (tb->get_etat() == ETATZERO) {
444 int base_r = b & MSQ_R ;
445 int base_p = b & MSQ_P ;
446 b = base_r | base_p | T_COS_I ;
447 return ;
448 }
449
450 // Protection
451 assert(tb->get_etat() == ETATQCQ) ;
452
453 // Pour le confort
454 int nr = (tb->dim).dim[0] ; // Nombre
455 int nt = (tb->dim).dim[1] ; // de points
456 int np = (tb->dim).dim[2] ; // physiques REELS
457 np = np - 2 ; // Nombre de points physiques
458
459 // Variables statiques
460 static double* cx = 0 ;
461 static int nt_pre =0 ;
462
463 // Test sur np pour initialisation eventuelle
464 {
465 if (nt > nt_pre) {
466 nt_pre = nt ;
467 cx = reinterpret_cast<double*>(realloc(cx, nt * sizeof(double))) ;
468 for (int i=0 ; i<nt ; i++) {
469 cx[i] = (2*i+1) ;
470 }
471 }
472 } // Fin de region critique
473
474 // pt. sur le tableau de double resultat
475 double* xo = new double[(tb->dim).taille] ;
476
477 // Initialisation a zero :
478 for (int i=0; i<(tb->dim).taille; i++) {
479 xo[i] = 0 ;
480 }
481
482 // On y va...
483 double* xi = tb->t ;
484 double* xci = xi ; // Pointeurs
485 double* xco = xo ; // courants
486
487 // k = 0
488 for (int j=0 ; j<nt ; j++) {
489 for (int i=0 ; i<nr ; i++ ) {
490 *xco = cx[j] * (*xci) ;
491 xci++ ;
492 xco++ ;
493 } // Fin de la boucle sur r
494 } // Fin de la boucle sur theta
495
496 // k = 1
497 xci += nr*nt ;
498 xco += nr*nt ;
499
500 // k >=2
501 for (int k=2 ; k<np+1 ; k++) {
502 for (int j=0 ; j<nt ; j++) {
503 for (int i=0 ; i<nr ; i++ ) {
504 *xco = cx[j] * (*xci) ;
505 xci++ ;
506 xco++ ;
507 } // Fin de la boucle sur r
508 } // Fin de la boucle sur theta
509 } // Fin de la boucle sur phi
510
511 // On remet les choses la ou il faut
512 delete [] tb->t ;
513 tb->t = xo ;
514
515 // base de developpement
516 int base_r = b & MSQ_R ;
517 int base_p = b & MSQ_P ;
518 b = base_r | base_p | T_COS_I ;
519}
520
521// cas T_COS_I
522//------------
523void _dsdtet_t_cos_i(Tbl* tb, int & b)
524{
525
526 // Peut-etre rien a faire ?
527 if (tb->get_etat() == ETATZERO) {
528 int base_r = b & MSQ_R ;
529 int base_p = b & MSQ_P ;
530 b = base_r | base_p | T_SIN_I ;
531 return ;
532 }
533
534 // Protection
535 assert(tb->get_etat() == ETATQCQ) ;
536
537 // Pour le confort
538 int nr = (tb->dim).dim[0] ; // Nombre
539 int nt = (tb->dim).dim[1] ; // de points
540 int np = (tb->dim).dim[2] ; // physiques REELS
541 np = np - 2 ; // Nombre de points physiques
542
543 // Variables statiques
544 static double* cx = 0 ;
545 static int nt_pre =0 ;
546
547 // Test sur np pour initialisation eventuelle
548 {
549 if (nt > nt_pre) {
550 nt_pre = nt ;
551 cx = reinterpret_cast<double*>(realloc(cx, nt * sizeof(double))) ;
552 for (int i=0 ; i<nt ; i++) {
553 cx[i] = - (2*i+1) ;
554 }
555 }
556 } // Fin de region critique
557
558 // pt. sur le tableau de double resultat
559 double* xo = new double[(tb->dim).taille] ;
560
561 // Initialisation a zero :
562 for (int i=0; i<(tb->dim).taille; i++) {
563 xo[i] = 0 ;
564 }
565
566 // On y va...
567 double* xi = tb->t ;
568 double* xci = xi ; // Pointeurs
569 double* xco = xo ; // courants
570
571 // k = 0
572 for (int j=0 ; j<nt ; j++) {
573 for (int i=0 ; i<nr ; i++ ) {
574 *xco = cx[j] * (*xci) ;
575 xci++ ;
576 xco++ ;
577 } // Fin de la boucle sur r
578 } // Fin de la boucle sur theta
579
580 // k = 1
581 xci += nr*nt ;
582 xco += nr*nt ;
583
584 // k >=2
585 for (int k=2 ; k<np+1 ; k++) {
586 for (int j=0 ; j<nt ; j++) {
587 for (int i=0 ; i<nr ; i++ ) {
588 *xco = cx[j] * (*xci) ;
589 xci++ ;
590 xco++ ;
591 } // Fin de la boucle sur r
592 } // Fin de la boucle sur theta
593 } // Fin de la boucle sur phi
594
595 // On remet les choses la ou il faut
596 delete [] tb->t ;
597 tb->t = xo ;
598
599 // base de developpement
600 int base_r = b & MSQ_R ;
601 int base_p = b & MSQ_P ;
602 b = base_r | base_p | T_SIN_I ;
603}
604
605// cas T_COSSIN_CP
606//----------------
607void _dsdtet_t_cossin_cp(Tbl* tb, int & b)
608{
609
610 // Peut-etre rien a faire ?
611 if (tb->get_etat() == ETATZERO) {
612 int base_r = b & MSQ_R ;
613 int base_p = b & MSQ_P ;
614 b = base_r | base_p | T_COSSIN_SP ;
615 return ;
616 }
617
618 // Protection
619 assert(tb->get_etat() == ETATQCQ) ;
620
621 // Pour le confort
622 int nr = (tb->dim).dim[0] ; // Nombre
623 int nt = (tb->dim).dim[1] ; // de points
624 int np = (tb->dim).dim[2] ; // physiques REELS
625 np = np - 2 ; // Nombre de points physiques
626
627 // Variables statiques
628 static double* cxp = 0 ;
629 static double* cxi = 0 ;
630 static int nt_pre = 0 ;
631
632 // Test sur np pour initialisation eventuelle
633 {
634 if (nt > nt_pre) {
635 nt_pre = nt ;
636 cxp = reinterpret_cast<double*>(realloc(cxp, nt * sizeof(double))) ;
637 cxi = reinterpret_cast<double*>(realloc(cxi, nt * sizeof(double))) ;
638 for (int i=0 ; i<nt ; i++) {
639 cxp[i] = - (2*i) ;
640 cxi[i] = (2*i+1) ;
641 }
642 }
643 } // Fin de region critique
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 double* cx[2] ; // Tableau des Pointeur de coefficient
658
659 // Initialisation des pointeurs de coefficients
660 cx[0] = cxp ; // cos pairs pour m pair
661 cx[1] = cxi ; // sin impair pour m impair
662
663 // k = 0
664 // Choix de la parite
665 double* cxl = cx[0] ; // Pointeur de coefficients local
666 for (int j=0 ; j<nt ; j++) {
667 for (int i=0 ; i<nr ; i++) {
668 *xco = cxl[j] * (*xci) ;
669 xci++ ;
670 xco++ ;
671 } // Fin de la Boucle sur r
672 } // Fin de la boucle sur theta
673
674 // k = 1
675 xci += nr*nt ;
676 xco += nr*nt ;
677
678 // k >= 2
679 for (int k=2 ; k<np+1 ; k++) {
680 // Choix de la parite
681 int m = (k/2) % 2 ;
682 cxl = cx[m] ; // Pointeur de coefficients local
683 for (int j=0 ; j<nt ; j++) {
684 for (int i=0 ; i<nr ; i++) {
685 *xco = cxl[j] * (*xci) ;
686 xci++ ;
687 xco++ ;
688 } // Fin de la Boucle sur r
689 } // Fin de la boucle sur theta
690 } // Fin de la boucle sur phi
691
692 // On remet les choses la ou il faut
693 delete [] tb->t ;
694 tb->t = xo ;
695
696 // base de developpement
697 int base_r = b & MSQ_R ;
698 int base_p = b & MSQ_P ;
699 b = base_r | base_p | T_COSSIN_SP ;
700}
701
702// cas T_COSSIN_SP
703//----------------
704void _dsdtet_t_cossin_sp(Tbl* tb, int & b)
705{
706
707 // Peut-etre rien a faire ?
708 if (tb->get_etat() == ETATZERO) {
709 int base_r = b & MSQ_R ;
710 int base_p = b & MSQ_P ;
711 b = base_r | base_p | T_COSSIN_CP ;
712 return ;
713 }
714
715 // Protection
716 assert(tb->get_etat() == ETATQCQ) ;
717
718 // Pour le confort
719 int nr = (tb->dim).dim[0] ; // Nombre
720 int nt = (tb->dim).dim[1] ; // de points
721 int np = (tb->dim).dim[2] ; // physiques REELS
722 np = np - 2 ; // Nombre de points physiques
723
724 // Variables statiques
725 static double* cxp = 0 ;
726 static double* cxi = 0 ;
727 static int nt_pre =0 ;
728
729 // Test sur np pour initialisation eventuelle
730 {
731 if (nt > nt_pre) {
732 nt_pre = nt ;
733 cxp = reinterpret_cast<double*>(realloc(cxp, nt * sizeof(double))) ;
734 cxi = reinterpret_cast<double*>(realloc(cxi, nt * sizeof(double))) ;
735 for (int i=0 ; i<nt ; i++) {
736 cxp[i] = (2*i) ;
737 cxi[i] = - (2*i+1) ;
738 }
739 }
740 } // Fin de region critique
741
742 // pt. sur le tableau de double resultat
743 double* xo = new double[(tb->dim).taille] ;
744
745 // Initialisation a zero :
746 for (int i=0; i<(tb->dim).taille; i++) {
747 xo[i] = 0 ;
748 }
749
750 // On y va...
751 double* xi = tb->t ;
752 double* xci = xi ; // Pointeurs
753 double* xco = xo ; // courants
754 double* cx[2] ; // Tableau des Pointeur de coefficient
755
756 // Initialisation des pointeurs de coefficients
757 cx[0] = cxp ; // sin pairs pour m pair
758 cx[1] = cxi ; // cos impair pour m impair
759
760 // k = 0
761 // Choix de la parite
762 double* cxl = cx[0] ; // Pointeur de coefficients local
763 for (int j=0 ; j<nt ; j++) {
764 for (int i=0 ; i<nr ; i++) {
765 *xco = cxl[j] * (*xci) ;
766 xci++ ;
767 xco++ ;
768 } // Fin de la Boucle sur r
769 } // Fin de la boucle sur theta
770
771 // k = 1
772 xci += nr*nt ;
773 xco += nr*nt ;
774
775 // k >= 2
776 for (int k=2 ; k<np+1 ; k++) {
777 // Choix de la parite
778 int m = (k/2) % 2 ;
779 cxl = cx[m] ; // Pointeur de coefficients local
780 for (int j=0 ; j<nt ; j++) {
781 for (int i=0 ; i<nr ; i++) {
782 *xco = cxl[j] * (*xci) ;
783 xci++ ;
784 xco++ ;
785 } // Fin de la Boucle sur r
786 } // Fin de la boucle sur theta
787 } // Fin de la boucle sur phi
788
789 // On remet les choses la ou il faut
790 delete [] tb->t ;
791 tb->t = xo ;
792
793 // base de developpement
794 int base_r = b & MSQ_R ;
795 int base_p = b & MSQ_P ;
796 b = base_r | base_p | T_COSSIN_CP ;
797}
798
799// cas T_COSSIN_CI
800//----------------
801void _dsdtet_t_cossin_ci(Tbl* tb, int & b)
802{
803
804 // Peut-etre rien a faire ?
805 if (tb->get_etat() == ETATZERO) {
806 int base_r = b & MSQ_R ;
807 int base_p = b & MSQ_P ;
808 b = base_r | base_p | T_COSSIN_SI ;
809 return ;
810 }
811
812 // Protection
813 assert(tb->get_etat() == ETATQCQ) ;
814
815 // Pour le confort
816 int nr = (tb->dim).dim[0] ; // Nombre
817 int nt = (tb->dim).dim[1] ; // de points
818 int np = (tb->dim).dim[2] ; // physiques REELS
819 np = np - 2 ; // Nombre de points physiques
820
821 // Variables statiques
822 static double* cxp = 0 ;
823 static double* cxi = 0 ;
824 static int nt_pre =0 ;
825
826 // Test sur np pour initialisation eventuelle
827 {
828 if (nt > nt_pre) {
829 nt_pre = nt ;
830 cxp = reinterpret_cast<double*>(realloc(cxp, nt * sizeof(double))) ;
831 cxi = reinterpret_cast<double*>(realloc(cxi, nt * sizeof(double))) ;
832 for (int i=0 ; i<nt ; i++) {
833 cxp[i] = (2*i) ;
834 cxi[i] = - (2*i+1) ;
835 }
836 }
837 } // Fin de region critique
838
839 // pt. sur le tableau de double resultat
840 double* xo = new double[(tb->dim).taille] ;
841
842 // Initialisation a zero :
843 for (int i=0; i<(tb->dim).taille; i++) {
844 xo[i] = 0 ;
845 }
846
847 // On y va...
848 double* xi = tb->t ;
849 double* xci = xi ; // Pointeurs
850 double* xco = xo ; // courants
851 double* cx[2] ; // Tableau des Pointeur de coefficient
852
853 // Initialisation des pointeurs de coefficients
854 cx[0] = cxi ; // cos impairs pour m pair
855 cx[1] = cxp ; // sin pair pour m impair
856
857 // k = 0
858 // Choix de la parite
859 double* cxl = cx[0] ; // Pointeur de coefficients local
860 for (int j=0 ; j<nt ; j++) {
861 for (int i=0 ; i<nr ; i++) {
862 *xco = cxl[j] * (*xci) ;
863 xci++ ;
864 xco++ ;
865 } // Fin de la Boucle sur r
866 } // Fin de la boucle sur theta
867
868 // k = 1
869 xci += nr*nt ;
870 xco += nr*nt ;
871
872 // k >= 2
873 for (int k=2 ; k<np+1 ; k++) {
874 // Choix de la parite
875 int m = (k/2) % 2 ;
876 cxl = cx[m] ; // Pointeur de coefficients local
877 for (int j=0 ; j<nt ; j++) {
878 for (int i=0 ; i<nr ; i++) {
879 *xco = cxl[j] * (*xci) ;
880 xci++ ;
881 xco++ ;
882 } // Fin de la Boucle sur r
883 } // Fin de la boucle sur theta
884 } // Fin de la boucle sur phi
885
886 // On remet les choses la ou il faut
887 delete [] tb->t ;
888 tb->t = xo ;
889
890 // base de developpement
891 int base_r = b & MSQ_R ;
892 int base_p = b & MSQ_P ;
893 b = base_r | base_p | T_COSSIN_SI ;
894}
895
896// cas T_COSSIN_SI
897//----------------
898void _dsdtet_t_cossin_si(Tbl* tb, int & b)
899{
900
901 // Peut-etre rien a faire ?
902 if (tb->get_etat() == ETATZERO) {
903 int base_r = b & MSQ_R ;
904 int base_p = b & MSQ_P ;
905 b = base_r | base_p | T_COSSIN_CI ;
906 return ;
907 }
908
909 // Protection
910 assert(tb->get_etat() == ETATQCQ) ;
911
912 // Pour le confort
913 int nr = (tb->dim).dim[0] ; // Nombre
914 int nt = (tb->dim).dim[1] ; // de points
915 int np = (tb->dim).dim[2] ; // physiques REELS
916 np = np - 2 ; // Nombre de points physiques
917
918 // Variables statiques
919 static double* cxp = 0 ;
920 static double* cxi = 0 ;
921 static int nt_pre =0 ;
922
923 // Test sur np pour initialisation eventuelle
924 {
925 if (nt > nt_pre) {
926 nt_pre = nt ;
927 cxp = reinterpret_cast<double*>(realloc(cxp, nt * sizeof(double))) ;
928 cxi = reinterpret_cast<double*>(realloc(cxi, nt * sizeof(double))) ;
929 for (int i=0 ; i<nt ; i++) {
930 cxp[i] = - (2*i) ;
931 cxi[i] = (2*i+1) ;
932 }
933 }
934 } // Fin de region critique
935
936 // pt. sur le tableau de double resultat
937 double* xo = new double[(tb->dim).taille] ;
938
939 // Initialisation a zero :
940 for (int i=0; i<(tb->dim).taille; i++) {
941 xo[i] = 0 ;
942 }
943
944 // On y va...
945 double* xi = tb->t ;
946 double* xci = xi ; // Pointeurs
947 double* xco = xo ; // courants
948 double* cx[2] ; // Tableau des Pointeur de coefficient
949
950 // Initialisation des pointeurs de coefficients
951 cx[0] = cxi ; // sin impair pour m pair
952 cx[1] = cxp ; // cos pairs pour m impair
953
954 // k = 0
955 // Choix de la parite
956 double* cxl = cx[0] ; // Pointeur de coefficients local
957 for (int j=0 ; j<nt ; j++) {
958 for (int i=0 ; i<nr ; i++) {
959 *xco = cxl[j] * (*xci) ;
960 xci++ ;
961 xco++ ;
962 } // Fin de la Boucle sur r
963 } // Fin de la boucle sur theta
964
965 // k = 1
966 xci += nr*nt ;
967 xco += nr*nt ;
968
969 // k >= 2
970 for (int k=2 ; k<np+1 ; k++) {
971 // Choix de la parite
972 int m = (k/2) % 2 ;
973 cxl = cx[m] ; // Pointeur de coefficients local
974 for (int j=0 ; j<nt ; j++) {
975 for (int i=0 ; i<nr ; i++) {
976 *xco = cxl[j] * (*xci) ;
977 xci++ ;
978 xco++ ;
979 } // Fin de la Boucle sur r
980 } // Fin de la boucle sur theta
981 } // Fin de la boucle sur phi
982
983 // On remet les choses la ou il faut
984 delete [] tb->t ;
985 tb->t = xo ;
986
987 // base de developpement
988 int base_r = b & MSQ_R ;
989 int base_p = b & MSQ_P ;
990 b = base_r | base_p | T_COSSIN_CI ;
991}
992
993// cas T_COSSIN_C
994//----------------
995void _dsdtet_t_cossin_c(Tbl* tb, int & b)
996{
997
998 // Peut-etre rien a faire ?
999 if (tb->get_etat() == ETATZERO) {
1000 int base_r = b & MSQ_R ;
1001 int base_p = b & MSQ_P ;
1002 b = base_r | base_p | T_COSSIN_S ;
1003 return ;
1004 }
1005
1006 // Protection
1007 assert(tb->get_etat() == ETATQCQ) ;
1008
1009 // Pour le confort
1010 int nr = (tb->dim).dim[0] ; // Nombre
1011 int nt = (tb->dim).dim[1] ; // de points
1012 int np = (tb->dim).dim[2] ; // physiques REELS
1013 np = np - 2 ; // Nombre de points physiques
1014
1015 // Variables statiques
1016 static double* cxp = 0 ;
1017 static double* cxi = 0 ;
1018 static int nt_pre = 0 ;
1019
1020 // Test sur np pour initialisation eventuelle
1021 {
1022 if (nt > nt_pre) {
1023 nt_pre = nt ;
1024 cxp = reinterpret_cast<double*>(realloc(cxp, nt * sizeof(double))) ;
1025 cxi = reinterpret_cast<double*>(realloc(cxi, nt * sizeof(double))) ;
1026 for (int i=0 ; i<nt ; i++) {
1027 cxp[i] = - i ;
1028 cxi[i] = i ;
1029 }
1030 }
1031 } // Fin de region critique
1032
1033 // pt. sur le tableau de double resultat
1034 double* xo = new double[(tb->dim).taille] ;
1035
1036 // Initialisation a zero :
1037 for (int i=0; i<(tb->dim).taille; i++) {
1038 xo[i] = 0 ;
1039 }
1040
1041 // On y va...
1042 double* xi = tb->t ;
1043 double* xci = xi ; // Pointeurs
1044 double* xco = xo ; // courants
1045 double* cx[2] ; // Tableau des Pointeur de coefficient
1046
1047 // Initialisation des pointeurs de coefficients
1048 cx[0] = cxp ; // cos pairs pour m pair
1049 cx[1] = cxi ; // sin impair pour m impair
1050
1051 // k = 0
1052 // Choix de la parite
1053 double* cxl = cx[0] ; // Pointeur de coefficients local
1054 for (int j=0 ; j<nt ; j++) {
1055 for (int i=0 ; i<nr ; i++) {
1056 *xco = cxl[j] * (*xci) ;
1057 xci++ ;
1058 xco++ ;
1059 } // Fin de la Boucle sur r
1060 } // Fin de la boucle sur theta
1061
1062 // k = 1
1063 xci += nr*nt ;
1064 xco += nr*nt ;
1065
1066 // k >= 2
1067 for (int k=2 ; k<np+1 ; k++) {
1068 // Choix de la parite
1069 int m = (k/2) % 2 ;
1070 cxl = cx[m] ; // Pointeur de coefficients local
1071 for (int j=0 ; j<nt ; j++) {
1072 for (int i=0 ; i<nr ; i++) {
1073 *xco = cxl[j] * (*xci) ;
1074 xci++ ;
1075 xco++ ;
1076 } // Fin de la Boucle sur r
1077 } // Fin de la boucle sur theta
1078 } // Fin de la boucle sur phi
1079
1080 // On remet les choses la ou il faut
1081 delete [] tb->t ;
1082 tb->t = xo ;
1083
1084 // base de developpement
1085 int base_r = b & MSQ_R ;
1086 int base_p = b & MSQ_P ;
1087 b = base_r | base_p | T_COSSIN_S ;
1088}
1089
1090// cas T_COSSIN_S
1091//----------------
1092void _dsdtet_t_cossin_s(Tbl* tb, int & b)
1093{
1094
1095 // Peut-etre rien a faire ?
1096 if (tb->get_etat() == ETATZERO) {
1097 int base_r = b & MSQ_R ;
1098 int base_p = b & MSQ_P ;
1099 b = base_r | base_p | T_COSSIN_C ;
1100 return ;
1101 }
1102
1103 // Protection
1104 assert(tb->get_etat() == ETATQCQ) ;
1105
1106 // Pour le confort
1107 int nr = (tb->dim).dim[0] ; // Nombre
1108 int nt = (tb->dim).dim[1] ; // de points
1109 int np = (tb->dim).dim[2] ; // physiques REELS
1110 np = np - 2 ; // Nombre de points physiques
1111
1112 // Variables statiques
1113 static double* cxp = 0 ;
1114 static double* cxi = 0 ;
1115 static int nt_pre =0 ;
1116
1117 // Test sur np pour initialisation eventuelle
1118 {
1119 if (nt > nt_pre) {
1120 nt_pre = nt ;
1121 cxp = reinterpret_cast<double*>(realloc(cxp, nt * sizeof(double))) ;
1122 cxi = reinterpret_cast<double*>(realloc(cxi, nt * sizeof(double))) ;
1123 for (int i=0 ; i<nt ; i++) {
1124 cxp[i] = i ;
1125 cxi[i] = - i ;
1126 }
1127 }
1128 } // Fin de region critique
1129
1130 // pt. sur le tableau de double resultat
1131 double* xo = new double[(tb->dim).taille] ;
1132
1133 // Initialisation a zero :
1134 for (int i=0; i<(tb->dim).taille; i++) {
1135 xo[i] = 0 ;
1136 }
1137
1138 // On y va...
1139 double* xi = tb->t ;
1140 double* xci = xi ; // Pointeurs
1141 double* xco = xo ; // courants
1142 double* cx[2] ; // Tableau des Pointeur de coefficient
1143
1144 // Initialisation des pointeurs de coefficients
1145 cx[0] = cxp ; // sin pairs pour m pair
1146 cx[1] = cxi ; // cos impair pour m impair
1147
1148 // k = 0
1149 // Choix de la parite
1150 double* cxl = cx[0] ; // Pointeur de coefficients local
1151 for (int j=0 ; j<nt ; j++) {
1152 for (int i=0 ; i<nr ; i++) {
1153 *xco = cxl[j] * (*xci) ;
1154 xci++ ;
1155 xco++ ;
1156 } // Fin de la Boucle sur r
1157 } // Fin de la boucle sur theta
1158
1159 // k = 1
1160 xci += nr*nt ;
1161 xco += nr*nt ;
1162
1163 // k >= 2
1164 for (int k=2 ; k<np+1 ; k++) {
1165 // Choix de la parite
1166 int m = (k/2) % 2 ;
1167 cxl = cx[m] ; // Pointeur de coefficients local
1168 for (int j=0 ; j<nt ; j++) {
1169 for (int i=0 ; i<nr ; i++) {
1170 *xco = cxl[j] * (*xci) ;
1171 xci++ ;
1172 xco++ ;
1173 } // Fin de la Boucle sur r
1174 } // Fin de la boucle sur theta
1175 } // Fin de la boucle sur phi
1176
1177 // On remet les choses la ou il faut
1178 delete [] tb->t ;
1179 tb->t = xo ;
1180
1181 // base de developpement
1182 int base_r = b & MSQ_R ;
1183 int base_p = b & MSQ_P ;
1184 b = base_r | base_p | T_COSSIN_C ;
1185}
1186}
Lorene::Tbl
Basic array class.
Definition tbl.h:161
T_COSSIN_SP
#define T_COSSIN_SP
sin pair-cos impair alternes, sin pour m=0
Definition type_parite.h:210
T_SIN_P
#define T_SIN_P
dev. sin seulement, harmoniques paires
Definition type_parite.h:202
MSQ_R
#define MSQ_R
Extraction de l'info sur R.
Definition type_parite.h:152
T_COSSIN_S
#define T_COSSIN_S
dev. cos-sin alternes, sin pour m=0
Definition type_parite.h:194
T_COSSIN_SI
#define T_COSSIN_SI
sin impair-cos pair alternes, sin pour m=0
Definition type_parite.h:214
T_COS_P
#define T_COS_P
dev. cos seulement, harmoniques paires
Definition type_parite.h:200
T_COSSIN_CI
#define T_COSSIN_CI
cos impair-sin pair alternes, cos pour m=0
Definition type_parite.h:212
T_COSSIN_CP
#define T_COSSIN_CP
cos pair-sin impair alternes, cos pour m=0
Definition type_parite.h:208
T_SIN_I
#define T_SIN_I
dev. sin seulement, harmoniques impaires
Definition type_parite.h:206
T_COS
#define T_COS
dev. cos seulement
Definition type_parite.h:196
MSQ_P
#define MSQ_P
Extraction de l'info sur Phi.
Definition type_parite.h:156
T_SIN
#define T_SIN
dev. sin seulement
Definition type_parite.h:198
T_COS_I
#define T_COS_I
dev. cos seulement, harmoniques impaires
Definition type_parite.h:204
T_COSSIN_C
#define T_COSSIN_C
dev. cos-sin alternes, cos pour m=0
Definition type_parite.h:192
Lorene
Lorene prototypes.
Definition app_hor.h:67