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 char 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 
83 namespace Lorene {
84 void _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 
95 void _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 
254 void _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 
408 void _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 
525 void _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 
649 void _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  //---------------------
774 void _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  //----------------------
890 void _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  //----------------
1028 void _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  //----------------
1181 void _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  //----------------
1323 void _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  //----------------------
1482 void _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 
1654 void _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
Definition: type_parite.h:210
#define T_SIN_P
dev. sin seulement, harmoniques paires
Definition: type_parite.h:202
#define MSQ_R
Extraction de l'info sur R.
Definition: type_parite.h:152
#define T_COSSIN_S
dev. cos-sin alternes, sin pour m=0
Definition: type_parite.h:194
#define T_COSSIN_SI
sin impair-cos pair alternes, sin pour m=0
Definition: type_parite.h:214
#define R_CHEBPI_I
Cheb. pair-impair suivant l impair pour l=0.
Definition: type_parite.h:174
#define T_COS_P
dev. cos seulement, harmoniques paires
Definition: type_parite.h:200
#define T_COSSIN_CI
cos impair-sin pair alternes, cos pour m=0
Definition: type_parite.h:212
#define T_COSSIN_CP
cos pair-sin impair alternes, cos pour m=0
Definition: type_parite.h:208
#define T_SIN_I
dev. sin seulement, harmoniques impaires
Definition: type_parite.h:206
#define T_COS
dev. cos seulement
Definition: type_parite.h:196
#define MSQ_P
Extraction de l'info sur Phi.
Definition: type_parite.h:156
#define T_SIN
dev. sin seulement
Definition: type_parite.h:198
#define T_COS_I
dev. cos seulement, harmoniques impaires
Definition: type_parite.h:204
#define T_COSSIN_C
dev. cos-sin alternes, cos pour m=0
Definition: type_parite.h:192
#define R_CHEBPI_P
Cheb. pair-impair suivant l pair pour l=0.
Definition: type_parite.h:172
Lorene prototypes.
Definition: app_hor.h:64