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