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