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
23
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.5 2017/02/24 15:34:59 j_novak Exp $
48 * $Log: som_symy.C,v $
49 * Revision 1.5 2017/02/24 15:34:59 j_novak
50 * Removal of spurious comments
51 *
52 * Revision 1.4 2016/12/05 16:18:08 j_novak
53 * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
54 *
55 * Revision 1.3 2014/10/13 08:53:27 j_novak
56 * Lorene classes and functions now belong to the namespace Lorene.
57 *
58 * Revision 1.2 2014/10/06 15:16:07 j_novak
59 * Modified #include directives to use c++ syntax.
60 *
61 * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
62 * LORENE
63 *
64 * Revision 2.0 2000/03/06 10:27:53 eric
65 * *** empty log message ***
66 *
67 *
68 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/som_symy.C,v 1.5 2017/02/24 15:34:59 j_novak Exp $
69 *
70 */
71
72
73// Headers C
74#include <cmath>
75
76namespace Lorene {
77
78//****************************************************************************
79// Sommation en r
80//****************************************************************************
81
82
84 // Cas R_CHEB //
86
87void som_r_cheb_symy
88 (double* ti, const int nr, const int nt, const int np, const double x,
89 double* trtp) {
90
91// Variables diverses
92int i, j, k ;
93double* pi = ti ; // pointeur courant sur l'entree
94double* po = trtp ; // pointeur courant sur la sortie
95
96 // Valeurs des polynomes de Chebyshev au point x demande
97 double* cheb = new double [nr] ;
98 cheb[0] = 1. ;
99 cheb[1] = x ;
100 for (i=2; i<nr; i++) {
101 cheb[i] = 2*x* cheb[i-1] - cheb[i-2] ;
102 }
103
104 // Sommation pour le premier phi, k=0
105 for (j=0 ; j<nt ; j++) {
106 *po = 0 ;
107 // Sommation sur r
108 for (i=0 ; i<nr ; i++) {
109 *po += (*pi) * cheb[i] ;
110 pi++ ; // R suivant
111 }
112 po++ ; // Theta suivant
113 }
114
115 if (np > 1) {
116
117 // On saute le deuxieme phi (sin(0)), k=1
118 pi += nr*nt ;
119 po += nt ;
120
121 // Sommation sur les phi suivants (pour k=2,...,np)
122 // en sautant les sinus (d'ou le k+=2)
123 for (k=2 ; k<np+1 ; k+=2) {
124 for (j=0 ; j<nt ; j++) {
125 // Sommation sur r
126 *po = 0 ;
127 for (i=0 ; i<nr ; i++) {
128 *po += (*pi) * cheb[i] ;
129 pi++ ; // R suivant
130 }
131 po++ ; // Theta suivant
132 }
133 // On saute le sin(k*phi) :
134 pi += nr*nt ;
135 po += nt ;
136 }
137
138 } // fin du cas np > 1
139
140 // Menage
141 delete [] cheb ;
142}
143
144
145
147 // Cas R_CHEBU //
149
150void som_r_chebu_symy
151 (double* ti, const int nr, const int nt, const int np, const double x, double* trtp) {
152
153// Variables diverses
154int i, j, k ;
155double* pi = ti ; // pointeur courant sur l'entree
156double* po = trtp ; // pointeur courant sur la sortie
157
158 // Valeurs des polynomes de Chebyshev au point x demande
159 double* cheb = new double [nr] ;
160 cheb[0] = 1. ;
161 cheb[1] = x ;
162 for (i=2; i<nr; i++) {
163 cheb[i] = 2*x* cheb[i-1] - cheb[i-2] ;
164 }
165
166 // Sommation pour le premier phi, k=0
167 for (j=0 ; j<nt ; j++) {
168 *po = 0 ;
169 // Sommation sur r
170 for (i=0 ; i<nr ; i++) {
171 *po += (*pi) * cheb[i] ;
172 pi++ ; // R suivant
173 }
174 po++ ; // Theta suivant
175 }
176
177 if (np > 1) {
178
179 // On saute le deuxieme phi (sin(0)), k=1
180 pi += nr*nt ;
181 po += nt ;
182
183 // Sommation sur les phi suivants (pour k=2,...,np)
184 // en sautant les sinus (d'ou le k+=2)
185 for (k=2 ; k<np+1 ; k+=2) {
186 for (j=0 ; j<nt ; j++) {
187 // Sommation sur r
188 *po = 0 ;
189 for (i=0 ; i<nr ; i++) {
190 *po += (*pi) * cheb[i] ;
191 pi++ ; // R suivant
192 }
193 po++ ; // Theta suivant
194 }
195 // On saute le sin(k*phi) :
196 pi += nr*nt ;
197 po += nt ;
198 }
199
200 } // fin du cas np > 1
201
202 // Menage
203 delete [] cheb ;
204
205}
206
208 // Cas R_CHEBPIM_P //
210
211void som_r_chebpim_p_symy
212 (double* ti, const int nr, const int nt, const int np, const double x, double* trtp) {
213
214// Variables diverses
215int i, j, k ;
216double* pi = ti ; // pointeur courant sur l'entree
217double* po = trtp ; // pointeur courant sur la sortie
218
219 // Valeurs des polynomes de Chebyshev au point x demande
220 double* cheb = new double [2*nr] ;
221 cheb[0] = 1. ;
222 cheb[1] = x ;
223 for (i=2 ; i<2*nr ; i++) {
224 cheb[i] = 2*x* cheb[i-1] - cheb[i-2] ;
225 }
226
227 // Sommation pour le premier phi, k=0
228 int m = 0;
229 for (j=0 ; j<nt ; j++) {
230 *po = 0 ;
231 // Sommation sur r
232 for (i=0 ; i<nr ; i++) {
233 *po += (*pi) * cheb[2*i] ;
234 pi++ ; // R suivant
235 }
236 po++ ; // Theta suivant
237 }
238
239 if (np > 1) {
240
241 // On saute le deuxieme phi (sin(0)), k=1
242 pi += nr*nt ;
243 po += nt ;
244
245 // Sommation sur les phi suivants (pour k=2,...,np)
246 // en sautant les sinus (d'ou le k+=2)
247 for (k=2 ; k<np+1 ; k+=2) {
248 m = (k/2) % 2 ; // parite: 0 <-> 1
249 for (j=0 ; j<nt ; j++) {
250 // Sommation sur r
251 *po = 0 ;
252 for (i=0 ; i<nr ; i++) {
253 *po += (*pi) * cheb[2*i + m] ;
254 pi++ ; // R suivant
255 }
256 po++ ; // Theta suivant
257 }
258 // On saute le sin(k*phi) :
259 pi += nr*nt ;
260 po += nt ;
261 }
262
263 } // fin du cas np > 1
264
265 // Menage
266 delete [] cheb ;
267
268}
269
271 // Cas R_CHEBPIM_I //
273
274void som_r_chebpim_i_symy
275 (double* ti, const int nr, const int nt, const int np, const double x, double* trtp) {
276
277// Variables diverses
278int i, j, k ;
279double* pi = ti ; // pointeur courant sur l'entree
280double* po = trtp ; // pointeur courant sur la sortie
281
282 // Valeurs des polynomes de Chebyshev au point x demande
283 double* cheb = new double [2*nr] ;
284 cheb[0] = 1. ;
285 cheb[1] = x ;
286 for (i=2 ; i<2*nr ; i++) {
287 cheb[i] = 2*x* cheb[i-1] - cheb[i-2] ;
288 }
289
290 // Sommation pour le premier phi, k=0
291 int m = 0;
292 for (j=0 ; j<nt ; j++) {
293 *po = 0 ;
294 // Sommation sur r
295 for (i=0 ; i<nr ; i++) {
296 *po += (*pi) * cheb[2*i+1] ;
297 pi++ ; // R suivant
298 }
299 po++ ; // Theta suivant
300 }
301
302 if (np > 1) {
303
304 // On saute le deuxieme phi (sin(0)), k=1
305 pi += nr*nt ;
306 po += nt ;
307
308 // Sommation sur les phi suivants (pour k=2,...,np)
309 // en sautant les sinus (d'ou le k+=2)
310 for (k=2 ; k<np+1 ; k+=2) {
311 m = (k/2) % 2 ; // parite: 0 <-> 1
312 for (j=0 ; j<nt ; j++) {
313 // Sommation sur r
314 *po = 0 ;
315 for (i=0 ; i<nr ; i++) {
316 *po += (*pi) * cheb[2*i + 1 - m] ;
317 pi++ ; // R suivant
318 }
319 po++ ; // Theta suivant
320 }
321 // On saute le sin(k*phi) :
322 pi += nr*nt ;
323 po += nt ;
324 }
325
326 } // fin du cas np > 1
327
328 // Menage
329 delete [] cheb ;
330
331}
332
333
334
335//****************************************************************************
336// Sommation en theta
337//****************************************************************************
338
340 // Cas T_COSSIN_CP //
342
343void som_tet_cossin_cp_symy
344 (double* ti, const int nt, const int np,
345 const double tet, double* to) {
346
347// Variables diverses
348int j, k ;
349double* pi = ti ; // Pointeur courant sur l'entree
350double* po = to ; // Pointeur courant sur la sortie
351
352 // Initialisation des tables trigo
353 double* cossin = new double [2*nt] ;
354 for (j=0 ; j<2*nt ; j +=2) {
355 cossin[j] = cos(j * tet) ;
356 cossin[j+1] = sin((j+1) * tet) ;
357 }
358
359 // Sommation sur le premier phi -> cosinus, k=0
360 *po = 0 ;
361 for (j=0 ; j<nt ; j++) {
362 *po += (*pi) * cossin[2*j] ;
363 pi++ ; // Theta suivant
364 }
365 po++ ; // Phi suivant
366
367 if (np > 1) {
368
369 // On saute le phi suivant (sin(0)), k=1
370 pi += nt ;
371 po++ ;
372
373 // Sommation sur le reste des phi (pour k=2,...,np), suivant parite de m
374 for (k=2 ; k<np+1 ; k+=2) {
375 int m = (k/2) % 2 ; // parite: 0 <-> 1
376 (*po) = 0 ;
377 for (j=0 ; j<nt ; j++) {
378 *po += (*pi) * cossin[2*j + m] ;
379 pi++ ; // Theta suivant
380 }
381 po++ ; // Phi suivant
382
383 // On saute le sin(k*phi) :
384 pi += nt ;
385 po++ ;
386
387 }
388 } // fin du cas np > 1
389
390 // Menage
391 delete [] cossin ;
392}
393
395 // Cas T_COSSIN_CI //
397
398void som_tet_cossin_ci_symy
399 (double* ti, const int nt, const int np,
400 const double tet, double* to) {
401
402// Variables diverses
403int j, k ;
404double* pi = ti ; // Pointeur courant sur l'entree
405double* po = to ; // Pointeur courant sur la sortie
406
407 // Initialisation des tables trigo
408 double* cossin = new double [2*nt] ;
409 for (j=0 ; j<2*nt ; j +=2) {
410 cossin[j] = cos((j+1) * tet) ;
411 cossin[j+1] = sin(j * tet) ;
412 }
413
414 // Sommation sur le premier phi -> cosinus, k=0
415 *po = 0 ;
416 for (j=0 ; j<nt ; j++) {
417 *po += (*pi) * cossin[2*j] ;
418 pi++ ; // Theta suivant
419 }
420 po++ ; // Phi suivant
421
422 if (np > 1) {
423
424 // On saute le phi suivant (sin(0)), k=1
425 pi += nt ;
426 po++ ;
427
428 // Sommation sur le reste des phi (pour k=2,...,np), suivant parite de m
429 for (k=2 ; k<np+1 ; k+=2) {
430 int m = (k/2) % 2 ; // parite: 0 <-> 1
431 (*po) = 0 ;
432 for (j=0 ; j<nt ; j++) {
433 *po += (*pi) * cossin[2*j + m] ;
434 pi++ ; // Theta suivant
435 }
436 po++ ; // Phi suivant
437
438 // On saute le sin(k*phi) :
439 pi += nt ;
440 po++ ;
441
442 }
443 } // fin du cas np > 1
444
445 // Menage
446 delete [] cossin ;
447}
448
449
450//****************************************************************************
451// Sommation en phi
452//****************************************************************************
453
454void som_phi_cossin_symy
455 (double* ti, const int np, const double phi, double* xo) {
456
457 *xo = ti[0] ; // premier element
458
459 // Sommation sur les cosinus seulement
460 for (int k=2 ; k<np-1 ; k +=2 ) {
461 int m = k/2 ;
462 *xo += ti[k] * cos(m * phi) ;
463 }
464 *xo += ti[np] * cos(np/2 * phi) ;
465}
466
467}
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::phi
Coord phi
coordinate centered on the grid
Definition map.h:732
Lorene::tet
Coord tet
coordinate centered on the grid
Definition map.h:731
Lorene::x
Coord x
x coordinate centered on the grid
Definition map.h:738