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