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 
23 char 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 
70 namespace Lorene {
71 
72 //****************************************************************************
73 // Sommation en r
74 //****************************************************************************
75 
79 
80 void 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
84 int i, j, k ;
85 double* pi = ti ; // pointeur courant sur l'entree
86 double* 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 
128 void 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
132 int i, j, k ;
133 double* pi = ti ; // pointeur courant sur l'entree
134 double* 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 
177 void 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
181 int i, j, k ;
182 double* pi = ti ; // pointeur courant sur l'entree
183 double* 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 
227 void 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
231 int i, j, k ;
232 double* pi = ti ; // pointeur courant sur l'entree
233 double* 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 
281 void som_tet_cossin_cp_asymy(double* ti, const int nt, const int np,
282  const double tet, double* to) {
283 
284 // Variables diverses
285 int j, k ;
286 double* pi = ti ; // Pointeur courant sur l'entree
287 double* 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 
325 void som_tet_cossin_ci_asymy(double* ti, const int nt, const int np,
326  const double tet, double* to) {
327 
328 // Variables diverses
329 int j, k ;
330 double* pi = ti ; // Pointeur courant sur l'entree
331 double* 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 
370 void 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 }
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
Lorene prototypes.
Definition: app_hor.h:64
Cmp cos(const Cmp &)
Cosine.
Definition: cmp_math.C:94
void som_r_chebu_asymy(double *, const int, const int, const int, const double, double *)
Cas R_CHEBU ///.
Definition: som_asymy.C:128
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_chebpim_i_asymy(double *, const int, const int, const int, const double, double *)
Cas R_CHEBPIM_I ///.
Definition: som_asymy.C:227
void som_r_cheb_asymy(double *, const int, const int, const int, const double, double *)
Cas R_CHEB ///.
Definition: som_asymy.C:80
Cmp sin(const Cmp &)
Sine.
Definition: cmp_math.C:69