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