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