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 char 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 
70 namespace 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
86 int i, j, k ;
87 double* pi = ti ; // pointeur courant sur l'entree
88 double* 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
148 int i, j, k ;
149 double* pi = ti ; // pointeur courant sur l'entree
150 double* 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
209 int i, j, k ;
210 double* pi = ti ; // pointeur courant sur l'entree
211 double* 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
272 int i, j, k ;
273 double* pi = ti ; // pointeur courant sur l'entree
274 double* 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
342 int j, k ;
343 double* pi = ti ; // Pointeur courant sur l'entree
344 double* 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
397 int j, k ;
398 double* pi = ti ; // Pointeur courant sur l'entree
399 double* 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 
448 void 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