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