LORENE
som_symy.C
1/*
2 * Copyright (c) 2000-2001 Eric Gourgoulhon
3 *
4 * This file is part of LORENE.
5 *
6 * LORENE is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * LORENE is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with LORENE; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 *
20 */
21
22
23char som_symy_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/som_symy.C,v 1.3 2014/10/13 08:53:27 j_novak Exp $" ;
24
25/*
26 * Ensemble des routine pour la sommation directe en r, theta et phi dans
27 * le cas symetrie equatoriale (plan z=0) + symetrie par rapport au plan y=0
28 *
29 * SYNOPSYS:
30 * double som_r_XX_symy
31 * (double* ti, int nr, int nt, int np, double x, double* trtp)
32 *
33 * x est l'argument du polynome de Chebychev: x in [0, 1] ou x in [-1, 1]
34 *
35 *
36 * double som_tet_XX_symy
37 * (double* ti, int nt, int np, double tet, double* to)
38 *
39 * double som_phi_XX_symy
40 * (double* ti, int np, double phi, double* xo)
41 *
42 * ATTENTION: np est suppose etre le nombre de points (sans le +2)
43 * on suppose que trtp tient compte de ca.
44 */
45
46/*
47 * $Id: som_symy.C,v 1.3 2014/10/13 08:53:27 j_novak Exp $
48 * $Log: som_symy.C,v $
49 * Revision 1.3 2014/10/13 08:53:27 j_novak
50 * Lorene classes and functions now belong to the namespace Lorene.
51 *
52 * Revision 1.2 2014/10/06 15:16:07 j_novak
53 * Modified #include directives to use c++ syntax.
54 *
55 * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
56 * LORENE
57 *
58 * Revision 2.0 2000/03/06 10:27:53 eric
59 * *** empty log message ***
60 *
61 *
62 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/som_symy.C,v 1.3 2014/10/13 08:53:27 j_novak Exp $
63 *
64 */
65
66
67// Headers C
68#include <cmath>
69
70namespace Lorene {
71
72//****************************************************************************
73// Sommation en r
74//****************************************************************************
75
76
80
82 (double* ti, const int nr, const int nt, const int np, const double x,
83 double* trtp) {
84
85// Variables diverses
86int i, j, k ;
87double* pi = ti ; // pointeur courant sur l'entree
88double* po = trtp ; // pointeur courant sur la sortie
89
90 // Valeurs des polynomes de Chebyshev au point x demande
91 double* cheb = new double [nr] ;
92 cheb[0] = 1. ;
93 cheb[1] = x ;
94 for (i=2; i<nr; i++) {
95 cheb[i] = 2*x* cheb[i-1] - cheb[i-2] ;
96 }
97
98 // Sommation pour le premier phi, k=0
99 for (j=0 ; j<nt ; j++) {
100 *po = 0 ;
101 // Sommation sur r
102 for (i=0 ; i<nr ; i++) {
103 *po += (*pi) * cheb[i] ;
104 pi++ ; // R suivant
105 }
106 po++ ; // Theta suivant
107 }
108
109 if (np > 1) {
110
111 // On saute le deuxieme phi (sin(0)), k=1
112 pi += nr*nt ;
113 po += nt ;
114
115 // Sommation sur les phi suivants (pour k=2,...,np)
116 // en sautant les sinus (d'ou le k+=2)
117 for (k=2 ; k<np+1 ; k+=2) {
118 for (j=0 ; j<nt ; j++) {
119 // Sommation sur r
120 *po = 0 ;
121 for (i=0 ; i<nr ; i++) {
122 *po += (*pi) * cheb[i] ;
123 pi++ ; // R suivant
124 }
125 po++ ; // Theta suivant
126 }
127 // On saute le sin(k*phi) :
128 pi += nr*nt ;
129 po += nt ;
130 }
131
132 } // fin du cas np > 1
133
134 // Menage
135 delete [] cheb ;
136}
137
138
139
143
145 (double* ti, const int nr, const int nt, const int np, const double x, double* trtp) {
146
147// Variables diverses
148int i, j, k ;
149double* pi = ti ; // pointeur courant sur l'entree
150double* po = trtp ; // pointeur courant sur la sortie
151
152 // Valeurs des polynomes de Chebyshev au point x demande
153 double* cheb = new double [nr] ;
154 cheb[0] = 1. ;
155 cheb[1] = x ;
156 for (i=2; i<nr; i++) {
157 cheb[i] = 2*x* cheb[i-1] - cheb[i-2] ;
158 }
159
160 // Sommation pour le premier phi, k=0
161 for (j=0 ; j<nt ; j++) {
162 *po = 0 ;
163 // Sommation sur r
164 for (i=0 ; i<nr ; i++) {
165 *po += (*pi) * cheb[i] ;
166 pi++ ; // R suivant
167 }
168 po++ ; // Theta suivant
169 }
170
171 if (np > 1) {
172
173 // On saute le deuxieme phi (sin(0)), k=1
174 pi += nr*nt ;
175 po += nt ;
176
177 // Sommation sur les phi suivants (pour k=2,...,np)
178 // en sautant les sinus (d'ou le k+=2)
179 for (k=2 ; k<np+1 ; k+=2) {
180 for (j=0 ; j<nt ; j++) {
181 // Sommation sur r
182 *po = 0 ;
183 for (i=0 ; i<nr ; i++) {
184 *po += (*pi) * cheb[i] ;
185 pi++ ; // R suivant
186 }
187 po++ ; // Theta suivant
188 }
189 // On saute le sin(k*phi) :
190 pi += nr*nt ;
191 po += nt ;
192 }
193
194 } // fin du cas np > 1
195
196 // Menage
197 delete [] cheb ;
198
199}
200
204
206 (double* ti, const int nr, const int nt, const int np, const double x, double* trtp) {
207
208// Variables diverses
209int i, j, k ;
210double* pi = ti ; // pointeur courant sur l'entree
211double* po = trtp ; // pointeur courant sur la sortie
212
213 // Valeurs des polynomes de Chebyshev au point x demande
214 double* cheb = new double [2*nr] ;
215 cheb[0] = 1. ;
216 cheb[1] = x ;
217 for (i=2 ; i<2*nr ; i++) {
218 cheb[i] = 2*x* cheb[i-1] - cheb[i-2] ;
219 }
220
221 // Sommation pour le premier phi, k=0
222 int m = 0;
223 for (j=0 ; j<nt ; j++) {
224 *po = 0 ;
225 // Sommation sur r
226 for (i=0 ; i<nr ; i++) {
227 *po += (*pi) * cheb[2*i] ;
228 pi++ ; // R suivant
229 }
230 po++ ; // Theta suivant
231 }
232
233 if (np > 1) {
234
235 // On saute le deuxieme phi (sin(0)), k=1
236 pi += nr*nt ;
237 po += nt ;
238
239 // Sommation sur les phi suivants (pour k=2,...,np)
240 // en sautant les sinus (d'ou le k+=2)
241 for (k=2 ; k<np+1 ; k+=2) {
242 m = (k/2) % 2 ; // parite: 0 <-> 1
243 for (j=0 ; j<nt ; j++) {
244 // Sommation sur r
245 *po = 0 ;
246 for (i=0 ; i<nr ; i++) {
247 *po += (*pi) * cheb[2*i + m] ;
248 pi++ ; // R suivant
249 }
250 po++ ; // Theta suivant
251 }
252 // On saute le sin(k*phi) :
253 pi += nr*nt ;
254 po += nt ;
255 }
256
257 } // fin du cas np > 1
258
259 // Menage
260 delete [] cheb ;
261
262}
263
267
269 (double* ti, const int nr, const int nt, const int np, const double x, double* trtp) {
270
271// Variables diverses
272int i, j, k ;
273double* pi = ti ; // pointeur courant sur l'entree
274double* po = trtp ; // pointeur courant sur la sortie
275
276 // Valeurs des polynomes de Chebyshev au point x demande
277 double* cheb = new double [2*nr] ;
278 cheb[0] = 1. ;
279 cheb[1] = x ;
280 for (i=2 ; i<2*nr ; i++) {
281 cheb[i] = 2*x* cheb[i-1] - cheb[i-2] ;
282 }
283
284 // Sommation pour le premier phi, k=0
285 int m = 0;
286 for (j=0 ; j<nt ; j++) {
287 *po = 0 ;
288 // Sommation sur r
289 for (i=0 ; i<nr ; i++) {
290 *po += (*pi) * cheb[2*i+1] ;
291 pi++ ; // R suivant
292 }
293 po++ ; // Theta suivant
294 }
295
296 if (np > 1) {
297
298 // On saute le deuxieme phi (sin(0)), k=1
299 pi += nr*nt ;
300 po += nt ;
301
302 // Sommation sur les phi suivants (pour k=2,...,np)
303 // en sautant les sinus (d'ou le k+=2)
304 for (k=2 ; k<np+1 ; k+=2) {
305 m = (k/2) % 2 ; // parite: 0 <-> 1
306 for (j=0 ; j<nt ; j++) {
307 // Sommation sur r
308 *po = 0 ;
309 for (i=0 ; i<nr ; i++) {
310 *po += (*pi) * cheb[2*i + 1 - m] ;
311 pi++ ; // R suivant
312 }
313 po++ ; // Theta suivant
314 }
315 // On saute le sin(k*phi) :
316 pi += nr*nt ;
317 po += nt ;
318 }
319
320 } // fin du cas np > 1
321
322 // Menage
323 delete [] cheb ;
324
325}
326
327
328
329//****************************************************************************
330// Sommation en theta
331//****************************************************************************
332
336
338 (double* ti, const int nt, const int np,
339 const double tet, double* to) {
340
341// Variables diverses
342int j, k ;
343double* pi = ti ; // Pointeur courant sur l'entree
344double* po = to ; // Pointeur courant sur la sortie
345
346 // Initialisation des tables trigo
347 double* cossin = new double [2*nt] ;
348 for (j=0 ; j<2*nt ; j +=2) {
349 cossin[j] = cos(j * tet) ;
350 cossin[j+1] = sin((j+1) * tet) ;
351 }
352
353 // Sommation sur le premier phi -> cosinus, k=0
354 *po = 0 ;
355 for (j=0 ; j<nt ; j++) {
356 *po += (*pi) * cossin[2*j] ;
357 pi++ ; // Theta suivant
358 }
359 po++ ; // Phi suivant
360
361 if (np > 1) {
362
363 // On saute le phi suivant (sin(0)), k=1
364 pi += nt ;
365 po++ ;
366
367 // Sommation sur le reste des phi (pour k=2,...,np), suivant parite de m
368 for (k=2 ; k<np+1 ; k+=2) {
369 int m = (k/2) % 2 ; // parite: 0 <-> 1
370 (*po) = 0 ;
371 for (j=0 ; j<nt ; j++) {
372 *po += (*pi) * cossin[2*j + m] ;
373 pi++ ; // Theta suivant
374 }
375 po++ ; // Phi suivant
376
377 // On saute le sin(k*phi) :
378 pi += nt ;
379 po++ ;
380
381 }
382 } // fin du cas np > 1
383
384 // Menage
385 delete [] cossin ;
386}
387
391
393 (double* ti, const int nt, const int np,
394 const double tet, double* to) {
395
396// Variables diverses
397int j, k ;
398double* pi = ti ; // Pointeur courant sur l'entree
399double* po = to ; // Pointeur courant sur la sortie
400
401 // Initialisation des tables trigo
402 double* cossin = new double [2*nt] ;
403 for (j=0 ; j<2*nt ; j +=2) {
404 cossin[j] = cos((j+1) * tet) ;
405 cossin[j+1] = sin(j * tet) ;
406 }
407
408 // Sommation sur le premier phi -> cosinus, k=0
409 *po = 0 ;
410 for (j=0 ; j<nt ; j++) {
411 *po += (*pi) * cossin[2*j] ;
412 pi++ ; // Theta suivant
413 }
414 po++ ; // Phi suivant
415
416 if (np > 1) {
417
418 // On saute le phi suivant (sin(0)), k=1
419 pi += nt ;
420 po++ ;
421
422 // Sommation sur le reste des phi (pour k=2,...,np), suivant parite de m
423 for (k=2 ; k<np+1 ; k+=2) {
424 int m = (k/2) % 2 ; // parite: 0 <-> 1
425 (*po) = 0 ;
426 for (j=0 ; j<nt ; j++) {
427 *po += (*pi) * cossin[2*j + m] ;
428 pi++ ; // Theta suivant
429 }
430 po++ ; // Phi suivant
431
432 // On saute le sin(k*phi) :
433 pi += nt ;
434 po++ ;
435
436 }
437 } // fin du cas np > 1
438
439 // Menage
440 delete [] cossin ;
441}
442
443
444//****************************************************************************
445// Sommation en phi
446//****************************************************************************
447
448void som_phi_cossin_symy
449 (double* ti, const int np, const double phi, double* xo) {
450
451 *xo = ti[0] ; // premier element
452
453 // Sommation sur les cosinus seulement
454 for (int k=2 ; k<np-1 ; k +=2 ) {
455 int m = k/2 ;
456 *xo += ti[k] * cos(m * phi) ;
457 }
458 *xo += ti[np] * cos(np/2 * phi) ;
459}
460
461}
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
void som_r_cheb_symy(double *, const int, const int, const int, const double, double *)
Cas R_CHEB ///.
Definition som_symy.C:82
void som_tet_cossin_ci_symy(double *, const int, const int, const double, double *)
Cas T_COSSIN_CI ///.
Definition som_symy.C:393
void som_r_chebu_symy(double *, const int, const int, const int, const double, double *)
Cas R_CHEBU ///.
Definition som_symy.C:145
void som_r_chebpim_i_symy(double *, const int, const int, const int, const double, double *)
Cas R_CHEBPIM_I ///.
Definition som_symy.C:269
void som_tet_cossin_cp_symy(double *, const int, const int, const double, double *)
Cas T_COSSIN_CP ///.
Definition som_symy.C:338
void som_r_chebpim_p_symy(double *, const int, const int, const int, const double, double *)
Cas R_CHEBPIM_P ///.
Definition som_symy.C:206