LORENE
FFT991/cfrchebpii.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 cfrchebpii_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/FFT991/cfrchebpii.C,v 1.3 2014/10/13 08:53:15 j_novak Exp $" ;
24
25
26/*
27 * Transformation de Tchebyshev (cas rare) sur le troisieme indice (indice
28 * correspondant a r) d'un tableau 3-D decrivant une fonction quelconque.
29 * Utilise la routine FFT Fortran FFT991
30 *
31 *
32 * Entree:
33 * -------
34 * int* deg : tableau du nombre effectif de degres de liberte dans chacune
35 * des 3 dimensions: le nombre de points de collocation
36 * en r est nr = deg[2] et doit etre de la forme
37 * nr = 2^p 3^q 5^r + 1
38 * int* dimf : tableau du nombre d'elements de ff dans chacune des trois
39 * dimensions.
40 * On doit avoir dimf[2] >= deg[2] = nr.
41 * NB: pour dimf[0] = 1 (un seul point en phi), la transformation
42 * est bien effectuee.
43 * pour dimf[0] > 1 (plus d'un point en phi), la
44 * transformation n'est effectuee que pour les indices (en phi)
45 * j != 1 et j != dimf[0]-1 (cf. commentaires sur borne_phi).
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 * Si l impair f(x) = som_{i=0}^{nr-1} c_i T_{2i}(x) ,
71 * Si l pair f(x) = som_{i=0}^{nr-1} c_i T_{2i+1}(x)
72 *
73 * ou T_{i}(x) designe le polynome de Tchebyshev de degre i.
74 * Les coefficients c_i (0 <= i <= nr-1) sont stokes dans
75 * le tableau cf comme suit
76 * c_i = cf[ dimc[1]*dimc[2] * j + dimc[2] * k + i ]
77 * ou j et k sont les indices correspondant a phi et theta
78 * respectivement.
79 * L'espace memoire correspondant a ce pointeur doit etre
80 * dimc[0]*dimc[1]*dimc[2] et doit avoir ete alloue avant
81 * l'appel a la routine.
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 * $Id: cfrchebpii.C,v 1.3 2014/10/13 08:53:15 j_novak Exp $
89 * $Log: cfrchebpii.C,v $
90 * Revision 1.3 2014/10/13 08:53:15 j_novak
91 * Lorene classes and functions now belong to the namespace Lorene.
92 *
93 * Revision 1.2 2014/10/06 15:18:44 j_novak
94 * Modified #include directives to use c++ syntax.
95 *
96 * Revision 1.1 2004/12/21 17:06:01 j_novak
97 * Added all files for using fftw3.
98 *
99 * Revision 1.1 2004/11/23 15:13:50 m_forot
100 * Added the bases for the cases without any equatorial symmetry
101 * (T_COSSIN_C, T_COSSIN_S, T_LEG, R_CHEBPI_P, R_CHEBPI_I).
102 *
103 *
104 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/FFT991/cfrchebpii.C,v 1.3 2014/10/13 08:53:15 j_novak Exp $
105 *
106 */
107
108
109// headers du C
110#include <cassert>
111#include <cstdlib>
112
113// Prototypes of F77 subroutines
114#include "headcpp.h"
115#include "proto_f77.h"
116
117// Prototypage des sous-routines utilisees:
118namespace Lorene {
119int* facto_ini(int ) ;
120double* trigo_ini(int ) ;
121double* cheb_ini(const int) ;
122double* chebimp_ini(const int ) ;
123
124//*****************************************************************************
125
126void cfrchebpii(const int* deg, const int* dimf, double* ff, const int* dimc,
127 double* cf)
128
129{
130
131int i, j, k ;
132
133// Dimensions des tableaux ff et cf :
134 int n1f = dimf[0] ;
135 int n2f = dimf[1] ;
136 int n3f = dimf[2] ;
137 int n1c = dimc[0] ;
138 int n2c = dimc[1] ;
139 int n3c = dimc[2] ;
140
141// Nombres de degres de liberte en r :
142 int nr = deg[2] ;
143
144// Tests de dimension:
145 if (nr > n3f) {
146 cout << "cfrchebpii: nr > n3f : nr = " << nr << " , n3f = "
147 << n3f << endl ;
148 abort () ;
149 exit(-1) ;
150 }
151 if (nr > n3c) {
152 cout << "cfrchebpii: nr > n3c : nr = " << nr << " , n3c = "
153 << n3c << endl ;
154 abort () ;
155 exit(-1) ;
156 }
157 if (n1f > n1c) {
158 cout << "cfrchebpii: n1f > n1c : n1f = " << n1f << " , n1c = "
159 << n1c << endl ;
160 abort () ;
161 exit(-1) ;
162 }
163 if (n2f > n2c) {
164 cout << "cfrchebpii: n2f > n2c : n2f = " << n2f << " , n2c = "
165 << n2c << endl ;
166 abort () ;
167 exit(-1) ;
168 }
169
170// Nombre de points pour la FFT:
171 int nm1 = nr - 1;
172 int nm1s2 = nm1 / 2;
173
174// Recherche des tables pour la FFT:
175 int* facto = facto_ini(nm1) ;
176 double* trigo = trigo_ini(nm1) ;
177
178// Recherche de la table des sin(psi) :
179 double* sinp = cheb_ini(nr);
180
181// Recherche de la table des points de collocations x_k :
182 double* x = chebimp_ini(nr);
183
184 // tableau de travail G et t1
185 // (la dimension nm1+2 = nr+1 est exigee par la routine fft991)
186 double* g = (double*)( malloc( (nm1+2)*sizeof(double) ) );
187 double* t1 = (double*)( malloc( (nm1+2)*sizeof(double) ) ) ;
188
189// Parametres pour la routine FFT991
190 int jump = 1 ;
191 int inc = 1 ;
192 int lot = 1 ;
193 int isign = - 1 ;
194
195// boucle sur phi et theta
196
197 int n2n3f = n2f * n3f ;
198 int n2n3c = n2c * n3c ;
199
200/*
201 * Borne de la boucle sur phi:
202 * si n1f = 1, on effectue la boucle une fois seulement.
203 * si n1f > 1, on va jusqu'a j = n1f-2 en sautant j = 1 (les coefficients
204 * j=n1f-1 et j=0 ne sont pas consideres car nuls).
205 */
206 int borne_phi = ( n1f > 1 ) ? n1f-1 : 1 ;
207
208 for (j=0; j< borne_phi; j++) {
209
210 if (j==1) continue ; // on ne traite pas le terme en sin(0 phi)
211
212 /************** Cas l impair ***************/
213
214 for (k=1; k<n2f; k+=2) {
215
216 int i0 = n2n3f * j + n3f * k ; // indice de depart
217 double* ff0 = ff + i0 ; // tableau des donnees a transformer
218
219 i0 = n2n3c * j + n3c * k ; // indice de depart
220 double* cf0 = cf + i0 ; // tableau resultat
221
222
223// Valeur en theta=0 de la partie antisymetrique de F, F_ :
224 double fmoins0 = 0.5 * ( ff0[nm1] - ff0[0] );
225
226// Fonction G(theta) = F+(theta) + F_(theta) sin(theta)
227//---------------------------------------------
228 for ( i = 1; i < nm1s2 ; i++ ) {
229// ... indice (dans le tableau g) du pt symetrique de theta par rapport a pi/2:
230 int isym = nm1 - i ;
231// ... indice (dans le tableau ff0) du point x correspondant a theta
232 int ix = nm1 - i ;
233// ... indice (dans le tableau ff0) du point x correspondant a sym(theta)
234 int ixsym = nm1 - isym ;
235
236// ... F+(psi)
237 double fp = 0.5 * ( ff0[ix] + ff0[ixsym] ) ;
238// ... F_(psi) sin(psi)
239 double fms = 0.5 * ( ff0[ix] - ff0[ixsym] ) * sinp[i] ;
240 g[i] = fp + fms ;
241 g[isym] = fp - fms ;
242 }
243//... cas particuliers:
244 g[0] = 0.5 * ( ff0[nm1] + ff0[0] );
245 g[nm1s2] = ff0[nm1s2];
246
247// Developpement de G en series de Fourier par une FFT
248//----------------------------------------------------
249
250 F77_fft991( g, t1, trigo, facto, &inc, &jump, &nm1, &lot, &isign) ;
251
252// Coefficients pairs du developmt. de Tchebyshev de f
253//----------------------------------------------------
254// Ces coefficients sont egaux aux coefficients en cosinus du developpement
255// de G en series de Fourier (le facteur 2 vient de la normalisation
256// de fft991) :
257
258 cf0[0] = g[0] ;
259 for (i=2; i<nm1; i += 2 ) cf0[i] = 2. * g[i] ;
260 cf0[nm1] = g[nm1] ;
261
262// Coefficients impairs du developmt. de Tchebyshev de f
263//------------------------------------------------------
264// 1. Coef. c'_k (recurrence amorcee a partir de zero)
265// Le +4. en facteur de g[i] est du a la normalisation de fft991
266// (si fft991 donnait reellement les coef. en sinus, il faudrait le
267// remplacer par un -2.)
268 cf0[1] = 0 ;
269 double som = 0;
270 for ( i = 3; i < nr; i += 2 ) {
271 cf0[i] = cf0[i-2] + 4. * g[i] ;
272 som += cf0[i] ;
273 }
274
275// 2. Calcul de c_1 :
276 double c1 = ( fmoins0 - som ) / nm1s2 ;
277
278// 3. Coef. c_k avec k impair:
279 cf0[1] = c1 ;
280 for ( i = 3; i < nr; i += 2 ) cf0[i] += c1 ;
281
282
283 } // fin de la boucle sur theta
284
285 /************** Cas l pair ***************/
286
287 for (k=0; k<n2f; k+=2) {
288 int i0 = n2n3f * j + n3f * k ; // indice de depart
289 double* ff0 = ff + i0 ; // tableau des donnees a transformer
290
291 i0 = n2n3c * j + n3c * k ; // indice de depart
292 double* cf0 = cf + i0 ; // tableau resultat
293
294// Multiplication de la fonction par x (pour la rendre paire)
295// NB: dans les commentaires qui suivent, on note h(x) = x f(x).
296// (Les valeurs de h dans l'espace des configurations sont stokees dans le
297// tableau cf0).
298 cf0[0] = 0 ;
299 for (i=1; i<nr; i++) cf0[i] = x[i] * ff0[i] ;
300
301// Valeur en psi=0 de la partie antisymetrique de F, F_ :
302 double fmoins0 = 0.5 * ( cf0[nm1] - cf0[0] );
303
304// Fonction G(theta) = F+(theta) + F_(theta) sin(theta)
305//---------------------------------------------
306 for ( i = 1; i < nm1s2 ; i++ ) {
307// ... indice (dans le tableau g) du pt symetrique de theta par rapport a pi/2:
308 int isym = nm1 - i ;
309// ... indice (dans le tableau cf0) du point x correspondant a theta
310 int ix = nm1 - i ;
311// ... indice (dans le tableau cf0) du point x correspondant a sym(theta)
312 int ixsym = nm1 - isym ;
313
314// ... F+(psi)
315 double fp = 0.5 * ( cf0[ix] + cf0[ixsym] ) ;
316// ... F_(psi) sin(psi)
317 double fms = 0.5 * ( cf0[ix] - cf0[ixsym] ) * sinp[i] ;
318 g[i] = fp + fms ;
319 g[isym] = fp - fms ;
320 }
321//... cas particuliers:
322 g[0] = 0.5 * ( cf0[nm1] + cf0[0] );
323 g[nm1s2] = cf0[nm1s2];
324
325// Developpement de G en series de Fourier par une FFT
326//----------------------------------------------------
327
328 F77_fft991( g, t1, trigo, facto, &inc, &jump, &nm1, &lot, &isign) ;
329
330// Coefficients pairs du developmt. de Tchebyshev de h
331//----------------------------------------------------
332// Ces coefficients sont egaux aux coefficients en cosinus du developpement
333// de G en series de Fourier (le facteur 2. vient de la normalisation
334// de fft991; si fft991 donnait reellement les coef. en cosinus, il faudrait le
335// remplacer par un +1.) :
336
337 cf0[0] = g[0] ;
338 for (i=2; i<nm1; i += 2 ) cf0[i] = 2. * g[i] ;
339 cf0[nm1] = g[nm1] ;
340
341// Coefficients impairs du developmt. de Tchebyshev de h
342//------------------------------------------------------
343// 1. Coef. c'_k (recurrence amorcee a partir de zero)
344// Le +4. en facteur de g[i] est du a la normalisation de fft991
345// (si fft991 donnait reellement les coef. en sinus, il faudrait le
346// remplacer par un -2.)
347 cf0[1] = 0 ;
348 double som = 0;
349 for ( i = 3; i < nr; i += 2 ) {
350 cf0[i] = cf0[i-2] + 4. * g[i] ;
351 som += cf0[i] ;
352 }
353
354// 2. Calcul de c_1 :
355 double c1 = ( fmoins0 - som ) / nm1s2 ;
356
357// 3. Coef. c_k avec k impair:
358 cf0[1] = c1 ;
359 for ( i = 3; i < nr; i += 2 ) cf0[i] += c1 ;
360
361// Coefficients de f en fonction de ceux de h
362//-------------------------------------------
363
364 cf0[0] = 2* cf0[0] ;
365 for (i=1; i<nm1; i++) {
366 cf0[i] = 2 * cf0[i] - cf0[i-1] ;
367 }
368 cf0[nm1] = 0 ;
369
370
371 } // fin de la boucle sur theta
372 } // fin de la boucle sur phi
373
374 // Menage
375 free (t1) ;
376 free (g) ;
377
378}
379}
380
Lorene prototypes.
Definition app_hor.h:64