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