LORENE
op_ssint.C
1/*
2 * Copyright (c) 1999-2000 Jean-Alain Marck
3 * Copyright (c) 1999-2001 Eric Gourgoulhon
4 * Copyright (c) 2004 Michael Forot
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 division par rapport a sint theta
29 * (Utilisation interne)
30 *
31 * void _ssint_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_ssint.C,v 1.9 2016/12/05 16:18:08 j_novak Exp $
39 * $Log: op_ssint.C,v $
40 * Revision 1.9 2016/12/05 16:18:08 j_novak
41 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
42 *
43 * Revision 1.8 2014/10/13 08:53:26 j_novak
44 * Lorene classes and functions now belong to the namespace Lorene.
45 *
46 * Revision 1.7 2009/10/10 18:28:11 j_novak
47 * New bases T_COS and T_SIN.
48 *
49 * Revision 1.6 2007/12/21 10:43:38 j_novak
50 * Corrected some bugs in the case nt=1 (1 point in theta).
51 *
52 * Revision 1.5 2007/10/05 12:37:20 j_novak
53 * Corrected a few errors in the theta-nonsymmetric case (bases T_COSSIN_C and
54 * T_COSSIN_S).
55 *
56 * Revision 1.4 2005/02/16 15:29:48 m_forot
57 * Correct T_COSSIN_S and T_COSSIN_C cases
58 *
59 * Revision 1.3 2004/12/17 13:20:18 m_forot
60 * Modify the case T_COSSIN_C and T_COSSIN_S
61 *
62 * Revision 1.2 2004/11/23 15:16:01 m_forot
63 *
64 * Added the bases for the cases without any equatorial symmetry
65 * (T_COSSIN_C, T_COSSIN_S, T_LEG, R_CHEBPI_P, R_CHEBPI_I).
66 *
67 * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
68 * LORENE
69 *
70 * Revision 2.4 2000/02/24 16:42:49 eric
71 * Initialisation a zero du tableau xo avant le calcul.
72 *
73 *
74 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/op_ssint.C,v 1.9 2016/12/05 16:18:08 j_novak Exp $
75 *
76 */
77
78// Fichier includes
79#include "tbl.h"
80
81
82 //-----------------------------------
83 // Routine pour les cas non prevus --
84 //-----------------------------------
85
86namespace Lorene {
87void _ssint_pas_prevu(Tbl * tb, int& base) {
88 cout << "ssint pas prevu..." << endl ;
89 cout << "Tbl: " << tb << " base: " << base << endl ;
90 abort () ;
91 exit(-1) ;
92}
93
94 //--------------
95 // cas T_COS
96 //--------------
97
98void _ssint_t_cos(Tbl* tb, int & b)
99{
100
101 // Peut-etre rien a faire ?
102 if (tb->get_etat() == ETATZERO) {
103 int base_r = b & MSQ_R ;
104 int base_p = b & MSQ_P ;
105 switch(base_r){
106 case(R_CHEBPI_P):
107 b = R_CHEBPI_I | base_p | T_SIN ;
108 break ;
109 case(R_CHEBPI_I):
110 b = R_CHEBPI_P | base_p | T_SIN ;
111 break ;
112 default:
113 b = base_r | base_p | T_SIN ;
114 break;
115 }
116 return ;
117 }
118
119 // Protection
120 assert(tb->get_etat() == ETATQCQ) ;
121
122 // Pour le confort : nbre de points reels.
123 int nr = (tb->dim).dim[0] ;
124 int nt = (tb->dim).dim[1] ;
125 int np = (tb->dim).dim[2] ;
126 np = np - 2 ;
127
128 // zone de travail
129 double* somP = new double [nr] ;
130 double* somN = new double [nr] ;
131
132 // pt. sur le tableau de double resultat
133 double* xo = new double[(tb->dim).taille] ;
134
135 // Initialisation a zero :
136 for (int i=0; i<(tb->dim).taille; i++) {
137 xo[i] = 0 ;
138 }
139
140 // On y va...
141 double* xi = tb->t ;
142 double* xci = xi ; // Pointeurs
143 double* xco = xo ; // courants
144
145 double cx ;
146
147 // k = 0
148
149 // Dernier j: j = nt-1
150 // Positionnement
151 xci += nr * (nt-1) ;
152 cx = -2 ;
153 xco += nr * (nt-1) ;
154
155 // Initialisation de som
156 for (int i=0 ; i<nr ; i++) {
157 somP[i] = 0. ;
158 somN[i] = 0. ;
159 xco[i] = 0. ; // mise a zero du dernier coef en theta
160 }
161
162 // j suivants
163 for (int j=nt-2 ; j >0 ; j--) {
164 int l = j % 2 ;
165 // Positionnement
166 xci -= nr ;
167 xco -= nr ;
168 // Calcul
169 for (int i=0 ; i<nr ; i++ ) {
170 if(l==1) somN[i] += cx * xci[i] ;
171 else somP[i] += cx * xci[i] ;
172 xco[i] = somN[i]*(1-l)+somP[i]*l ;
173 } // Fin de la boucle sur r
174 } // Fin de la boucle sur theta
175 // j = 0 sin(0*theta)
176 xci -= nr ;
177 xco -= nr ;
178 for (int i=0 ; i<nr ; i++) {
179 xco[i] = 0 ;
180 }
181 // Positionnement phi suivant
182
183 xci += nr*nt ;
184 xco += nr*nt ;
185
186 // k = 1
187 xci += nr*nt ;
188 xco += nr*nt ;
189
190 // k >= 2
191 for (int k=2 ; k<np+1 ; k++) {
192
193 // Dernier j: j = nt-1
194 // Positionnement
195 cx = -2 ;
196 xci += nr * (nt-1) ;
197 xco += nr * (nt-1) ;
198
199 // Initialisation de som
200 for (int i=0 ; i<nr ; i++) {
201 somP[i] = 0. ;
202 somN[i] = 0. ;
203 xco[i] = 0. ; // mise a zero dui dernier coefficient en theta.
204 }
205
206 // j suivants
207 for (int j=nt-2 ; j >= 0 ; j--) {
208 int l = j% 2 ;
209 // Positionnement
210 xci -= nr ;
211 xco -= nr ;
212 // Calcul
213 for (int i=0 ; i<nr ; i++ ) {
214 if(l==1) somN[i] += -2 * xci[i] ;
215 else somP[i] += -2 * xci[i] ;
216 xco[i] = somN[i]*(1-l)+somP[i]*l ;
217 } // Fin de la boucle sur r
218 } // Fin de la boucle sur theta
219 for (int i=0 ; i<nr ; i++) {
220 xco[i] = 0. ;
221 }
222
223 // Positionnement phi suivant
224 xci += nr*nt ;
225 xco += nr*nt ;
226
227 } // Fin de la boucle sur phi
228
229 // On remet les choses la ou il faut
230 delete [] tb->t ;
231 tb->t = xo ;
232
233 //Menage
234 delete [] somP ;
235 delete [] somN ;
236
237 // base de developpement
238 int base_r = b & MSQ_R ;
239 int base_p = b & MSQ_P ;
240 switch(base_r){
241 case(R_CHEBPI_P):
242 b = R_CHEBPI_I | base_p | T_SIN ;
243 break ;
244 case(R_CHEBPI_I):
245 b = R_CHEBPI_P | base_p | T_SIN ;
246 break ;
247 default:
248 b = base_r | base_p | T_SIN ;
249 break;
250 }
251}
252
253 //------------
254 // cas T_SIN
255 //------------
256
257void _ssint_t_sin(Tbl* tb, int & b)
258{
259
260
261 // Peut-etre rien a faire ?
262 if (tb->get_etat() == ETATZERO) {
263 int base_r = b & MSQ_R ;
264 int base_p = b & MSQ_P ;
265 switch(base_r){
266 case(R_CHEBPI_P):
267 b = R_CHEBPI_I | base_p | T_COS ;
268 break ;
269 case(R_CHEBPI_I):
270 b = R_CHEBPI_P | base_p | T_COS ;
271 break ;
272 default:
273 b = base_r | base_p | T_COS ;
274 break;
275 }
276 return ;
277 }
278
279 // Protection
280 assert(tb->get_etat() == ETATQCQ) ;
281
282 // Pour le confort : nbre de points reels.
283 int nr = (tb->dim).dim[0] ;
284 int nt = (tb->dim).dim[1] ;
285 int np = (tb->dim).dim[2] ;
286 np = np - 2 ;
287
288 // zone de travail
289 double* somP = new double [nr] ;
290 double* somN = new double [nr] ;
291
292 // pt. sur le tableau de double resultat
293 double* xo = new double[(tb->dim).taille] ;
294
295 // Initialisation a zero :
296 for (int i=0; i<(tb->dim).taille; i++) {
297 xo[i] = 0 ;
298 }
299
300 // On y va...
301 double* xi = tb->t ;
302 double* xci = xi ; // Pointeurs
303 double* xco = xo ; // courants
304
305 // k = 0
306
307 // Dernier j: j = nt-1
308 // Positionnement
309 xci += nr * (nt-1) ;
310 xco += nr * (nt-1) ;
311
312 // Initialisation de som
313 for (int i=0 ; i<nr ; i++) {
314 somP[i] = 0. ;
315 somN[i] = 0. ;
316 xco[i] = 0. ;
317 }
318
319 // j suivants
320 for (int j=nt-2 ; j >= 0 ; j--) {
321 int l = j % 2 ;
322 // Positionnement
323 xci -= nr ;
324 xco -= nr ;
325 // Calcul
326 for (int i=0 ; i<nr ; i++ ) {
327 if(l==1) somN[i] += 2 * xci[i] ;
328 else somP[i] += 2 * xci[i] ;
329 xco[i] = somN[i]*(1-l)+somP[i]*l ;
330 } // Fin de la boucle sur r
331 } // Fin de la boucle sur theta
332 for (int i=0 ; i<nr ; i++) {
333 xco[i] *= .5 ;
334 }
335 // Positionnement phi suivant
336 xci += nr*nt ;
337 xco += nr*nt ;
338
339 // k = 1
340 xci += nr*nt ;
341 xco += nr*nt ;
342
343 // k >= 2
344 for (int k=2 ; k<np+1 ; k++) {
345
346 // Dernier j: j = nt-1
347 // Positionnement
348 xci += nr * (nt-1) ;
349 xco += nr * (nt-1) ;
350
351 // Initialisation de som
352 for (int i=0 ; i<nr ; i++) {
353 somP[i] = 0. ;
354 somN[i] = 0. ;
355 xco[i] = 0. ;
356 }
357
358 // j suivants
359 for (int j=nt-2 ; j >= 0 ; j--) {
360 int l = j % 2 ;
361 // Positionnement
362 xci -= nr ;
363 xco -= nr ;
364 // Calcul
365 for (int i=0 ; i<nr ; i++ ) {
366 if(l==1) somN[i] += 2 * xci[i] ;
367 else somP[i] += 2 * xci[i] ;
368 xco[i] = somN[i]*(1-l)+somP[i]*l ;
369 } // Fin de la boucle sur r
370 } // Fin de la boucle sur theta
371
372 // Normalisation du premier theta
373 for (int i=0 ; i<nr ; i++) {
374 xco[i] *= .5 ;
375 }
376
377 // Positionnement phi suivant
378
379 xci += nr*nt ;
380 xco += nr*nt ;
381
382 } // Fin de la boucle sur phi
383
384 // On remet les choses la ou il faut
385 delete [] tb->t ;
386 tb->t = xo ;
387
388 //Menage
389 delete [] somP ;
390 delete [] somN ;
391
392 // base de developpement
393 int base_r = b & MSQ_R ;
394 int base_p = b & MSQ_P ;
395 switch(base_r){
396 case(R_CHEBPI_P):
397 b = R_CHEBPI_I | base_p | T_COS ;
398 break ;
399 case(R_CHEBPI_I):
400 b = R_CHEBPI_P | base_p | T_COS ;
401 break ;
402 default:
403 b = base_r | base_p | T_COS ;
404 break;
405 }
406}
407 //----------------
408 // cas T_COS_P
409 //----------------
410
411void _ssint_t_cos_p(Tbl* tb, int & b)
412{
413
414 // Peut-etre rien a faire ?
415 if (tb->get_etat() == ETATZERO) {
416 int base_r = b & MSQ_R ;
417 int base_p = b & MSQ_P ;
418 b = base_r | base_p | T_SIN_I ;
419 return ;
420 }
421
422 // Protection
423 assert(tb->get_etat() == ETATQCQ) ;
424
425 // Pour le confort : nbre de points reels.
426 int nr = (tb->dim).dim[0] ;
427 int nt = (tb->dim).dim[1] ;
428 int np = (tb->dim).dim[2] ;
429 np = np - 2 ;
430
431 // zone de travail
432 double* som = new double [nr] ;
433
434 // pt. sur le tableau de double resultat
435 double* xo = new double[(tb->dim).taille] ;
436
437 // Initialisation a zero :
438 for (int i=0; i<(tb->dim).taille; i++) {
439 xo[i] = 0 ;
440 }
441
442 // On y va...
443 double* xi = tb->t ;
444 double* xci = xi ; // Pointeurs
445 double* xco = xo ; // courants
446
447 // k = 0
448
449 // Dernier j: j = nt-1
450 // Positionnement
451 xci += nr * (nt) ;
452 xco += nr * (nt-1) ;
453
454 // Initialisation de som
455 for (int i=0 ; i<nr ; i++) {
456 som[i] = 0. ;
457 xco[i] = 0. ; // mise a zero du dernier coef en theta
458 }
459
460 // j suivants
461 for (int j=nt-2 ; j >= 0 ; j--) {
462 // Positionnement
463 xci -= nr ;
464 xco -= nr ;
465 // Calcul
466 for (int i=0 ; i<nr ; i++ ) {
467 som[i] -= 2*xci[i] ;
468 xco[i] = som[i] ;
469 } // Fin de la boucle sur r
470 } // Fin de la boucle sur theta
471 // Positionnement phi suivant
472 xci -= nr ;
473 xci += nr*nt ;
474 xco += nr*nt ;
475
476 // k = 1
477 xci += nr*nt ;
478 xco += nr*nt ;
479
480 // k >= 2
481 for (int k=2 ; k<np+1 ; k++) {
482
483 // Dernier j: j = nt-1
484 // Positionnement
485 xci += nr * (nt) ;
486 xco += nr * (nt-1) ;
487
488 // Initialisation de som
489 for (int i=0 ; i<nr ; i++) {
490 som[i] = 0. ;
491 xco[i] = 0. ; // mise a zero dui dernier coefficient en theta.
492 }
493
494 // j suivants
495 for (int j=nt-2 ; j >= 0 ; j--) {
496 // Positionnement
497 xci -= nr ;
498 xco -= nr ;
499 // Calcul
500 for (int i=0 ; i<nr ; i++ ) {
501 som[i] -= 2*xci[i] ;
502 xco[i] = som[i] ;
503 } // Fin de la boucle sur r
504 } // Fin de la boucle sur theta
505 // Positionnement phi suivant
506 xci -= nr ;
507 xci += nr*nt ;
508 xco += nr*nt ;
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 //Menage
516 delete [] som ;
517
518 // base de developpement
519 int base_r = b & MSQ_R ;
520 int base_p = b & MSQ_P ;
521 b = base_r | base_p | T_SIN_I ;
522}
523
524 //--------------
525 // cas T_SIN_P
526 //--------------
527
528void _ssint_t_sin_p(Tbl* tb, int & b)
529{
530
531
532 // Peut-etre rien a faire ?
533 if (tb->get_etat() == ETATZERO) {
534 int base_r = b & MSQ_R ;
535 int base_p = b & MSQ_P ;
536 b = base_r | base_p | T_COS_I ;
537 return ;
538 }
539
540 // Protection
541 assert(tb->get_etat() == ETATQCQ) ;
542
543 // Pour le confort : nbre de points reels.
544 int nr = (tb->dim).dim[0] ;
545 int nt = (tb->dim).dim[1] ;
546 int np = (tb->dim).dim[2] ;
547 np = np - 2 ;
548
549 // zone de travail
550 double* som = new double [nr] ;
551
552 // pt. sur le tableau de double resultat
553 double* xo = new double[(tb->dim).taille] ;
554
555 // Initialisation a zero :
556 for (int i=0; i<(tb->dim).taille; i++) {
557 xo[i] = 0 ;
558 }
559
560 // On y va...
561 double* xi = tb->t ;
562 double* xci = xi ; // Pointeurs
563 double* xco = xo ; // courants
564
565 // k = 0
566
567 // Dernier j: j = nt-1
568 // Positionnement
569 xci += nr * (nt-1) ;
570 xco += nr * (nt-1) ;
571
572 // Initialisation de som
573 for (int i=0 ; i<nr ; i++) {
574 som[i] = 0. ;
575 xco[i] = 0. ;
576 }
577 if (nt > 1) {
578 xco -= nr ;
579 for (int i=0 ; i<nr ; i++) xco[i] = 0 ;
580 }
581
582 // j suivants
583 for (int j=nt-2 ; j >= 1 ; j--) {
584 // Positionnement
585 xci -= nr ;
586 xco -= nr ;
587 // Calcul
588 for (int i=0 ; i<nr ; i++ ) {
589 som[i] += 2*xci[i] ;
590 xco[i] = som[i] ;
591 } // Fin de la boucle sur r
592 } // Fin de la boucle sur theta
593 // Positionnement phi suivant
594 xci -= nr ;
595 xci += nr*nt ;
596 xco += nr*nt ;
597
598 // k = 1
599 xci += nr*nt ;
600 xco += nr*nt ;
601
602 // k >= 2
603 for (int k=2 ; k<np+1 ; k++) {
604
605 // Dernier j: j = nt-1
606 // Positionnement
607 xci += nr * (nt-1) ;
608 xco += nr * (nt-1) ;
609
610 // Initialisation de som
611 for (int i=0 ; i<nr ; i++) {
612 som[i] = 0. ;
613 xco[i] = 0. ;
614 }
615 xco -= nr ;
616 for (int i=0 ; i<nr ; i++) xco[i] = 0 ;
617
618 // j suivants
619 for (int j=nt-2 ; j >= 1 ; j--) {
620 // Positionnement
621 xci -= nr ;
622 xco -= nr ;
623 // Calcul
624 for (int i=0 ; i<nr ; i++ ) {
625 som[i] += 2*xci[i] ;
626 xco[i] = som[i] ;
627 } // Fin de la boucle sur r
628 } // Fin de la boucle sur theta
629 // Positionnement phi suivant
630 xci -= nr ;
631 xci += nr*nt ;
632 xco += nr*nt ;
633 } // Fin de la boucle sur phi
634
635 // On remet les choses la ou il faut
636 delete [] tb->t ;
637 tb->t = xo ;
638
639 //Menage
640 delete [] som ;
641
642 // base de developpement
643 int base_r = b & MSQ_R ;
644 int base_p = b & MSQ_P ;
645 b = base_r | base_p | T_COS_I ;
646
647}
648 //--------------
649 // cas T_SIN_I
650 //--------------
651
652void _ssint_t_sin_i(Tbl* tb, int & b)
653{
654
655
656 // Peut-etre rien a faire ?
657 if (tb->get_etat() == ETATZERO) {
658 int base_r = b & MSQ_R ;
659 int base_p = b & MSQ_P ;
660 b = base_r | base_p | T_COS_P ;
661 return ;
662 }
663
664 // Protection
665 assert(tb->get_etat() == ETATQCQ) ;
666
667 // Pour le confort : nbre de points reels.
668 int nr = (tb->dim).dim[0] ;
669 int nt = (tb->dim).dim[1] ;
670 int np = (tb->dim).dim[2] ;
671 np = np - 2 ;
672
673 // zone de travail
674 double* som = new double [nr] ;
675
676 // pt. sur le tableau de double resultat
677 double* xo = new double[(tb->dim).taille] ;
678
679 // Initialisation a zero :
680 for (int i=0; i<(tb->dim).taille; i++) {
681 xo[i] = 0 ;
682 }
683
684 // On y va...
685 double* xi = tb->t ;
686 double* xci = xi ; // Pointeurs
687 double* xco = xo ; // courants
688
689
690 // k = 0
691
692 // Dernier j: j = nt-1
693 // Positionnement
694 xci += nr * (nt-1) ;
695 xco += nr * (nt-1) ;
696
697 // Initialisation de som
698 for (int i=0 ; i<nr ; i++) {
699 xco[i] = 0. ;
700 som[i] = 0. ;
701 }
702
703 // j suivants
704 for (int j=nt-2 ; j >= 0 ; j--) {
705 // Positionnement
706 xci -= nr ;
707 xco -= nr ;
708 // Calcul
709 for (int i=0 ; i<nr ; i++ ) {
710 som[i] += 2*xci[i] ;
711 xco[i] = som[i] ;
712 } // Fin de la boucle sur r
713 } // Fin de la boucle sur theta
714 // Normalisation du premier theta
715 for (int i=0 ; i<nr ; i++) {
716 xco[i] *= .5 ;
717 }
718 // Positionnement phi suivant
719 xci += nr*nt ;
720 xco += nr*nt ;
721
722 // k = 1
723 xci += nr*nt ;
724 xco += nr*nt ;
725
726 // k >= 2
727 for (int k=2 ; k<np+1 ; k++) {
728
729 // Dernier j: j = nt-1
730 // Positionnement
731 xci += nr * (nt-1) ;
732 xco += nr * (nt-1) ;
733
734 // Initialisation de som
735 for (int i=0 ; i<nr ; i++) {
736 xco[i] = 0. ;
737 som[i] = 0. ;
738 }
739
740 // j suivants
741 for (int j=nt-2 ; j >= 0 ; j--) {
742 // Positionnement
743 xci -= nr ;
744 xco -= nr ;
745 // Calcul
746 for (int i=0 ; i<nr ; i++ ) {
747 som[i] += 2*xci[i] ;
748 xco[i] = som[i] ;
749 } // Fin de la boucle sur r
750 } // Fin de la boucle sur theta
751 // Normalisation du premier theta
752 for (int i=0 ; i<nr ; i++) {
753 xco[i] *= .5 ;
754 }
755 // Positionnement phi suivant
756 xci += nr*nt ;
757 xco += nr*nt ;
758 } // Fin de la boucle sur phi
759
760
761 // On remet les choses la ou il faut
762 delete [] tb->t ;
763 tb->t = xo ;
764
765 //Menage
766 delete [] som ;
767
768 // base de developpement
769 int base_r = b & MSQ_R ;
770 int base_p = b & MSQ_P ;
771 b = base_r | base_p | T_COS_P ;
772
773}
774 //---------------------
775 // cas T_COS_I
776 //---------------------
777void _ssint_t_cos_i(Tbl* tb, int & b)
778{
779
780 // Peut-etre rien a faire ?
781 if (tb->get_etat() == ETATZERO) {
782 int base_r = b & MSQ_R ;
783 int base_p = b & MSQ_P ;
784 b = base_r | base_p | T_SIN_P ;
785 return ;
786 }
787
788 // Protection
789 assert(tb->get_etat() == ETATQCQ) ;
790
791 // Pour le confort : nbre de points reels.
792 int nr = (tb->dim).dim[0] ;
793 int nt = (tb->dim).dim[1] ;
794 int np = (tb->dim).dim[2] ;
795 np = np - 2 ;
796
797 // zone de travail
798 double* som = new double [nr] ;
799
800 // pt. sur le tableau de double resultat
801 double* xo = new double[(tb->dim).taille] ;
802
803 // Initialisation a zero :
804 for (int i=0; i<(tb->dim).taille; i++) {
805 xo[i] = 0 ;
806 }
807
808 // On y va...
809 double* xi = tb->t ;
810 double* xci = xi ; // Pointeurs
811 double* xco = xo ; // courants
812
813 // k = 0
814
815 // Dernier j: j = nt-1
816 // Positionnement
817 xci += nr * (nt) ;
818 xco += nr * (nt-1) ;
819
820 // Initialisation de som
821 for (int i=0 ; i<nr ; i++) {
822 som[i] = 0. ;
823 xco[i] = 0. ; // mise a zero du dernier coef en theta
824 }
825
826 // j suivants
827 for (int j=nt-1 ; j >= 0 ; j--) {
828 // Positionnement
829 xci -= nr ;
830 // Calcul
831 for (int i=0 ; i<nr ; i++ ) {
832 som[i] -= 2*xci[i] ;
833 xco[i] = som[i] ;
834 } // Fin de la boucle sur r
835 xco -= nr ;
836 } // Fin de la boucle sur theta
837 // Positionnement phi suivant
838 xci += nr*nt ;
839 xco += nr ;
840 xco += nr*nt ;
841
842 // k = 1
843 xci += nr*nt ;
844 xco += nr*nt ;
845
846 // k >= 2
847 for (int k=2 ; k<np+1 ; k++) {
848
849 // Dernier j: j = nt-1
850 // Positionnement
851 xci += nr * (nt) ;
852 xco += nr * (nt-1) ;
853
854 // Initialisation de som
855 for (int i=0 ; i<nr ; i++) {
856 som[i] = 0. ;
857 xco[i] = 0. ; // mise a zero du dernier coef en theta
858 }
859
860 // j suivants
861 for (int j=nt-1 ; j >= 0 ; j--) {
862 // Positionnement
863 xci -= nr ;
864 // Calcul
865 for (int i=0 ; i<nr ; i++ ) {
866 som[i] -= 2*xci[i] ;
867 xco[i] = som[i] ;
868 } // Fin de la boucle sur r
869 xco -= nr ;
870 } // Fin de la boucle sur theta
871 // Positionnement phi suivant
872 xci += nr*nt ;
873 xco += nr ;
874 xco += nr*nt ;
875 } // Fin de la boucle sur phi
876
877 // On remet les choses la ou il faut
878 delete [] tb->t ;
879 tb->t = xo ;
880
881 //Menage
882 delete [] som ;
883
884 // base de developpement
885 int base_r = b & MSQ_R ;
886 int base_p = b & MSQ_P ;
887 b = base_r | base_p | T_SIN_P ;
888
889}
890 //----------------------
891 // cas T_COSSIN_CP
892 //----------------------
893void _ssint_t_cossin_cp(Tbl* tb, int & b)
894{
895
896 // Peut-etre rien a faire ?
897 if (tb->get_etat() == ETATZERO) {
898 int base_r = b & MSQ_R ;
899 int base_p = b & MSQ_P ;
900 b = base_r | base_p | T_COSSIN_SI ;
901 return ;
902 }
903
904 // Protection
905 assert(tb->get_etat() == ETATQCQ) ;
906
907 // Pour le confort : nbre de points reels.
908 int nr = (tb->dim).dim[0] ;
909 int nt = (tb->dim).dim[1] ;
910 int np = (tb->dim).dim[2] ;
911 np = np - 2 ;
912
913 // zone de travail
914 double* som = new double [nr] ;
915
916 // pt. sur le tableau de double resultat
917 double* xo = new double[(tb->dim).taille] ;
918
919 // Initialisation a zero :
920 for (int i=0; i<(tb->dim).taille; i++) {
921 xo[i] = 0 ;
922 }
923
924 // On y va...
925 double* xi = tb->t ;
926 double* xci = xi ; // Pointeurs
927 double* xco = xo ; // courants
928
929 double cx ;
930
931 // k = 0
932 int m = 0 ; // Parite de l'harmonique en phi
933
934 // Dernier j: j = nt-1
935 // Positionnement
936 xci += nr * (nt) ;
937 cx = -2 ;
938 xco += nr * (nt-1) ;
939
940 // Initialisation de som
941 for (int i=0 ; i<nr ; i++) {
942 som[i] = 0. ;
943 xco[i] = 0. ;
944 }
945
946 // j suivants
947 for (int j=nt-2 ; j >= 0 ; j--) {
948 // Positionnement
949 xci -= nr ;
950 xco -= nr ;
951 // Calcul
952 for (int i=0 ; i<nr ; i++ ) {
953 som[i] += cx * xci[i] ;
954 xco[i] = som[i] ;
955 } // Fin de la boucle sur r
956 } // Fin de la boucle sur theta
957 // Positionnement phi suivant
958 if (m == 0) {
959 xci -= nr ;
960 }
961 xci += nr*nt ;
962 xco += nr*nt ;
963
964 // k=1
965 xci += nr*nt ;
966 xco += nr*nt ;
967
968 // k >= 2
969 for (int k=2 ; k<np+1 ; k++) {
970 m = (k/2) % 2 ; // Parite de l'harmonique en phi
971
972 // Dernier j: j = nt-1
973 // Positionnement
974 if (m == 0) {
975 xci += nr * (nt) ;
976 cx = -2 ;
977 }
978 else {
979 xci += nr * (nt-1) ;
980 cx = 2 ;
981 }
982 xco += nr * (nt-1) ;
983
984 // Initialisation de som
985 for (int i=0 ; i<nr ; i++) {
986 som[i] = 0. ;
987 xco[i] = 0. ;
988 }
989
990 // j suivants
991 for (int j=nt-2 ; j >= 0 ; j--) {
992 // Positionnement
993 xci -= nr ;
994 xco -= nr ;
995 // Calcul
996 for (int i=0 ; i<nr ; i++ ) {
997 som[i] += cx * xci[i] ;
998 xco[i] = som[i] ;
999 } // Fin de la boucle sur r
1000 } // Fin de la boucle sur theta
1001 // Normalisation du premier theta dans le cas sin(impair)
1002 if (m == 1) {
1003 for (int i=0 ; i<nr ; i++) {
1004 xco[i] *= .5 ;
1005 }
1006 }
1007 // Positionnement phi suivant
1008 if (m == 0) {
1009 xci -= nr ;
1010 }
1011 xci += nr*nt ;
1012 xco += nr*nt ;
1013 } // Fin de la boucle sur phi
1014
1015 // On remet les choses la ou il faut
1016 delete [] tb->t ;
1017 tb->t = xo ;
1018
1019 //Menage
1020 delete [] som ;
1021
1022 // base de developpement
1023 int base_r = b & MSQ_R ;
1024 int base_p = b & MSQ_P ;
1025 b = base_r | base_p | T_COSSIN_SI ;
1026}
1027
1028 //------------------
1029 // cas T_COSSIN_CI
1030 //----------------
1031void _ssint_t_cossin_ci(Tbl* tb, int & b)
1032{
1033 // Peut-etre rien a faire ?
1034 if (tb->get_etat() == ETATZERO) {
1035 int base_r = b & MSQ_R ;
1036 int base_p = b & MSQ_P ;
1037 b = base_r | base_p | T_COSSIN_SP ;
1038 return ;
1039 }
1040
1041 // Protection
1042 assert(tb->get_etat() == ETATQCQ) ;
1043
1044 // Pour le confort : nbre de points reels.
1045 int nr = (tb->dim).dim[0] ;
1046 int nt = (tb->dim).dim[1] ;
1047 int np = (tb->dim).dim[2] ;
1048 np = np - 2 ;
1049
1050 // zone de travail
1051 double* som = new double [nr] ;
1052
1053 // pt. sur le tableau de double resultat
1054 double* xo = new double[(tb->dim).taille] ;
1055
1056 // Initialisation a zero :
1057 for (int i=0; i<(tb->dim).taille; i++) {
1058 xo[i] = 0 ;
1059 }
1060
1061 // On y va...
1062 double* xi = tb->t ;
1063 double* xci = xi ; // Pointeurs
1064 double* xco = xo ; // courants
1065
1066 // k = 0
1067 int m = 0 ; // Parite de l'harmonique en phi
1068
1069
1070 // Dernier j: j = nt-1
1071 // Positionnement
1072 xci += nr * (nt) ;
1073 xco += nr * (nt-1) ;
1074
1075 // Initialisation de som
1076 for (int i=0 ; i<nr ; i++) {
1077 som[i] = 0. ;
1078 xco[i] = 0. ; // mise a zero du dernier coef en theta
1079 }
1080
1081 // j suivants
1082 for (int j=nt-1 ; j >= 0 ; j--) {
1083 // Positionnement
1084 xci -= nr ;
1085 // Calcul
1086 for (int i=0 ; i<nr ; i++ ) {
1087 som[i] -= 2*xci[i] ;
1088 xco[i] = som[i] ;
1089 } // Fin de la boucle sur r
1090 xco -= nr ;
1091 } // Fin de la boucle sur theta
1092 // Positionnement phi suivant
1093 xci += nr*nt ;
1094 xco += nr ;
1095 xco += nr*nt ;
1096
1097 // k=1
1098 xci += nr*nt ;
1099 xco += nr*nt ;
1100
1101 // k >= 2
1102 for (int k=2 ; k<np+1 ; k++) {
1103 m = (k/2) % 2 ; // Parite de l'harmonique en phi
1104
1105 switch(m) {
1106 case 0: // m pair: cos(impair)
1107 // Dernier j: j = nt-1
1108 // Positionnement
1109 xci += nr * (nt) ;
1110 xco += nr * (nt-1) ;
1111
1112 // Initialisation de som
1113 for (int i=0 ; i<nr ; i++) {
1114 som[i] = 0. ;
1115 xco[i] = 0. ; // mise a zero du dernier coef en theta
1116 }
1117
1118 // j suivants
1119 for (int j=nt-1 ; j >= 0 ; j--) {
1120 // Positionnement
1121 xci -= nr ;
1122 // Calcul
1123 for (int i=0 ; i<nr ; i++ ) {
1124 som[i] -= 2*xci[i] ;
1125 xco[i] = som[i] ;
1126 } // Fin de la boucle sur r
1127 xco -= nr ;
1128 } // Fin de la boucle sur theta
1129 // Positionnement phi suivant
1130 xci += nr*nt ;
1131 xco += nr ;
1132 xco += nr*nt ;
1133 break ;
1134
1135 case 1: // m impair: sin(impair)
1136 // Dernier j: j = nt-1
1137 // Positionnement
1138 xci += nr * (nt-1) ;
1139 xco += nr * (nt-1) ;
1140
1141 // Initialisation de som
1142 for (int i=0 ; i<nr ; i++) {
1143 som[i] = 0. ;
1144 xco[i] = 0. ;
1145 }
1146 xco -= nr ;
1147 for (int i=0 ; i<nr ; i++) xco[i] = 0 ;
1148
1149 // j suivants
1150 for (int j=nt-2 ; j >= 1 ; j--) {
1151 // Positionnement
1152 xci -= nr ;
1153 xco -= nr ;
1154 // Calcul
1155 for (int i=0 ; i<nr ; i++ ) {
1156 som[i] += 2*xci[i] ;
1157 xco[i] = som[i] ;
1158 } // Fin de la boucle sur r
1159 } // Fin de la boucle sur theta
1160 // Positionnement phi suivant
1161 xci -= nr ;
1162 xci += nr*nt ;
1163 xco += nr*nt ;
1164 break ;
1165 } // Fin du switch
1166 } // Fin de la boucle en phi
1167
1168 // On remet les choses la ou il faut
1169 delete [] tb->t ;
1170 tb->t = xo ;
1171
1172 //Menage
1173 delete [] som ;
1174
1175 // base de developpement
1176 int base_r = b & MSQ_R ;
1177 int base_p = b & MSQ_P ;
1178 b = base_r | base_p | T_COSSIN_SP ;
1179}
1180
1181 //---------------------
1182 // cas T_COSSIN_SI
1183 //----------------
1184void _ssint_t_cossin_si(Tbl* tb, int & b)
1185{
1186 // Peut-etre rien a faire ?
1187 if (tb->get_etat() == ETATZERO) {
1188 int base_r = b & MSQ_R ;
1189 int base_p = b & MSQ_P ;
1190 b = base_r | base_p | T_COSSIN_CP ;
1191 return ;
1192 }
1193
1194 // Protection
1195 assert(tb->get_etat() == ETATQCQ) ;
1196
1197 // Pour le confort : nbre de points reels.
1198 int nr = (tb->dim).dim[0] ;
1199 int nt = (tb->dim).dim[1] ;
1200 int np = (tb->dim).dim[2] ;
1201 np = np - 2 ;
1202
1203 // zone de travail
1204 double* som = new double [nr] ;
1205
1206 // pt. sur le tableau de double resultat
1207 double* xo = new double[(tb->dim).taille] ;
1208
1209 // Initialisation a zero :
1210 for (int i=0; i<(tb->dim).taille; i++) {
1211 xo[i] = 0 ;
1212 }
1213
1214 // On y va...
1215 double* xi = tb->t ;
1216 double* xci = xi ; // Pointeurs
1217 double* xco = xo ; // courants
1218
1219 double cx ;
1220
1221 // k = 0
1222 int m = 0 ; // Parite de l'harmonique en phi
1223
1224
1225 // Dernier j: j = nt-1
1226 // Positionnement
1227 xci += nr * (nt-1) ;
1228 cx = 2 ;
1229 xco += nr * (nt-1) ;
1230
1231
1232 // Initialisation de som
1233 for (int i=0 ; i<nr ; i++) {
1234 som[i] = 0. ;
1235 xco[i] = 0. ;
1236 }
1237
1238 // j suivants
1239 for (int j=nt-2 ; j >= 0 ; j--) {
1240 // Positionnement
1241 xci -= nr ;
1242 xco -= nr ;
1243 // Calcul
1244 for (int i=0 ; i<nr ; i++ ) {
1245 som[i] += cx * xci[i] ;
1246 xco[i] = som[i] ;
1247 } // Fin de la boucle sur r
1248 } // Fin de la boucle sur theta
1249 // Normalisation du premier theta dans le cas sin(impair)
1250 for (int i=0 ; i<nr ; i++) {
1251 xco[i] *= .5 ;
1252 }
1253
1254 xci += nr*nt ;
1255 xco += nr*nt ;
1256
1257 // k=1
1258 xci += nr*nt ;
1259 xco += nr*nt ;
1260
1261 // k >= 2
1262 for (int k=2 ; k<np+1 ; k++) {
1263 m = (k/2) % 2 ; // Parite de l'harmonique en phi
1264
1265 // Dernier j: j = nt-1
1266 // Positionnement
1267 if (m == 1) {
1268 xci += nr * (nt) ;
1269 cx = -2 ;
1270 }
1271 else {
1272 xci += nr * (nt-1) ;
1273 cx = 2 ;
1274 }
1275 xco += nr * (nt-1) ;
1276
1277 // Initialisation de som
1278 for (int i=0 ; i<nr ; i++) {
1279 som[i] = 0. ;
1280 xco[i] = 0. ;
1281 }
1282
1283 // j suivants
1284 for (int j=nt-2 ; j >= 0 ; j--) {
1285 // Positionnement
1286 xci -= nr ;
1287 xco -= nr ;
1288 // Calcul
1289 for (int i=0 ; i<nr ; i++ ) {
1290 som[i] += cx * xci[i] ;
1291 xco[i] = som[i] ;
1292 } // Fin de la boucle sur r
1293 } // Fin de la boucle sur theta
1294 // Normalisation du premier theta dans le cas sin(impair)
1295 if (m == 0) {
1296 for (int i=0 ; i<nr ; i++) {
1297 xco[i] *= .5 ;
1298 }
1299 }
1300 // Positionnement phi suivant
1301 if (m == 1) {
1302 xci -= nr ;
1303 }
1304 xci += nr*nt ;
1305 xco += nr*nt ;
1306 } // Fin de la boucle sur phi
1307
1308
1309 // On remet les choses la ou il faut
1310 delete [] tb->t ;
1311 tb->t = xo ;
1312
1313 //Menage
1314 delete [] som ;
1315
1316 // base de developpement
1317 int base_r = b & MSQ_R ;
1318 int base_p = b & MSQ_P ;
1319 b = base_r | base_p | T_COSSIN_CP ;
1320
1321
1322}
1323 //---------------------
1324 // cas T_COSSIN_SP
1325 //----------------
1326void _ssint_t_cossin_sp(Tbl* tb, int & b)
1327{
1328 // Peut-etre rien a faire ?
1329 if (tb->get_etat() == ETATZERO) {
1330 int base_r = b & MSQ_R ;
1331 int base_p = b & MSQ_P ;
1332 b = base_r | base_p | T_COSSIN_CI ;
1333 return ;
1334 }
1335
1336 // Protection
1337 assert(tb->get_etat() == ETATQCQ) ;
1338
1339 // Pour le confort : nbre de points reels.
1340 int nr = (tb->dim).dim[0] ;
1341 int nt = (tb->dim).dim[1] ;
1342 int np = (tb->dim).dim[2] ;
1343 np = np - 2 ;
1344
1345 // zone de travail
1346 double* som = new double [nr] ;
1347
1348 // pt. sur le tableau de double resultat
1349 double* xo = new double[(tb->dim).taille] ;
1350
1351 // Initialisation a zero :
1352 for (int i=0; i<(tb->dim).taille; i++) {
1353 xo[i] = 0 ;
1354 }
1355
1356 // On y va...
1357 double* xi = tb->t ;
1358 double* xci = xi ; // Pointeurs
1359 double* xco = xo ; // courants
1360
1361 double cx ;
1362
1363 // k = 0
1364 int m = 0 ; // Parite de l'harmonique en phi
1365
1366
1367 // Dernier j: j = nt-1
1368 // Positionnement
1369 xci += nr * (nt-1) ;
1370 cx = 2 ;
1371 xco += nr * (nt-1) ;
1372
1373
1374 // Initialisation de som
1375 for (int i=0 ; i<nr ; i++) {
1376 som[i] = 0. ;
1377 xco[i] = 0. ;
1378 }
1379 xco -= nr ;
1380 for (int i=0 ; i<nr ; i++) xco[i] = 0 ;
1381
1382 // j suivants
1383 for (int j=nt-2 ; j >= 1 ; j--) {
1384 // Positionnement
1385 xci -= nr ;
1386 xco -= nr ;
1387 // Calcul
1388 for (int i=0 ; i<nr ; i++ ) {
1389 som[i] += cx * xci[i] ;
1390 xco[i] = som[i] ;
1391 } // Fin de la boucle sur r
1392 } // Fin de la boucle sur theta
1393 xci += nr*nt ;
1394 xci -= nr ;
1395 xco += nr*nt ;
1396
1397 // k=1
1398 xci += nr*nt ;
1399 xco += nr*nt ;
1400
1401 for (int k=2 ; k<np+1 ; k++) {
1402 m = (k/2) % 2 ; // Parite de l'harmonique en phi
1403
1404 switch(m) {
1405 case 1: // m impair: cos(impair)
1406 // Dernier j: j = nt-1
1407 // Positionnement
1408 xci += nr * (nt) ;
1409 xco += nr * (nt-1) ;
1410
1411 // Initialisation de som
1412 for (int i=0 ; i<nr ; i++) {
1413 som[i] = 0. ;
1414 xco[i] = 0. ; // mise a zero du dernier coef en theta
1415 }
1416
1417 // j suivants
1418 for (int j=nt-1 ; j >= 0 ; j--) {
1419 // Positionnement
1420 xci -= nr ;
1421 // Calcul
1422 for (int i=0 ; i<nr ; i++ ) {
1423 som[i] -= 2*xci[i] ;
1424 xco[i] = som[i] ;
1425 } // Fin de la boucle sur r
1426 xco -= nr ;
1427 } // Fin de la boucle sur theta
1428 // Positionnement phi suivant
1429 xci += nr*nt ;
1430 xco += nr ;
1431 xco += nr*nt ;
1432 break ;
1433
1434 case 0: // m pair: sin(pair)
1435 // Dernier j: j = nt-1
1436 // Positionnement
1437 xci += nr * (nt-1) ;
1438 xco += nr * (nt-1) ;
1439
1440 // Initialisation de som
1441 for (int i=0 ; i<nr ; i++) {
1442 som[i] = 0. ;
1443 xco[i] = 0. ;
1444 }
1445 xco -= nr ;
1446 for (int i=0 ; i<nr ; i++) xco[i] = 0 ;
1447
1448 // j suivants
1449 for (int j=nt-2 ; j >= 1 ; j--) {
1450 // Positionnement
1451 xci -= nr ;
1452 xco -= nr ;
1453 // Calcul
1454 for (int i=0 ; i<nr ; i++ ) {
1455 som[i] += 2*xci[i] ;
1456 xco[i] = som[i] ;
1457 } // Fin de la boucle sur r
1458 } // Fin de la boucle sur theta
1459 // Positionnement phi suivant
1460 xci -= nr ;
1461 xci += nr*nt ;
1462 xco += nr*nt ;
1463 break ;
1464 } // Fin du switch
1465 } // Fin de la boucle en phi
1466
1467 // On remet les choses la ou il faut
1468 delete [] tb->t ;
1469 tb->t = xo ;
1470
1471 //Menage
1472 delete [] som ;
1473
1474 // base de developpement
1475 int base_r = b & MSQ_R ;
1476 int base_p = b & MSQ_P ;
1477 b = base_r | base_p | T_COSSIN_CI ;
1478
1479
1480}
1481
1482 //----------------------
1483 // cas T_COSSIN_C
1484 //----------------------
1485void _ssint_t_cossin_c(Tbl* tb, int & b)
1486{
1487
1488
1489 // Peut-etre rien a faire ?
1490 if (tb->get_etat() == ETATZERO) {
1491 int base_r = b & MSQ_R ;
1492 int base_p = b & MSQ_P ;
1493 switch(base_r){
1494 case(R_CHEBPI_P):
1495 b = R_CHEBPI_I | base_p | T_COSSIN_S ;
1496 break ;
1497 case(R_CHEBPI_I):
1498 b = R_CHEBPI_P | base_p | T_COSSIN_S ;
1499 break ;
1500 default:
1501 b = base_r | base_p | T_COSSIN_S ;
1502 break;
1503 }
1504 return ;
1505 }
1506
1507 // Protection
1508 assert(tb->get_etat() == ETATQCQ) ;
1509
1510 // Pour le confort : nbre de points reels.
1511 int nr = (tb->dim).dim[0] ;
1512 int nt = (tb->dim).dim[1] ;
1513 int np = (tb->dim).dim[2] ;
1514 np = np - 2 ;
1515
1516 // zone de travail
1517 double* somP = new double [nr] ;
1518 double* somN = new double [nr] ;
1519
1520 // pt. sur le tableau de double resultat
1521 double* xo = new double[(tb->dim).taille] ;
1522
1523 // Initialisation a zero :
1524 for (int i=0; i<(tb->dim).taille; i++) {
1525 xo[i] = 0 ;
1526 }
1527
1528 // On y va...
1529 double* xi = tb->t ;
1530 double* xci = xi ; // Pointeurs
1531 double* xco = xo ; // courants
1532
1533 double cx ;
1534
1535 // k = 0
1536 int m = 0 ; // Parite de l'harmonique en phi
1537
1538 // Dernier j: j = nt-1
1539 // Positionnement
1540 xci += nr * (nt-1) ;
1541 cx = -2 ;
1542 xco += nr * (nt-1) ;
1543
1544 // Initialisation de som
1545 for (int i=0 ; i<nr ; i++) {
1546 somP[i] = 0. ;
1547 somN[i] = 0. ;
1548 xco[i] = 0. ;
1549 }
1550
1551 // j suivants
1552 for (int j=nt-2 ; j >0 ; j--) {
1553 int l = j % 2 ;
1554 // Positionnement
1555 xci -= nr ;
1556 xco -= nr ;
1557 // Calcul
1558 for (int i=0 ; i<nr ; i++ ) {
1559 if(l==1) somN[i] += cx * xci[i] ;
1560 else somP[i] += cx * xci[i] ;
1561 xco[i] = somN[i]*(1-l)+somP[i]*l ;
1562 } // Fin de la boucle sur r
1563 } // Fin de la boucle sur theta
1564 // j = 0 sin(0*theta)
1565 xci -= nr ;
1566 xco -= nr ;
1567 for (int i=0 ; i<nr ; i++) {
1568 xco[i] = 0 ;
1569 }
1570 // Positionnement phi suivant
1571
1572 xci += nr*nt ;
1573 xco += nr*nt ;
1574
1575 // k=1
1576 xci += nr*nt ;
1577 xco += nr*nt ;
1578
1579 // k >= 2
1580 for (int k=2 ; k<np+1 ; k++) {
1581 m = (k/2) % 2 ; // Parite de l'harmonique en phi
1582
1583 // Dernier j: j = nt-1
1584 // Positionnement
1585 double fac_j0 = 0 ;
1586 if (m == 0) {
1587 cx = -2 ;
1588 }
1589 else {
1590 cx = 2 ;
1591 fac_j0 = 0.5 ;
1592 }
1593 xco += nr * (nt-1) ;
1594 xci += nr * (nt-1) ;
1595
1596 // Initialisation de som
1597 for (int i=0 ; i<nr ; i++) {
1598 somP[i] = 0. ;
1599 somN[i] = 0. ;
1600 xco[i] = 0. ;
1601 }
1602
1603 // j suivants
1604 for (int j=nt-2 ; j >= 0 ; j--) {
1605 int l = j % 2 ;
1606 // Positionnement
1607 xci -= nr ;
1608 xco -= nr ;
1609 // Calcul
1610 for (int i=0 ; i<nr ; i++ ) {
1611 if(l==1) somN[i] += cx * xci[i] ;
1612 else somP[i] += cx * xci[i] ;
1613 xco[i] = somN[i]*(1-l)+somP[i]*l ;
1614 } // Fin de la boucle sur r
1615 } // Fin de la boucle sur theta
1616
1617 // Normalisation du premier theta dans le cas sin, mise a zero dans le cas cos
1618 for (int i=0 ; i<nr ; i++) {
1619 xco[i] *= fac_j0 ;
1620 }
1621
1622 // Positionnement phi suivant
1623
1624 xci += nr*nt ;
1625 xco += nr*nt ;
1626
1627 } // Fin de la boucle sur phi
1628
1629 // On remet les choses la ou il faut
1630 delete [] tb->t ;
1631 tb->t = xo ;
1632
1633 //Menage
1634 delete [] somP ;
1635 delete [] somN ;
1636
1637 // base de developpement
1638 int base_r = b & MSQ_R ;
1639 int base_p = b & MSQ_P ;
1640 switch(base_r){
1641 case(R_CHEBPI_P):
1642 b = R_CHEBPI_I | base_p | T_COSSIN_S ;
1643 break ;
1644 case(R_CHEBPI_I):
1645 b = R_CHEBPI_P | base_p | T_COSSIN_S ;
1646 break ;
1647 default:
1648 b = base_r | base_p | T_COSSIN_S ;
1649 break;
1650 }
1651}
1652
1653 //---------------------
1654 // cas T_COSSIN_S
1655 //----------------
1656
1657void _ssint_t_cossin_s(Tbl* tb, int & b)
1658{
1659
1660
1661 // Peut-etre rien a faire ?
1662 if (tb->get_etat() == ETATZERO) {
1663 int base_r = b & MSQ_R ;
1664 int base_p = b & MSQ_P ;
1665 switch(base_r){
1666 case(R_CHEBPI_P):
1667 b = R_CHEBPI_I | base_p | T_COSSIN_C ;
1668 break ;
1669 case(R_CHEBPI_I):
1670 b = R_CHEBPI_P | base_p | T_COSSIN_C ;
1671 break ;
1672 default:
1673 b = base_r | base_p | T_COSSIN_C ;
1674 break;
1675 }
1676 return ;
1677 }
1678
1679 // Protection
1680 assert(tb->get_etat() == ETATQCQ) ;
1681
1682 // Pour le confort : nbre de points reels.
1683 int nr = (tb->dim).dim[0] ;
1684 int nt = (tb->dim).dim[1] ;
1685 int np = (tb->dim).dim[2] ;
1686 np = np - 2 ;
1687
1688 // zone de travail
1689 double* somP = new double [nr] ;
1690 double* somN = new double [nr] ;
1691
1692 // pt. sur le tableau de double resultat
1693 double* xo = new double[(tb->dim).taille] ;
1694
1695 // Initialisation a zero :
1696 for (int i=0; i<(tb->dim).taille; i++) {
1697 xo[i] = 0 ;
1698 }
1699
1700 // On y va...
1701 double* xi = tb->t ;
1702 double* xci = xi ; // Pointeurs
1703 double* xco = xo ; // courants
1704
1705 double cx ;
1706
1707 // k = 0
1708 int m = 0 ; // Parite de l'harmonique en phi
1709
1710 // Dernier j: j = nt-1
1711 // Positionnement
1712 xci += nr * (nt-1) ;
1713 cx = 2 ;
1714 xco += nr * (nt-1) ;
1715
1716 // Initialisation de som
1717 for (int i=0 ; i<nr ; i++) {
1718 somP[i] = 0. ;
1719 somN[i] = 0. ;
1720 xco[i] = 0. ;
1721 }
1722
1723 // j suivants
1724 for (int j=nt-2 ; j >= 0 ; j--) {
1725 int l = j % 2 ;
1726 // Positionnement
1727 xci -= nr ;
1728 xco -= nr ;
1729 // Calcul
1730 for (int i=0 ; i<nr ; i++ ) {
1731 if(l==1) somN[i] += cx * xci[i] ;
1732 else somP[i] += cx * xci[i] ;
1733 xco[i] = somN[i]*(1-l)+somP[i]*l ;
1734 } // Fin de la boucle sur r
1735 } // Fin de la boucle sur theta
1736 for (int i=0 ; i<nr ; i++) {
1737 xco[i] *= .5 ;
1738 }
1739 // Positionnement phi suivant
1740 xci += nr*nt ;
1741 xco += nr*nt ;
1742
1743 // k=1
1744 xci += nr*nt ;
1745 xco += nr*nt ;
1746
1747 // k >= 2
1748 for (int k=2 ; k<np+1 ; k++) {
1749 m = (k/2) % 2 ; // Parite de l'harmonique en phi
1750
1751 // Dernier j: j = nt-1
1752 // Positionnement
1753 if (m == 0) {
1754 xci += nr * (nt-1) ;
1755 cx = 2 ;
1756 }
1757 else {
1758 xci += nr * (nt-1) ;
1759 cx = -2 ;
1760 }
1761 xco += nr * (nt-1) ;
1762
1763 // Initialisation de som
1764 for (int i=0 ; i<nr ; i++) {
1765 somP[i] = 0. ;
1766 somN[i] = 0. ;
1767 xco[i] = 0. ;
1768 }
1769
1770 // j suivants
1771 for (int j=nt-2 ; j >= 0 ; j--) {
1772 int l = j % 2 ;
1773 // Positionnement
1774 xci -= nr ;
1775 xco -= nr ;
1776 // Calcul
1777 for (int i=0 ; i<nr ; i++ ) {
1778 if(l==1) somN[i] += cx * xci[i] ;
1779 else somP[i] += cx * xci[i] ;
1780 xco[i] = somN[i]*(1-l)+somP[i]*l ;
1781 } // Fin de la boucle sur r
1782 } // Fin de la boucle sur theta
1783
1784 // Normalisation du premier theta
1785 for (int i=0 ; i<nr ; i++) {
1786 xco[i] *= .5 ;
1787 }
1788
1789 // Positionnement phi suivant
1790
1791 xci += nr*nt ;
1792 xco += nr*nt ;
1793
1794 } // Fin de la boucle sur phi
1795
1796 // On remet les choses la ou il faut
1797 delete [] tb->t ;
1798 tb->t = xo ;
1799
1800 //Menage
1801 delete [] somP ;
1802 delete [] somN ;
1803
1804 // base de developpement
1805 int base_r = b & MSQ_R ;
1806 int base_p = b & MSQ_P ;
1807 switch(base_r){
1808 case(R_CHEBPI_P):
1809 b = R_CHEBPI_I | base_p | T_COSSIN_C ;
1810 break ;
1811 case(R_CHEBPI_I):
1812 b = R_CHEBPI_P | base_p | T_COSSIN_C ;
1813 break ;
1814 default:
1815 b = base_r | base_p | T_COSSIN_C ;
1816 break;
1817 }
1818}
1819}
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
R_CHEBPI_I
#define R_CHEBPI_I
Cheb. pair-impair suivant l impair pour l=0.
Definition type_parite.h:174
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
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