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