LORENE
som_tet.C
1 /*
2  * Copyright (c) 1999-2000 Jean-Alain Marck
3  * Copyright (c) 1999-2001 Philippe Grandclement
4  * Copyright (c) 1999-2002 Eric Gourgoulhon
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 som_tet_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/som_tet.C,v 1.8 2014/10/13 08:53:27 j_novak Exp $" ;
26 
27 /*
28  * Ensemble des routine pour la sommation directe en theta
29  *
30  * SYNOPSYS:
31  * double som_tet_XX
32  * (double* ti, int nt, int np, double tet, double* to)
33  *
34  * ATTENTION: np est le vrai nombre de points (sans le +2)
35  * on suppose que les tableaux tiennent compte de ca.
36  */
37 
38 
39 /*
40  * $Id: som_tet.C,v 1.8 2014/10/13 08:53:27 j_novak Exp $
41  * $Log: som_tet.C,v $
42  * Revision 1.8 2014/10/13 08:53:27 j_novak
43  * Lorene classes and functions now belong to the namespace Lorene.
44  *
45  * Revision 1.7 2014/10/06 15:16:07 j_novak
46  * Modified #include directives to use c++ syntax.
47  *
48  * Revision 1.6 2004/11/23 15:16:02 m_forot
49  *
50  * Added the bases for the cases without any equatorial symmetry
51  * (T_COSSIN_C, T_COSSIN_S, T_LEG, R_CHEBPI_P, R_CHEBPI_I).
52  *
53  * Revision 1.5 2002/10/16 14:36:59 j_novak
54  * Reorganization of #include instructions of standard C++, in order to
55  * use experimental version 3 of gcc.
56  *
57  * Revision 1.4 2002/05/11 12:36:54 e_gourgoulhon
58  * Added basis T_COSSIN_SI.
59  *
60  * Revision 1.3 2002/05/05 16:20:40 e_gourgoulhon
61  * Error message (in unknwon basis case) in English
62  * Added the basis T_COSSIN_SP
63  *
64  * Revision 1.2 2002/05/01 07:41:05 e_gourgoulhon
65  * Correction of an ERROR in som_tet_sin_p :
66  * sin(2*(j+1) * tet) --> sin(2*j * tet)
67  * idem in som_tet_sin:
68  * sin( (j+1) * tet) --> sin(j * tet)
69  *
70  * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
71  * LORENE
72  *
73  * Revision 2.5 2000/09/08 16:26:32 eric
74  * Ajout de la base T_SIN_I.
75  *
76  * Revision 2.4 2000/09/06 14:00:01 eric
77  * Ajout de la base T_COS_I.
78  *
79  * Revision 2.3 2000/03/06 09:34:47 eric
80  * Suppression des #include inutiles.
81  *
82  * Revision 2.2 1999/04/28 12:27:52 phil
83  * Correction tailles des tableaux
84  *
85  * Revision 2.1 1999/04/26 14:31:17 phil
86  * remplacement des malloc en new
87  *
88  * Revision 2.0 1999/04/12 15:42:03 phil
89  * *** empty log message ***
90  *
91  * Revision 1.1 1999/04/12 15:41:20 phil
92  * Initial revision
93  *
94  *
95  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/som_tet.C,v 1.8 2014/10/13 08:53:27 j_novak Exp $
96  *
97  */
98 // Headers C
99 #include <cstdlib>
100 #include <cmath>
101 
102 #include "headcpp.h"
103 
104 namespace Lorene {
105 
106  //--------------------
107  //- Cas Non-Prevu ---
108  //------------------
109 
110 void som_tet_pas_prevu
111  (double* , const int, const int, const double, double*) {
112  cout << "Mtbl_cf::val_point: theta basis not implemented yet ! "
113  << endl ;
114  abort () ;
115 }
116 
117  //----------------
118  // Cas T_COS ---
119  //----------------
120 
121 void som_tet_cos
122  (double* ti, const int nt, const int np,
123  const double tet, double* to) {
124 
125 // Variables diverses
126 int j, k ;
127 double* pi = ti ; // Pointeur courant sur l'entree
128 double* po = to ; // Pointeur courant sur la sortie
129 
130  // Initialisation des tables trigo
131  double* cosinus = new double [nt] ;
132  for (j=0 ; j<nt ; j++) {
133  cosinus[j] = cos(j * tet) ;
134  }
135 
136  // Sommation sur le premier phi, k=0
137  *po = 0 ;
138  for (j=0 ; j<nt ; j++) {
139  *po += (*pi) * cosinus[j] ;
140  pi++ ; // Theta suivant
141  }
142  po++ ; // Phi suivant
143 
144  if (np > 1) {
145 
146  // On saute le phi suivant (sin(0)), k=1
147  pi += nt ;
148  po++ ;
149 
150  // Sommation sur le reste des phi (pour k=2,...,np)
151  for (k=2 ; k<np+1 ; k++) {
152  (*po) = 0 ;
153  for (j=0 ; j<nt ; j++) {
154  *po += (*pi) * cosinus[j] ;
155  pi++ ; // Theta suivant
156  }
157  po++ ; // Phi suivant
158  }
159 
160  } // fin du cas np > 1
161 
162  // Menage
163  delete [] cosinus ;
164 
165 }
166 
167  //-------------------
168  // Cas T_COS_P ---
169  //------------------
170 
171 void som_tet_cos_p
172  (double* ti, const int nt, const int np,
173  const double tet, double* to) {
174 
175 // Variables diverses
176 int j, k ;
177 double* pi = ti ; // Pointeur courant sur l'entree
178 double* po = to ; // Pointeur courant sur la sortie
179 
180  // Initialisation des tables trigo
181  double* cosinus = new double [nt] ;
182  for (j=0 ; j<nt ; j++) {
183  cosinus[j] = cos(2*j * tet) ;
184  }
185 
186  // Sommation sur le premier phi, k=0
187  *po = 0 ;
188  for (j=0 ; j<nt ; j++) {
189  *po += (*pi) * cosinus[j] ;
190  pi++ ; // Theta suivant
191  }
192  po++ ; // Phi suivant
193 
194  if (np > 1) {
195 
196  // On saute le phi suivant (sin(0)), k=1
197  pi += nt ;
198  po++ ;
199 
200  // Sommation sur le reste des phi (pour k=2,...,np)
201  for (k=2 ; k<np+1 ; k++) {
202  (*po) = 0 ;
203  for (j=0 ; j<nt ; j++) {
204  *po += (*pi) * cosinus[j] ;
205  pi++ ; // Theta suivant
206  }
207  po++ ; // Phi suivant
208  }
209 
210  } // fin du cas np > 1
211  // Menage
212  delete [] cosinus ;
213 
214 }
215 
216  //----------------------
217  //- Cas T_COS_I ---
218  //---------------------
219 
220 void som_tet_cos_i
221  (double* ti, const int nt, const int np,
222  const double tet, double* to) {
223 
224 // Variables diverses
225 int j, k ;
226 double* pi = ti ; // Pointeur courant sur l'entree
227 double* po = to ; // Pointeur courant sur la sortie
228 
229  // Initialisation des tables trigo
230  double* cosinus = new double [nt] ;
231  for (j=0 ; j<nt-1 ; j++) {
232  cosinus[j] = cos((2*j+1) * tet) ;
233  }
234  cosinus[nt-1] = 0 ;
235 
236  // Sommation sur le premier phi, k=0
237  *po = 0 ;
238  for (j=0 ; j<nt ; j++) {
239  *po += (*pi) * cosinus[j] ;
240  pi++ ; // Theta suivant
241  }
242  po++ ; // Phi suivant
243 
244  if (np > 1) {
245 
246  // On saute le phi suivant (sin(0)), k=1
247  pi += nt ;
248  po++ ;
249 
250  // Sommation sur le reste des phi (pour k=2,...,np)
251  for (k=2 ; k<np+1 ; k++) {
252  (*po) = 0 ;
253  for (j=0 ; j<nt ; j++) {
254  *po += (*pi) * cosinus[j] ;
255  pi++ ; // Theta suivant
256  }
257  po++ ; // Phi suivant
258  }
259  } // fin du cas np > 1
260 
261  // Menage
262  delete [] cosinus ;
263 
264 }
265 
266 
267 
268 
269 
270  //---------------
271  // Cas T_SIN ---
272  //---------------
273 
274 void som_tet_sin
275  (double* ti, const int nt, const int np,
276  const double tet, double* to) {
277 
278 // Variables diverses
279 int j, k ;
280 double* pi = ti ; // Pointeur courant sur l'entree
281 double* po = to ; // Pointeur courant sur la sortie
282 
283  // Initialisation des tables trigo
284  double* sinus = new double [nt] ;
285  for (j=0 ; j<nt ; j++) {
286  sinus[j] = sin(j * tet) ;
287  }
288 
289  // Sommation sur le premier phi, k=0
290  *po = 0 ;
291  for (j=0 ; j<nt ; j++) {
292  *po += (*pi) * sinus[j] ;
293  pi++ ; // Theta suivant
294  }
295  po++ ; // Phi suivant
296 
297  if (np > 1) {
298 
299  // On saute le phi suivant (sin(0)), k=1
300  pi += nt ;
301  po++ ;
302 
303  // Sommation sur le reste des phi (pour k=2,...,np)
304  for (k=2 ; k<np+1 ; k++) {
305  (*po) = 0 ;
306  for (j=0 ; j<nt ; j++) {
307  *po += (*pi) * sinus[j] ;
308  pi++ ; // Theta suivant
309  }
310  po++ ; // Phi suivant
311  }
312  } // fin du cas np > 1
313 
314  // Menage
315  delete [] sinus ;
316 
317 }
318 
319  //------------------
320  // Cas T_SIN_P ---
321  //-----------------
322 
323 void som_tet_sin_p
324  (double* ti, const int nt, const int np,
325  const double tet, double* to) {
326 
327 // Variables diverses
328 int j, k ;
329 double* pi = ti ; // Pointeur courant sur l'entree
330 double* po = to ; // Pointeur courant sur la sortie
331 
332  // Initialisation des tables trigo
333  double* sinus = new double [nt] ;
334  for (j=0 ; j<nt-1 ; j++) {
335  sinus[j] = sin(2*j * tet) ;
336  }
337  sinus[nt-1] = 0 ;
338 
339  // Sommation sur le premier phi, k=0
340  *po = 0 ;
341  for (j=0 ; j<nt ; j++) {
342  *po += (*pi) * sinus[j] ;
343  pi++ ; // Theta suivant
344  }
345  po++ ; // Phi suivant
346 
347  if (np > 1) {
348 
349  // On saute le phi suivant (sin(0)), k=1
350  pi += nt ;
351  po++ ;
352 
353  // Sommation sur le reste des phi (pour k=2,...,np)
354  for (k=2 ; k<np+1 ; k++) {
355  (*po) = 0 ;
356  for (j=0 ; j<nt ; j++) {
357  *po += (*pi) * sinus[j] ;
358  pi++ ; // Theta suivant
359  }
360  po++ ; // Phi suivant
361  }
362  } // fin du cas np > 1
363 
364  // Menage
365  delete [] sinus ;
366 
367 }
368 
369  //------------------
370  // Cas T_SIN_I ---
371  //-----------------
372 
373 void som_tet_sin_i
374  (double* ti, const int nt, const int np,
375  const double tet, double* to) {
376 
377 // Variables diverses
378 int j, k ;
379 double* pi = ti ; // Pointeur courant sur l'entree
380 double* po = to ; // Pointeur courant sur la sortie
381 
382  // Initialisation des tables trigo
383  double* sinus = new double [nt] ;
384  for (j=0 ; j<nt-1 ; j++) {
385  sinus[j] = sin( (2*j+1) * tet) ;
386  }
387  sinus[nt-1] = 0 ;
388 
389  // Sommation sur le premier phi, k=0
390  *po = 0 ;
391  for (j=0 ; j<nt ; j++) {
392  *po += (*pi) * sinus[j] ;
393  pi++ ; // Theta suivant
394  }
395  po++ ; // Phi suivant
396 
397  if (np > 1) {
398 
399  // On saute le phi suivant (sin(0)), k=1
400  pi += nt ;
401  po++ ;
402 
403  // Sommation sur le reste des phi (pour k=2,...,np)
404  for (k=2 ; k<np+1 ; k++) {
405  (*po) = 0 ;
406  for (j=0 ; j<nt ; j++) {
407  *po += (*pi) * sinus[j] ;
408  pi++ ; // Theta suivant
409  }
410  po++ ; // Phi suivant
411  }
412  } // fin du cas np > 1
413 
414  // Menage
415  delete [] sinus ;
416 
417 }
418 
419 
420  //---------------------
421  // Cas T_COSSIN_CP ---
422  //---------------------
423 
424 void som_tet_cossin_cp
425  (double* ti, const int nt, const int np,
426  const double tet, double* to) {
427 
428 // Variables diverses
429 int j, k ;
430 double* pi = ti ; // Pointeur courant sur l'entree
431 double* po = to ; // Pointeur courant sur la sortie
432 
433  // Initialisation des tables trigo
434  double* cossin = new double [2*nt] ;
435  for (j=0 ; j<2*nt ; j +=2) {
436  cossin[j] = cos(j * tet) ;
437  cossin[j+1] = sin((j+1) * tet) ;
438  }
439 
440  // Sommation sur le premier phi -> cosinus, k=0
441  *po = 0 ;
442  for (j=0 ; j<nt ; j++) {
443  *po += (*pi) * cossin[2*j] ;
444  pi++ ; // Theta suivant
445  }
446  po++ ; // Phi suivant
447 
448  if (np > 1) {
449 
450  // On saute le phi suivant (sin(0)), k=1
451  pi += nt ;
452  po++ ;
453 
454  // Sommation sur le reste des phi (pour k=2,...,np), suivant parite de m
455  for (k=2 ; k<np+1 ; k++) {
456  int m = (k/2) % 2 ; // parite: 0 <-> 1
457  (*po) = 0 ;
458  for (j=0 ; j<nt ; j++) {
459  *po += (*pi) * cossin[2*j + m] ;
460  pi++ ; // Theta suivant
461  }
462  po++ ; // Phi suivant
463  }
464  } // fin du cas np > 1
465 
466  // Menage
467  delete [] cossin ;
468 
469 }
470 
471 
472  //----------------------
473  //- Cas T_COSSIN_CI ---
474  //---------------------
475 
476 void som_tet_cossin_ci
477  (double* ti, const int nt, const int np,
478  const double tet, double* to) {
479 
480 // Variables diverses
481 int j, k ;
482 double* pi = ti ; // Pointeur courant sur l'entree
483 double* po = to ; // Pointeur courant sur la sortie
484 
485  // Initialisation des tables trigo
486  double* cossin = new double [2*nt] ;
487  for (j=0 ; j<2*nt ; j +=2) {
488  cossin[j] = cos((j+1) * tet) ;
489  cossin[j+1] = sin(j * tet) ;
490  }
491 
492  // Sommation sur le premier phi -> cosinus, k=0
493  *po = 0 ;
494  for (j=0 ; j<nt ; j++) {
495  *po += (*pi) * cossin[2*j] ;
496  pi++ ; // Theta suivant
497  }
498  po++ ; // Phi suivant
499 
500  if (np > 1) {
501 
502  // On saute le phi suivant (sin(0)), k=1
503  pi += nt ;
504  po++ ;
505 
506  // Sommation sur le reste des phi (pour k=2,...,np), suivant parite de m
507  for (k=2 ; k<np+1 ; k++) {
508  int m = (k/2) % 2 ; // parite: 0 <-> 1
509  (*po) = 0 ;
510  for (j=0 ; j<nt ; j++) {
511  *po += (*pi) * cossin[2*j + m] ;
512  pi++ ; // Theta suivant
513  }
514  po++ ; // Phi suivant
515  }
516  } // fin du cas np > 1
517 
518  // Menage
519  delete [] cossin ;
520 
521 }
522 
523 
524  //---------------------
525  // Cas T_COSSIN_SP ---
526  //---------------------
527 
528 void som_tet_cossin_sp
529  (double* ti, const int nt, const int np,
530  const double tet, double* to) {
531 
532 // Variables diverses
533 int j, k ;
534 double* pi = ti ; // Pointeur courant sur l'entree
535 double* po = to ; // Pointeur courant sur la sortie
536 
537  // Initialisation des tables trigo
538  double* cossin = new double [2*nt] ;
539  for (j=0 ; j<2*nt ; j +=2) {
540  cossin[j] = sin(j * tet) ;
541  cossin[j+1] = cos((j+1) * tet) ;
542  }
543 
544  // Sommation sur le premier phi -> cosinus, k=0
545  *po = 0 ;
546  for (j=0 ; j<nt ; j++) {
547  *po += (*pi) * cossin[2*j] ;
548  pi++ ; // Theta suivant
549  }
550  po++ ; // Phi suivant
551 
552  if (np > 1) {
553 
554  // On saute le phi suivant (sin(0)), k=1
555  pi += nt ;
556  po++ ;
557 
558  // Sommation sur le reste des phi (pour k=2,...,np), suivant parite de m
559  for (k=2 ; k<np+1 ; k++) {
560  int m = (k/2) % 2 ; // parite: 0 <-> 1
561  (*po) = 0 ;
562  for (j=0 ; j<nt ; j++) {
563  *po += (*pi) * cossin[2*j + m] ;
564  pi++ ; // Theta suivant
565  }
566  po++ ; // Phi suivant
567  }
568  } // fin du cas np > 1
569 
570  // Menage
571  delete [] cossin ;
572 
573 }
574 
575 
576  //---------------------
577  // Cas T_COSSIN_SI ---
578  //---------------------
579 
580 void som_tet_cossin_si
581  (double* ti, const int nt, const int np,
582  const double tet, double* to) {
583 
584 // Variables diverses
585 int j, k ;
586 double* pi = ti ; // Pointeur courant sur l'entree
587 double* po = to ; // Pointeur courant sur la sortie
588 
589  // Initialisation des tables trigo
590  double* cossin = new double [2*nt] ;
591  for (j=0 ; j<2*nt ; j +=2) {
592  cossin[j] = sin((j+1) * tet) ;
593  cossin[j+1] = cos(j * tet) ;
594  }
595 
596  // Sommation sur le premier phi -> cosinus, k=0
597  *po = 0 ;
598  for (j=0 ; j<nt ; j++) {
599  *po += (*pi) * cossin[2*j] ;
600  pi++ ; // Theta suivant
601  }
602  po++ ; // Phi suivant
603 
604  if (np > 1) {
605 
606  // On saute le phi suivant (sin(0)), k=1
607  pi += nt ;
608  po++ ;
609 
610  // Sommation sur le reste des phi (pour k=2,...,np), suivant parite de m
611  for (k=2 ; k<np+1 ; k++) {
612  int m = (k/2) % 2 ; // parite: 0 <-> 1
613  (*po) = 0 ;
614  for (j=0 ; j<nt ; j++) {
615  *po += (*pi) * cossin[2*j + m] ;
616  pi++ ; // Theta suivant
617  }
618  po++ ; // Phi suivant
619  }
620  } // fin du cas np > 1
621 
622  // Menage
623  delete [] cossin ;
624 
625 }
626 
627  //---------------------
628  // Cas T_COSSIN_C ---
629  //---------------------
630 
631 void som_tet_cossin_c
632  (double* ti, const int nt, const int np,
633  const double tet, double* to) {
634 
635 // Variables diverses
636 int j, k ;
637 double* pi = ti ; // Pointeur courant sur l'entree
638 double* po = to ; // Pointeur courant sur la sortie
639 
640  // Initialisation des tables trigo
641  double* cossin = new double [2*nt] ;
642  for (j=0 ; j<2*nt ; j +=2) {
643  cossin[j] = cos(j/2 * tet) ;
644  cossin[j+1] = sin(j/2 * tet) ;
645  }
646 
647  // Sommation sur le premier phi -> cosinus, k=0
648  *po = 0 ;
649  for (j=0 ; j<nt ; j++) {
650  *po += (*pi) * cossin[2*j] ;
651  pi++ ; // Theta suivant
652  }
653  po++ ; // Phi suivant
654 
655  if (np > 1) {
656 
657  // On saute le phi suivant (sin(0)), k=1
658  pi += nt ;
659  po++ ;
660 
661  // Sommation sur le reste des phi (pour k=2,...,np), suivant parite de m
662  for (k=2 ; k<np+1 ; k++) {
663  int m = (k/2) % 2 ; // parite: 0 <-> 1
664  (*po) = 0 ;
665  for (j=0 ; j<nt ; j++) {
666  *po += (*pi) * cossin[2*j + m] ;
667  pi++ ; // Theta suivant
668  }
669  po++ ; // Phi suivant
670  }
671  } // fin du cas np > 1
672 
673  // Menage
674  delete [] cossin ;
675 
676 }
677 
678  //---------------------
679  // Cas T_COSSIN_S ---
680  //---------------------
681 
682 void som_tet_cossin_s
683  (double* ti, const int nt, const int np,
684  const double tet, double* to) {
685 
686 // Variables diverses
687 int j, k ;
688 double* pi = ti ; // Pointeur courant sur l'entree
689 double* po = to ; // Pointeur courant sur la sortie
690 
691  // Initialisation des tables trigo
692  double* cossin = new double [2*nt] ;
693  for (j=0 ; j<2*nt ; j +=2) {
694  cossin[j] = sin(j/2 * tet) ;
695  cossin[j+1] = cos(j/2 * tet) ;
696  }
697 
698  // Sommation sur le premier phi -> cosinus, k=0
699  *po = 0 ;
700  for (j=0 ; j<nt ; j++) {
701  *po += (*pi) * cossin[2*j] ;
702  pi++ ; // Theta suivant
703  }
704  po++ ; // Phi suivant
705 
706  if (np > 1) {
707 
708  // On saute le phi suivant (sin(0)), k=1
709  pi += nt ;
710  po++ ;
711 
712  // Sommation sur le reste des phi (pour k=2,...,np), suivant parite de m
713  for (k=2 ; k<np+1 ; k++) {
714  int m = (k/2) % 2 ; // parite: 0 <-> 1
715  (*po) = 0 ;
716  for (j=0 ; j<nt ; j++) {
717  *po += (*pi) * cossin[2*j + m] ;
718  pi++ ; // Theta suivant
719  }
720  po++ ; // Phi suivant
721  }
722  } // fin du cas np > 1
723 
724  // Menage
725  delete [] cossin ;
726 
727 }
728 
729 }
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::cos
Cmp cos(const Cmp &)
Cosine.
Definition: cmp_math.C:94
Lorene::sin
Cmp sin(const Cmp &)
Sine.
Definition: cmp_math.C:69