LORENE
som_asymy.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_asymy_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/som_asymy.C,v 1.3 2014/10/13 08:53:26 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) + antisymetrie par rapport au plan y=0
28 *
29 * SYNOPSYS:
30 * double som_r_XX_asymy
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_asymy
37 * (double* ti, int nt, int np, double tet, double* to)
38 *
39 * double som_phi_XX_asymy
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_asymy.C,v 1.3 2014/10/13 08:53:26 j_novak Exp $
48 * $Log: som_asymy.C,v $
49 * Revision 1.3 2014/10/13 08:53:26 j_novak
50 * Lorene classes and functions now belong to the namespace Lorene.
51 *
52 * Revision 1.2 2014/10/06 15:16:06 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:45 eric
59 * *** empty log message ***
60 *
61 *
62 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Operators/som_asymy.C,v 1.3 2014/10/13 08:53:26 j_novak Exp $
63 *
64 */
65
66
67// Headers C
68#include <cmath>
69
70namespace Lorene {
71
72//****************************************************************************
73// Sommation en r
74//****************************************************************************
75
79
80void som_r_cheb_asymy(double* ti, const int nr, const int nt, const int np,
81 const double x, double* trtp) {
82
83// Variables diverses
84int i, j, k ;
85double* pi = ti ; // pointeur courant sur l'entree
86double* po = trtp ; // pointeur courant sur la sortie
87
88 // Valeurs des polynomes de Chebyshev au point x demande
89 double* cheb = new double [nr] ;
90 cheb[0] = 1. ;
91 cheb[1] = x ;
92 for (i=2; i<nr; i++) {
93 cheb[i] = 2*x* cheb[i-1] - cheb[i-2] ;
94 }
95
96 // On saute les 3 premiers coef. en phi, qui sont respectivemment
97 // cos(0 phi), sin(0 phi) et cos(phi)
98 pi += 3*nr*nt ;
99 po += 3*nt ;
100
101 // Sommation sur les phi suivants (pour k=3,...,np)
102 // en sautant les cosinus (d'ou le k+=2)
103 for (k=3 ; k<np+1 ; k+=2) {
104 for (j=0 ; j<nt ; j++) {
105 // Sommation sur r
106 *po = 0 ;
107 for (i=0 ; i<nr ; i++) {
108 *po += (*pi) * cheb[i] ;
109 pi++ ; // R suivant
110 }
111 po++ ; // Theta suivant
112 }
113 // On saute le cos(m*phi) :
114 pi += nr*nt ;
115 po += nt ;
116 }
117
118 // Menage
119 delete [] cheb ;
120}
121
122
123
127
128void som_r_chebu_asymy(double* ti, const int nr, const int nt, const int np,
129 const double x, double* trtp) {
130
131// Variables diverses
132int i, j, k ;
133double* pi = ti ; // pointeur courant sur l'entree
134double* po = trtp ; // pointeur courant sur la sortie
135
136 // Valeurs des polynomes de Chebyshev au point x demande
137 double* cheb = new double [nr] ;
138 cheb[0] = 1. ;
139 cheb[1] = x ;
140 for (i=2; i<nr; i++) {
141 cheb[i] = 2*x* cheb[i-1] - cheb[i-2] ;
142 }
143
144 // On saute les 3 premiers coef. en phi, qui sont respectivemment
145 // cos(0 phi), sin(0 phi) et cos(phi)
146 pi += 3*nr*nt ;
147 po += 3*nt ;
148
149 // Sommation sur les phi suivants (pour k=3,...,np)
150 // en sautant les cosinus (d'ou le k+=2)
151 for (k=3 ; k<np+1 ; k+=2) {
152 for (j=0 ; j<nt ; j++) {
153 // Sommation sur r
154 *po = 0 ;
155 for (i=0 ; i<nr ; i++) {
156 *po += (*pi) * cheb[i] ;
157 pi++ ; // R suivant
158 }
159 po++ ; // Theta suivant
160 }
161 // On saute le cos(m*phi) :
162 pi += nr*nt ;
163 po += nt ;
164 }
165
166 // Menage
167 delete [] cheb ;
168
169}
170
171
172
176
177void som_r_chebpim_p_asymy(double* ti, const int nr, const int nt,
178 const int np, const double x, double* trtp) {
179
180// Variables diverses
181int i, j, k ;
182double* pi = ti ; // pointeur courant sur l'entree
183double* po = trtp ; // pointeur courant sur la sortie
184
185 // Valeurs des polynomes de Chebyshev au point x demande
186 double* cheb = new double [2*nr] ;
187 cheb[0] = 1. ;
188 cheb[1] = x ;
189 for (i=2 ; i<2*nr ; i++) {
190 cheb[i] = 2*x* cheb[i-1] - cheb[i-2] ;
191 }
192
193
194 // On saute les 3 premiers coef. en phi, qui sont respectivemment
195 // cos(0 phi), sin(0 phi) et cos(phi)
196 pi += 3*nr*nt ;
197 po += 3*nt ;
198
199 // Sommation sur les phi suivants (pour k=3,...,np)
200 // en sautant les cosinus (d'ou le k+=2)
201 for (k=3 ; k<np+1 ; k+=2) {
202 int m = (k/2) % 2 ; // parite: 0 <-> 1
203 for (j=0 ; j<nt ; j++) {
204 // Sommation sur r
205 *po = 0 ;
206 for (i=0 ; i<nr ; i++) {
207 *po += (*pi) * cheb[2*i + m] ;
208 pi++ ; // R suivant
209 }
210 po++ ; // Theta suivant
211 }
212 // On saute le cos(m*phi) :
213 pi += nr*nt ;
214 po += nt ;
215 }
216
217
218 // Menage
219 delete [] cheb ;
220
221}
222
226
227void som_r_chebpim_i_asymy(double* ti, const int nr, const int nt,
228 const int np, const double x, double* trtp) {
229
230// Variables diverses
231int i, j, k ;
232double* pi = ti ; // pointeur courant sur l'entree
233double* po = trtp ; // pointeur courant sur la sortie
234
235 // Valeurs des polynomes de Chebyshev au point x demande
236 double* cheb = new double [2*nr] ;
237 cheb[0] = 1. ;
238 cheb[1] = x ;
239 for (i=2 ; i<2*nr ; i++) {
240 cheb[i] = 2*x* cheb[i-1] - cheb[i-2] ;
241 }
242
243 // On saute les 3 premiers coef. en phi, qui sont respectivemment
244 // cos(0 phi), sin(0 phi) et cos(phi)
245 pi += 3*nr*nt ;
246 po += 3*nt ;
247
248 // Sommation sur les phi suivants (pour k=3,...,np)
249 // en sautant les cosinus (d'ou le k+=2)
250 for (k=3 ; k<np+1 ; k+=2) {
251 int m = (k/2) % 2 ; // parite: 0 <-> 1
252 for (j=0 ; j<nt ; j++) {
253 // Sommation sur r
254 *po = 0 ;
255 for (i=0 ; i<nr ; i++) {
256 *po += (*pi) * cheb[2*i + 1 - m] ;
257 pi++ ; // R suivant
258 }
259 po++ ; // Theta suivant
260 }
261 // On saute le cos(m*phi) :
262 pi += nr*nt ;
263 po += nt ;
264 }
265
266 // Menage
267 delete [] cheb ;
268
269}
270
271
272
273//****************************************************************************
274// Sommation en theta
275//****************************************************************************
276
280
281void som_tet_cossin_cp_asymy(double* ti, const int nt, const int np,
282 const double tet, double* to) {
283
284// Variables diverses
285int j, k ;
286double* pi = ti ; // Pointeur courant sur l'entree
287double* po = to ; // Pointeur courant sur la sortie
288
289 // Initialisation des tables trigo
290 double* cossin = new double [2*nt] ;
291 for (j=0 ; j<2*nt ; j +=2) {
292 cossin[j] = cos(j * tet) ;
293 cossin[j+1] = sin((j+1) * tet) ;
294 }
295
296 // On saute les 3 premiers coef. en phi, qui sont respectivemment
297 // cos(0 phi), sin(0 phi) et cos(phi)
298 pi += 3*nt ;
299 po += 3 ;
300
301 // Sommation sur le reste des phi (pour k=3,...,np), suivant parite de m
302 for (k=3 ; k<np+1 ; k+=2) {
303 int m = (k/2) % 2 ; // parite: 0 <-> 1
304 (*po) = 0 ;
305 for (j=0 ; j<nt ; j++) {
306 *po += (*pi) * cossin[2*j + m] ;
307 pi++ ; // Theta suivant
308 }
309 po++ ; // Phi suivant
310
311 // On saute le cos(m*phi) :
312 pi += nt ;
313 po++ ;
314
315 }
316
317 // Menage
318 delete [] cossin ;
319}
320
324
325void som_tet_cossin_ci_asymy(double* ti, const int nt, const int np,
326 const double tet, double* to) {
327
328// Variables diverses
329int j, k ;
330double* pi = ti ; // Pointeur courant sur l'entree
331double* po = to ; // Pointeur courant sur la sortie
332
333 // Initialisation des tables trigo
334 double* cossin = new double [2*nt] ;
335 for (j=0 ; j<2*nt ; j +=2) {
336 cossin[j] = cos((j+1) * tet) ;
337 cossin[j+1] = sin(j * tet) ;
338 }
339
340 // On saute les 3 premiers coef. en phi, qui sont respectivemment
341 // cos(0 phi), sin(0 phi) et cos(phi)
342 pi += 3*nt ;
343 po += 3 ;
344
345 // Sommation sur le reste des phi (pour k=3,...,np), suivant parite de m
346 for (k=3 ; k<np+1 ; k+=2) {
347 int m = (k/2) % 2 ; // parite: 0 <-> 1
348 (*po) = 0 ;
349 for (j=0 ; j<nt ; j++) {
350 *po += (*pi) * cossin[2*j + m] ;
351 pi++ ; // Theta suivant
352 }
353 po++ ; // Phi suivant
354
355 // On saute le cos(m*phi) :
356 pi += nt ;
357 po++ ;
358
359 }
360
361 // Menage
362 delete [] cossin ;
363}
364
365
366//****************************************************************************
367// Sommation en phi
368//****************************************************************************
369
370void som_phi_cossin_asymy(double* ti, const int np, const double phi,
371 double* xo) {
372
373 *xo = 0 ;
374
375 // Sommation sur les sinus seulement
376 for (int k=3 ; k<np+1 ; k +=2 ) {
377 int m = k/2 ;
378 *xo += ti[k] * sin(m * phi) ;
379 }
380
381}
382
383}
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_asymy(double *, const int, const int, const int, const double, double *)
Cas R_CHEB ///.
Definition som_asymy.C:80
void som_tet_cossin_cp_asymy(double *, const int, const int, const double, double *)
Cas T_COSSIN_CP ///.
Definition som_asymy.C:281
void som_r_chebpim_p_asymy(double *, const int, const int, const int, const double, double *)
Cas R_CHEBPIM_P ///.
Definition som_asymy.C:177
void som_r_chebpim_i_asymy(double *, const int, const int, const int, const double, double *)
Cas R_CHEBPIM_I ///.
Definition som_asymy.C:227
void som_tet_cossin_ci_asymy(double *, const int, const int, const double, double *)
Cas T_COSSIN_CI ///.
Definition som_asymy.C:325
void som_r_chebu_asymy(double *, const int, const int, const int, const double, double *)
Cas R_CHEBU ///.
Definition som_asymy.C:128