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 
86 namespace Lorene {
87 void _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 
98 void _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 
257 void _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 
411 void _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 
528 void _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 
652 void _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  //---------------------
777 void _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  //----------------------
893 void _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  //----------------
1031 void _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  //----------------
1184 void _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  //----------------
1326 void _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  //----------------------
1485 void _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 
1657 void _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 }
R_CHEBPI_I
#define R_CHEBPI_I
Cheb. pair-impair suivant l impair pour l=0.
Definition: type_parite.h:174
Lorene
Lorene prototypes.
Definition: app_hor.h:67
MSQ_P
#define MSQ_P
Extraction de l&#39;info sur Phi.
Definition: type_parite.h:156
T_COS
#define T_COS
dev. cos seulement
Definition: type_parite.h:196
T_COSSIN_SP
#define T_COSSIN_SP
sin pair-cos impair alternes, sin pour m=0
Definition: type_parite.h:210
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_COS_P
#define T_COS_P
dev. cos seulement, harmoniques paires
Definition: type_parite.h:200
T_COSSIN_C
#define T_COSSIN_C
dev. cos-sin alternes, cos pour m=0
Definition: type_parite.h:192
MSQ_R
#define MSQ_R
Extraction de l&#39;info sur R.
Definition: type_parite.h:152
T_SIN_P
#define T_SIN_P
dev. sin seulement, harmoniques paires
Definition: type_parite.h:202
T_COSSIN_SI
#define T_COSSIN_SI
sin impair-cos pair alternes, sin pour m=0
Definition: type_parite.h:214
R_CHEBPI_P
#define R_CHEBPI_P
Cheb. pair-impair suivant l pair pour l=0.
Definition: type_parite.h:172
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_COSSIN_S
#define T_COSSIN_S
dev. cos-sin alternes, sin pour m=0
Definition: type_parite.h:194