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
25char 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
104namespace Lorene {
105
106 //--------------------
107 //- Cas Non-Prevu ---
108 //------------------
109
110void 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
121void som_tet_cos
122 (double* ti, const int nt, const int np,
123 const double tet, double* to) {
124
125// Variables diverses
126int j, k ;
127double* pi = ti ; // Pointeur courant sur l'entree
128double* 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
171void som_tet_cos_p
172 (double* ti, const int nt, const int np,
173 const double tet, double* to) {
174
175// Variables diverses
176int j, k ;
177double* pi = ti ; // Pointeur courant sur l'entree
178double* 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
220void som_tet_cos_i
221 (double* ti, const int nt, const int np,
222 const double tet, double* to) {
223
224// Variables diverses
225int j, k ;
226double* pi = ti ; // Pointeur courant sur l'entree
227double* 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
274void som_tet_sin
275 (double* ti, const int nt, const int np,
276 const double tet, double* to) {
277
278// Variables diverses
279int j, k ;
280double* pi = ti ; // Pointeur courant sur l'entree
281double* 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
323void som_tet_sin_p
324 (double* ti, const int nt, const int np,
325 const double tet, double* to) {
326
327// Variables diverses
328int j, k ;
329double* pi = ti ; // Pointeur courant sur l'entree
330double* 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
373void som_tet_sin_i
374 (double* ti, const int nt, const int np,
375 const double tet, double* to) {
376
377// Variables diverses
378int j, k ;
379double* pi = ti ; // Pointeur courant sur l'entree
380double* 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
424void som_tet_cossin_cp
425 (double* ti, const int nt, const int np,
426 const double tet, double* to) {
427
428// Variables diverses
429int j, k ;
430double* pi = ti ; // Pointeur courant sur l'entree
431double* 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
476void som_tet_cossin_ci
477 (double* ti, const int nt, const int np,
478 const double tet, double* to) {
479
480// Variables diverses
481int j, k ;
482double* pi = ti ; // Pointeur courant sur l'entree
483double* 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
528void som_tet_cossin_sp
529 (double* ti, const int nt, const int np,
530 const double tet, double* to) {
531
532// Variables diverses
533int j, k ;
534double* pi = ti ; // Pointeur courant sur l'entree
535double* 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
580void som_tet_cossin_si
581 (double* ti, const int nt, const int np,
582 const double tet, double* to) {
583
584// Variables diverses
585int j, k ;
586double* pi = ti ; // Pointeur courant sur l'entree
587double* 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
631void som_tet_cossin_c
632 (double* ti, const int nt, const int np,
633 const double tet, double* to) {
634
635// Variables diverses
636int j, k ;
637double* pi = ti ; // Pointeur courant sur l'entree
638double* 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
682void som_tet_cossin_s
683 (double* ti, const int nt, const int np,
684 const double tet, double* to) {
685
686// Variables diverses
687int j, k ;
688double* pi = ti ; // Pointeur courant sur l'entree
689double* 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}
Cmp sin(const Cmp &)
Sine.
Definition cmp_math.C:69
Cmp cos(const Cmp &)
Cosine.
Definition cmp_math.C:94
Lorene prototypes.
Definition app_hor.h:64