LORENE
FFT991/cfrchebpimi.C
1 /*
2  * Copyright (c) 1999-2002 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 cfrchebpimi_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/FFT991/cfrchebpimi.C,v 1.4 2014/10/15 12:48:20 j_novak Exp $" ;
24 
25 /*
26  * Transformation de Tchebyshev T_{2k+1}/T_{2k} sur le troisieme indice (indice
27  * correspondant a r) d'un tableau 3-D decrivant par exemple d/dr d'une
28  * fonction symetrique par rapport au plan equatorial z = 0 et sans aucune
29  * autre symetrie, cad que l'on effectue
30  * 1/ un developpement en polynomes de Tchebyshev impairs pour m pair
31  * 2/ un developpement en polynomes de Tchebyshev pairs pour m impair
32  *
33  * Utilise la routine FFT Fortran FFT991
34  *
35  *
36  * Entree:
37  * -------
38  * int* deg : tableau du nombre effectif de degres de liberte dans chacune
39  * des 3 dimensions: le nombre de points de collocation
40  * en r est nr = deg[2] et doit etre de la forme
41  * nr = 2^p 3^q 5^r + 1
42  * int* dimf : tableau du nombre d'elements de ff dans chacune des trois
43  * dimensions.
44  * On doit avoir dimf[2] >= deg[2] = nr.
45  *
46  * double* ff : tableau des valeurs de la fonction aux nr points de
47  * de collocation
48  *
49  * x_i = sin( pi/2 i/(nr-1) ) 0 <= i <= nr-1
50  *
51  * Les valeurs de la fonction doivent etre stokees dans le
52  * tableau ff comme suit
53  * f( x_i ) = ff[ dimf[1]*dimf[2] * j + dimf[2] * k + i ]
54  * ou j et k sont les indices correspondant a phi et theta
55  * respectivement.
56  * L'espace memoire correspondant a ce pointeur doit etre
57  * dimf[0]*dimf[1]*dimf[2] et doit etre alloue avant l'appel a
58  * la routine.
59  *
60  * int* dimc : tableau du nombre d'elements de cf dans chacune des trois
61  * dimensions.
62  * On doit avoir dimc[2] >= deg[2] = nr.
63  *
64  * Sortie:
65  * -------
66  * double* cf : tableau des nr coefficients c_i de la fonction definis
67  * comme suit (a theta et phi fixes)
68  *
69  * -- pour m pair (i.e. j = 0, 1, 4, 5, 8, 9, ...) :
70  *
71  * f(x) = som_{i=0}^{nr-2} c_i T_{2i+1}(x)
72  *
73  * ou T_{2i+1}(x) designe le polynome de Tchebyshev de
74  * degre 2i+1.
75  *
76  * -- pour m impair (i.e. j = 2, 3, 6, 7, 10, 11, ...) :
77  *
78  * f(x) = som_{i=0}^{nr-1} c_i T_{2i}(x) ,
79  *
80  * ou T_{2i}(x) designe le polynome de Tchebyshev de
81  * degre 2i.
82  *
83  * Les coefficients c_i (0 <= i <= nr-1) sont stokes dans
84  * le tableau cf comme suit
85  * c_i = cf[ dimc[1]*dimc[2] * j + dimc[2] * k + i ]
86  * ou j et k sont les indices correspondant a phi et theta
87  * respectivement.
88  * L'espace memoire correspondant a ce pointeur doit etre
89  * dimc[0]*dimc[1]*dimc[2] et doit avoir ete alloue avant
90  * l'appel a la routine.
91  *
92  * NB: Si le pointeur ff est egal a cf, la routine ne travaille que sur un
93  * seul tableau, qui constitue une entree/sortie.
94  */
95 
96 /*
97  * $Id: cfrchebpimi.C,v 1.4 2014/10/15 12:48:20 j_novak Exp $
98  * $Log: cfrchebpimi.C,v $
99  * Revision 1.4 2014/10/15 12:48:20 j_novak
100  * Corrected namespace declaration.
101  *
102  * Revision 1.3 2014/10/13 08:53:15 j_novak
103  * Lorene classes and functions now belong to the namespace Lorene.
104  *
105  * Revision 1.2 2014/10/06 15:18:44 j_novak
106  * Modified #include directives to use c++ syntax.
107  *
108  * Revision 1.1 2004/12/21 17:06:01 j_novak
109  * Added all files for using fftw3.
110  *
111  * Revision 1.4 2003/01/31 10:31:23 e_gourgoulhon
112  * Suppressed the directive #include <malloc.h> for malloc is defined
113  * in <stdlib.h>
114  *
115  * Revision 1.3 2002/10/16 14:36:44 j_novak
116  * Reorganization of #include instructions of standard C++, in order to
117  * use experimental version 3 of gcc.
118  *
119  * Revision 1.2 2002/09/09 13:00:39 e_gourgoulhon
120  * Modification of declaration of Fortran 77 prototypes for
121  * a better portability (in particular on IBM AIX systems):
122  * All Fortran subroutine names are now written F77_* and are
123  * defined in the new file C++/Include/proto_f77.h.
124  *
125  * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
126  * LORENE
127  *
128  * Revision 2.0 1999/02/22 15:48:21 hyc
129  * *** empty log message ***
130  *
131  *
132  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/FFT991/cfrchebpimi.C,v 1.4 2014/10/15 12:48:20 j_novak Exp $
133  *
134  */
135 
136 // headers du C
137 #include <cassert>
138 #include <cstdlib>
139 
140 // Prototypes of F77 subroutines
141 #include "headcpp.h"
142 #include "proto_f77.h"
143 
144 // Prototypage des sous-routines utilisees:
145 namespace Lorene {
146 int* facto_ini(int ) ;
147 double* trigo_ini(int ) ;
148 double* cheb_ini(const int) ;
149 double* chebimp_ini(const int ) ;
150 
151 //*****************************************************************************
152 
153 void cfrchebpimi(const int* deg, const int* dimf, double* ff, const int* dimc,
154  double* cf)
155 
156 {
157 
158 int i, j, k ;
159 
160 // Dimensions des tableaux ff et cf :
161  int n1f = dimf[0] ;
162  int n2f = dimf[1] ;
163  int n3f = dimf[2] ;
164  int n1c = dimc[0] ;
165  int n2c = dimc[1] ;
166  int n3c = dimc[2] ;
167 
168 // Nombres de degres de liberte en r :
169  int nr = deg[2] ;
170 
171 // Tests de dimension:
172  if (nr > n3f) {
173  cout << "cfrchebpimi: nr > n3f : nr = " << nr << " , n3f = "
174  << n3f << endl ;
175  abort () ;
176  exit(-1) ;
177  }
178  if (nr > n3c) {
179  cout << "cfrchebpimi: nr > n3c : nr = " << nr << " , n3c = "
180  << n3c << endl ;
181  abort () ;
182  exit(-1) ;
183  }
184  if (n1f > n1c) {
185  cout << "cfrchebpimi: n1f > n1c : n1f = " << n1f << " , n1c = "
186  << n1c << endl ;
187  abort () ;
188  exit(-1) ;
189  }
190  if (n2f > n2c) {
191  cout << "cfrchebpimi: n2f > n2c : n2f = " << n2f << " , n2c = "
192  << n2c << endl ;
193  abort () ;
194  exit(-1) ;
195  }
196 
197 // Nombre de points pour la FFT:
198  int nm1 = nr - 1;
199  int nm1s2 = nm1 / 2;
200 
201 // Recherche des tables pour la FFT:
202  int* facto = facto_ini(nm1) ;
203  double* trigo = trigo_ini(nm1) ;
204 
205 // Recherche de la table des sin(psi) :
206  double* sinp = cheb_ini(nr);
207 
208 // Recherche de la table des points de collocations x_k :
209  double* x = chebimp_ini(nr);
210 
211  // tableau de travail G et t1
212  // (la dimension nm1+2 = nr+1 est exigee par la routine fft991)
213  double* g = (double*)( malloc( (nm1+2)*sizeof(double) ) );
214  double* t1 = (double*)( malloc( (nm1+2)*sizeof(double) ) ) ;
215 
216 // Parametres pour la routine FFT991
217  int jump = 1 ;
218  int inc = 1 ;
219  int lot = 1 ;
220  int isign = - 1 ;
221 
222 // boucle sur phi et theta
223 
224  int n2n3f = n2f * n3f ;
225  int n2n3c = n2c * n3c ;
226 
227 //=======================================================================
228 // Cas m pair
229 //=======================================================================
230 
231  j = 0 ;
232 
233  while (j<n1f-1) {
234 
235 //------------------------------------------------------------------------
236 // partie cos(m phi) avec m pair : developpement en T_{2i+1}(x)
237 //------------------------------------------------------------------------
238 
239  for (k=0; k<n2f; k++) {
240 
241  int i0 = n2n3f * j + n3f * k ; // indice de depart
242  double* ff0 = ff + i0 ; // tableau des donnees a transformer
243 
244  i0 = n2n3c * j + n3c * k ; // indice de depart
245  double* cf0 = cf + i0 ; // tableau resultat
246 
247 // Multiplication de la fonction par x (pour la rendre paire)
248 // NB: dans les commentaires qui suivent, on note h(x) = x f(x).
249 // (Les valeurs de h dans l'espace des configurations sont stokees dans le
250 // tableau cf0).
251  cf0[0] = 0 ;
252  for (i=1; i<nr; i++) cf0[i] = x[i] * ff0[i] ;
253 
254 /*
255  * NB: dans les commentaires qui suivent, psi designe la variable de [0, pi]
256  * reliee a x par x = cos(psi/2) et F(psi) = h(x(psi)).
257  */
258 
259 // Valeur en psi=0 de la partie antisymetrique de F, F_ :
260  double fmoins0 = 0.5 * ( cf0[nm1] - cf0[0] );
261 
262 // Fonction G(psi) = F+(psi) + F_(psi) sin(psi)
263 //---------------------------------------------
264  for ( i = 1; i < nm1s2 ; i++ ) {
265 // ... indice (dans le tableau g) du pt symetrique de psi par rapport a pi/2:
266  int isym = nm1 - i ;
267 // ... indice (dans le tableau cf0) du point x correspondant a psi
268  int ix = nm1 - i ;
269 // ... indice (dans le tableau cf0) du point x correspondant a sym(psi)
270  int ixsym = nm1 - isym ;
271 
272 // ... F+(psi)
273  double fp = 0.5 * ( cf0[ix] + cf0[ixsym] ) ;
274 // ... F_(psi) sin(psi)
275  double fms = 0.5 * ( cf0[ix] - cf0[ixsym] ) * sinp[i] ;
276  g[i] = fp + fms ;
277  g[isym] = fp - fms ;
278  }
279 //... cas particuliers:
280  g[0] = 0.5 * ( cf0[nm1] + cf0[0] );
281  g[nm1s2] = cf0[nm1s2];
282 
283 // Developpement de G en series de Fourier par une FFT
284 //----------------------------------------------------
285 
286  F77_fft991( g, t1, trigo, facto, &inc, &jump, &nm1, &lot, &isign) ;
287 
288 // Coefficients pairs du developmt. de Tchebyshev de h
289 //----------------------------------------------------
290 // Ces coefficients sont egaux aux coefficients en cosinus du developpement
291 // de G en series de Fourier (le facteur 2. vient de la normalisation
292 // de fft991; si fft991 donnait reellement les coef. en cosinus, il faudrait le
293 // remplacer par un +1.) :
294 
295  cf0[0] = g[0] ;
296  for (i=2; i<nm1; i += 2 ) cf0[i] = 2. * g[i] ;
297  cf0[nm1] = g[nm1] ;
298 
299 // Coefficients impairs du developmt. de Tchebyshev de h
300 //------------------------------------------------------
301 // 1. Coef. c'_k (recurrence amorcee a partir de zero)
302 // Le +4. en facteur de g[i] est du a la normalisation de fft991
303 // (si fft991 donnait reellement les coef. en sinus, il faudrait le
304 // remplacer par un -2.)
305  cf0[1] = 0 ;
306  double som = 0;
307  for ( i = 3; i < nr; i += 2 ) {
308  cf0[i] = cf0[i-2] + 4. * g[i] ;
309  som += cf0[i] ;
310  }
311 
312 // 2. Calcul de c_1 :
313  double c1 = ( fmoins0 - som ) / nm1s2 ;
314 
315 // 3. Coef. c_k avec k impair:
316  cf0[1] = c1 ;
317  for ( i = 3; i < nr; i += 2 ) cf0[i] += c1 ;
318 
319 // Coefficients de f en fonction de ceux de h
320 //-------------------------------------------
321 
322  cf0[0] = 2* cf0[0] ;
323  for (i=1; i<nm1; i++) {
324  cf0[i] = 2 * cf0[i] - cf0[i-1] ;
325  }
326  cf0[nm1] = 0 ;
327 
328 
329  } // fin de la boucle sur theta
330 
331 
332 //------------------------------------------------------------------------
333 // partie sin(m phi) avec m pair : developpement en T_{2i+1}(x)
334 //------------------------------------------------------------------------
335 
336  j++ ;
337 
338  if ( (j != 1) && (j != n1f-1) ) {
339 // on effectue le calcul seulement dans les cas ou les coef en phi ne sont
340 // pas nuls
341 
342  for (k=0; k<n2f; k++) {
343 
344  int i0 = n2n3f * j + n3f * k ; // indice de depart
345  double* ff0 = ff + i0 ; // tableau des donnees a transformer
346 
347  i0 = n2n3c * j + n3c * k ; // indice de depart
348  double* cf0 = cf + i0 ; // tableau resultat
349 
350 // Multiplication de la fonction par x (pour la rendre paire)
351 // NB: dans les commentaires qui suivent, on note h(x) = x f(x).
352 // (Les valeurs de h dans l'espace des configurations sont stokees dans le
353 // tableau cf0).
354  cf0[0] = 0 ;
355  for (i=1; i<nr; i++) cf0[i] = x[i] * ff0[i] ;
356 
357 /*
358  * NB: dans les commentaires qui suivent, psi designe la variable de [0, pi]
359  * reliee a x par x = cos(psi/2) et F(psi) = h(x(psi)).
360  */
361 
362 // Valeur en psi=0 de la partie antisymetrique de F, F_ :
363  double fmoins0 = 0.5 * ( cf0[nm1] - cf0[0] );
364 
365 // Fonction G(psi) = F+(psi) + F_(psi) sin(psi)
366 //---------------------------------------------
367  for ( i = 1; i < nm1s2 ; i++ ) {
368 // ... indice (dans le tableau g) du pt symetrique de psi par rapport a pi/2:
369  int isym = nm1 - i ;
370 // ... indice (dans le tableau cf0) du point x correspondant a psi
371  int ix = nm1 - i ;
372 // ... indice (dans le tableau cf0) du point x correspondant a sym(psi)
373  int ixsym = nm1 - isym ;
374 
375 // ... F+(psi)
376  double fp = 0.5 * ( cf0[ix] + cf0[ixsym] ) ;
377 // ... F_(psi) sin(psi)
378  double fms = 0.5 * ( cf0[ix] - cf0[ixsym] ) * sinp[i] ;
379  g[i] = fp + fms ;
380  g[isym] = fp - fms ;
381  }
382 //... cas particuliers:
383  g[0] = 0.5 * ( cf0[nm1] + cf0[0] );
384  g[nm1s2] = cf0[nm1s2];
385 
386 // Developpement de G en series de Fourier par une FFT
387 //----------------------------------------------------
388 
389  F77_fft991( g, t1, trigo, facto, &inc, &jump, &nm1, &lot, &isign) ;
390 
391 // Coefficients pairs du developmt. de Tchebyshev de h
392 //----------------------------------------------------
393 // Ces coefficients sont egaux aux coefficients en cosinus du developpement
394 // de G en series de Fourier (le facteur 2. vient de la normalisation
395 // de fft991; si fft991 donnait reellement les coef. en cosinus, il faudrait le
396 // remplacer par un +1.) :
397 
398  cf0[0] = g[0] ;
399  for (i=2; i<nm1; i += 2 ) cf0[i] = 2. * g[i] ;
400  cf0[nm1] = g[nm1] ;
401 
402 // Coefficients impairs du developmt. de Tchebyshev de h
403 //------------------------------------------------------
404 // 1. Coef. c'_k (recurrence amorcee a partir de zero)
405 // Le +4. en facteur de g[i] est du a la normalisation de fft991
406 // (si fft991 donnait reellement les coef. en sinus, il faudrait le
407 // remplacer par un -2.)
408  cf0[1] = 0 ;
409  double som = 0;
410  for ( i = 3; i < nr; i += 2 ) {
411  cf0[i] = cf0[i-2] + 4. * g[i] ;
412  som += cf0[i] ;
413  }
414 
415 // 2. Calcul de c_1 :
416  double c1 = ( fmoins0 - som ) / nm1s2 ;
417 
418 // 3. Coef. c_k avec k impair:
419  cf0[1] = c1 ;
420  for ( i = 3; i < nr; i += 2 ) cf0[i] += c1 ;
421 
422 // Coefficients de f en fonction de ceux de h
423 //-------------------------------------------
424 
425  cf0[0] = 2* cf0[0] ;
426  for (i=1; i<nm1; i++) {
427  cf0[i] = 2 * cf0[i] - cf0[i-1] ;
428  }
429  cf0[nm1] = 0 ;
430 
431  } // fin de la boucle sur theta
432 
433  } // fin du cas ou le calcul etait necessaire (i.e. ou les
434  // coef en phi n'etaient pas nuls)
435 
436 // On passe au cas m pair suivant:
437  j+=3 ;
438 
439  } // fin de la boucle sur les cas m pair
440 
441  if (n1f<=3) { // cas m=0 seulement (symetrie axiale)
442  free (t1) ;
443  free (g) ;
444  return ;
445  }
446 
447 //=======================================================================
448 // Cas m impair
449 //=======================================================================
450 
451  j = 2 ;
452 
453  while (j<n1f-1) {
454 
455 //--------------------------------------------------------------------
456 // partie cos(m phi) avec m impair : developpement en T_{2i}(x)
457 //--------------------------------------------------------------------
458 
459  for (k=0; k<n2f; k++) {
460 
461  int i0 = n2n3f * j + n3f * k ; // indice de depart
462  double* ff0 = ff + i0 ; // tableau des donnees a transformer
463 
464  i0 = n2n3c * j + n3c * k ; // indice de depart
465  double* cf0 = cf + i0 ; // tableau resultat
466 
467 /*
468  * NB: dans les commentaires qui suivent, psi designe la variable de [0, pi]
469  * reliee a x par x = cos(psi/2) et F(psi) = f(x(psi)).
470  */
471 
472 // Valeur en psi=0 de la partie antisymetrique de F, F_ :
473  double fmoins0 = 0.5 * ( ff0[nm1] - ff0[0] );
474 
475 // Fonction G(psi) = F+(psi) + F_(psi) sin(psi)
476 //---------------------------------------------
477  for ( i = 1; i < nm1s2 ; i++ ) {
478 // ... indice (dans le tableau g) du pt symetrique de psi par rapport a pi/2:
479  int isym = nm1 - i ;
480 // ... indice (dans le tableau ff0) du point x correspondant a psi
481  int ix = nm1 - i ;
482 // ... indice (dans le tableau ff0) du point x correspondant a sym(psi)
483  int ixsym = nm1 - isym ;
484 
485 // ... F+(psi)
486  double fp = 0.5 * ( ff0[ix] + ff0[ixsym] ) ;
487 // ... F_(psi) sin(psi)
488  double fms = 0.5 * ( ff0[ix] - ff0[ixsym] ) * sinp[i] ;
489  g[i] = fp + fms ;
490  g[isym] = fp - fms ;
491  }
492 //... cas particuliers:
493  g[0] = 0.5 * ( ff0[nm1] + ff0[0] );
494  g[nm1s2] = ff0[nm1s2];
495 
496 // Developpement de G en series de Fourier par une FFT
497 //----------------------------------------------------
498 
499  F77_fft991( g, t1, trigo, facto, &inc, &jump, &nm1, &lot, &isign) ;
500 
501 // Coefficients pairs du developmt. de Tchebyshev de f
502 //----------------------------------------------------
503 // Ces coefficients sont egaux aux coefficients en cosinus du developpement
504 // de G en series de Fourier (le facteur 2 vient de la normalisation
505 // de fft991) :
506 
507  cf0[0] = g[0] ;
508  for (i=2; i<nm1; i += 2 ) cf0[i] = 2. * g[i] ;
509  cf0[nm1] = g[nm1] ;
510 
511 // Coefficients impairs du developmt. de Tchebyshev de f
512 //------------------------------------------------------
513 // 1. Coef. c'_k (recurrence amorcee a partir de zero)
514 // Le +4. en facteur de g[i] est du a la normalisation de fft991
515 // (si fft991 donnait reellement les coef. en sinus, il faudrait le
516 // remplacer par un -2.)
517  cf0[1] = 0 ;
518  double som = 0;
519  for ( i = 3; i < nr; i += 2 ) {
520  cf0[i] = cf0[i-2] + 4. * g[i] ;
521  som += cf0[i] ;
522  }
523 
524 // 2. Calcul de c_1 :
525  double c1 = ( fmoins0 - som ) / nm1s2 ;
526 
527 // 3. Coef. c_k avec k impair:
528  cf0[1] = c1 ;
529  for ( i = 3; i < nr; i += 2 ) cf0[i] += c1 ;
530 
531  } // fin de la boucle sur theta
532 
533 
534 //--------------------------------------------------------------------
535 // partie sin(m phi) avec m impair : developpement en T_{2i}(x)
536 //--------------------------------------------------------------------
537 
538  j++ ;
539 
540  if ( j != n1f-1 ) {
541 // on effectue le calcul seulement dans les cas ou les coef en phi ne sont
542 // pas nuls
543 
544  for (k=0; k<n2f; k++) {
545 
546  int i0 = n2n3f * j + n3f * k ; // indice de depart
547  double* ff0 = ff + i0 ; // tableau des donnees a transformer
548 
549  i0 = n2n3c * j + n3c * k ; // indice de depart
550  double* cf0 = cf + i0 ; // tableau resultat
551 
552 /*
553  * NB: dans les commentaires qui suivent, psi designe la variable de [0, pi]
554  * reliee a x par x = cos(psi/2) et F(psi) = f(x(psi)).
555  */
556 
557 // Valeur en psi=0 de la partie antisymetrique de F, F_ :
558  double fmoins0 = 0.5 * ( ff0[nm1] - ff0[0] );
559 
560 // Fonction G(psi) = F+(psi) + F_(psi) sin(psi)
561 //---------------------------------------------
562  for ( i = 1; i < nm1s2 ; i++ ) {
563 // ... indice (dans le tableau g) du pt symetrique de psi par rapport a pi/2:
564  int isym = nm1 - i ;
565 // ... indice (dans le tableau ff0) du point x correspondant a psi
566  int ix = nm1 - i ;
567 // ... indice (dans le tableau ff0) du point x correspondant a sym(psi)
568  int ixsym = nm1 - isym ;
569 
570 // ... F+(psi)
571  double fp = 0.5 * ( ff0[ix] + ff0[ixsym] ) ;
572 // ... F_(psi) sin(psi)
573  double fms = 0.5 * ( ff0[ix] - ff0[ixsym] ) * sinp[i] ;
574  g[i] = fp + fms ;
575  g[isym] = fp - fms ;
576  }
577 //... cas particuliers:
578  g[0] = 0.5 * ( ff0[nm1] + ff0[0] );
579  g[nm1s2] = ff0[nm1s2];
580 
581 // Developpement de G en series de Fourier par une FFT
582 //----------------------------------------------------
583 
584  F77_fft991( g, t1, trigo, facto, &inc, &jump, &nm1, &lot, &isign) ;
585 
586 // Coefficients pairs du developmt. de Tchebyshev de f
587 //----------------------------------------------------
588 // Ces coefficients sont egaux aux coefficients en cosinus du developpement
589 // de G en series de Fourier (le facteur 2 vient de la normalisation
590 // de fft991) :
591 
592  cf0[0] = g[0] ;
593  for (i=2; i<nm1; i += 2 ) cf0[i] = 2. * g[i] ;
594  cf0[nm1] = g[nm1] ;
595 
596 // Coefficients impairs du developmt. de Tchebyshev de f
597 //------------------------------------------------------
598 // 1. Coef. c'_k (recurrence amorcee a partir de zero)
599 // Le +4. en facteur de g[i] est du a la normalisation de fft991
600 // (si fft991 donnait reellement les coef. en sinus, il faudrait le
601 // remplacer par un -2.)
602  cf0[1] = 0 ;
603  double som = 0;
604  for ( i = 3; i < nr; i += 2 ) {
605  cf0[i] = cf0[i-2] + 4. * g[i] ;
606  som += cf0[i] ;
607  }
608 
609 // 2. Calcul de c_1 :
610  double c1 = ( fmoins0 - som ) / nm1s2 ;
611 
612 // 3. Coef. c_k avec k impair:
613  cf0[1] = c1 ;
614  for ( i = 3; i < nr; i += 2 ) cf0[i] += c1 ;
615 
616  } // fin de la boucle sur theta
617 
618  } // fin du cas ou le calcul etait necessaire (i.e. ou les
619  // coef en phi n'etaient pas nuls)
620 
621 // On passe au cas m impair suivant:
622  j+=3 ;
623 
624  } // fin de la boucle sur les cas m impair
625 
626  // Menage
627  free (t1) ;
628  free (g) ;
629 }
630 }
Lorene
Lorene prototypes.
Definition: app_hor.h:64