LORENE
FFTW3/cfpcossin.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 cfpcossin_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/FFTW3/cfpcossin.C,v 1.3 2014/10/13 08:53:18 j_novak Exp $" ;
24
25/*
26 * Transformation de Fourier sur le premier indice d'un tableau 3-D
27 * par le biais de la bibliotheque fftw.
28 *
29 * Entree:
30 * -------
31 * int* deg : tableau du nombre effectif de degres de liberte dans chacune
32 * des 3 dimensions: le nombre de points de collocation
33 * en phi est np = deg[0] et doit etre de la forme
34 * np = 2*p
35 * int* dim : tableau du nombre d'elements de ff dans chacune des trois
36 * dimensions.
37 * On doit avoir dim[0] >= deg[0] + 2 = np + 2.
38 *
39 * Entree/Sortie :
40 * ---------------
41 * double* cf : entree: tableau des valeurs de la fonction f aux np points de
42 * de collocation
43 * phi_k = T k/np 0 <= k <= np-1
44 * T etant la periode de f. La convention de stokage
45 * utilisee est la suivante
46 * cf[ dim[2]*dim[1]*k + dim[2]*j + i ] = f(phi_k) 0 <= k <= np-1,
47 * les indices j et i correspondant respectivement
48 * a theta et a r.
49 * L'espace memoire correspondant au
50 * pointeur cf doit etre dim[0]*dim[1]*dim[2] et doit
51 * etre alloue avant l'appel a la routine.
52 *
53 * sortie: tableau des coefficients de la fonction suivant
54 * la convention de stokage
55 * cf[ dim[2]*dim[1]*k + dim[2]*j + i ] = c_k 0<= k <= np,
56 * ou les indices j et i correspondent respectivement
57 * a theta et a r et ou les c_k sont les coefficients
58 * du developpement de f en series de Fourier:
59 *
60 * f(phi) = som_{l=0}^{np/2} c_{2l} cos( 2 pi/T l phi )
61 * + c_{2l+1} sin( 2 pi/T l phi ),
62 *
63 * ou T est la periode de f.
64 * En particulier cf[1] = 0.
65 * De plus, cf[np+1] n'est pas egal a c_{np+1}
66 * mais a zero.
67 *
68 */
69
70/*
71 * $Id: cfpcossin.C,v 1.3 2014/10/13 08:53:18 j_novak Exp $
72 * $Log: cfpcossin.C,v $
73 * Revision 1.3 2014/10/13 08:53:18 j_novak
74 * Lorene classes and functions now belong to the namespace Lorene.
75 *
76 * Revision 1.2 2014/10/06 15:18:47 j_novak
77 * Modified #include directives to use c++ syntax.
78 *
79 * Revision 1.1 2004/12/21 17:06:02 j_novak
80 * Added all files for using fftw3.
81 *
82 * Revision 1.4 2003/01/31 10:31:23 e_gourgoulhon
83 * Suppressed the directive #include <malloc.h> for malloc is defined
84 * in <stdlib.h>
85 *
86 * Revision 1.3 2002/10/16 14:36:43 j_novak
87 * Reorganization of #include instructions of standard C++, in order to
88 * use experimental version 3 of gcc.
89 *
90 * Revision 1.2 2002/09/09 13:00:39 e_gourgoulhon
91 * Modification of declaration of Fortran 77 prototypes for
92 * a better portability (in particular on IBM AIX systems):
93 * All Fortran subroutine names are now written F77_* and are
94 * defined in the new file C++/Include/proto_f77.h.
95 *
96 * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
97 * LORENE
98 *
99 * Revision 2.0 1999/02/22 15:48:58 hyc
100 * *** empty log message ***
101 *
102 *
103 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/FFTW3/cfpcossin.C,v 1.3 2014/10/13 08:53:18 j_novak Exp $
104 *
105 */
106// headers du C
107#include <cstdlib>
108#include <fftw3.h>
109
110//Lorene prototypes
111#include "tbl.h"
112
113// Prototypage des sous-routines utilisees:
114namespace Lorene {
115fftw_plan prepare_fft(int, Tbl*&) ;
116//*****************************************************************************
117
118void cfpcossin(const int* deg, const int* dim, double* cf)
119{
120// Dimensions du tableau cf :
121 int n1 = dim[0] ;
122 int n2 = dim[1] ;
123 int n3 = dim[2] ;
124
125// Nombres de degres de liberte en phi :
126 int np = deg[0] ;
127
128// Tests de dimension:
129 if (np+2 > n1) {
130 cout << "cfpcossin: np+2 > n1 : np+2 = " << np+2 << " , n1 = "
131 << n1 << endl ;
132 abort () ;
133 exit(-1) ;
134 }
135
136// Recherche des tables
137 Tbl* pg = 0x0 ;
138 fftw_plan p = prepare_fft(np, pg) ;
139 Tbl& g = *pg ;
140
141 int index = 0 ;
142 int n2n3 = n2 * n3 ;
143 double fac = 2./double(np) ;
144
145// boucle sur theta et r
146 for (int j=0; j<n2; j++) {
147 for (int k=0; k<n3; k++) {
148 index = n3 * j + k ;
149// FFT
150 double* debut = cf + index ;
151 double* tab = g.t ;
152 for (int i=0; i<np; i++) {
153 *tab = *debut ;
154 tab++;
155 debut += n2n3 ;
156 }
157
158 fftw_execute(p) ;
159
160 debut = cf+index ;
161 double* pcos = g.t ;
162 double* psin = g.t + np - 1 ;
163 (*debut) = (*pcos)/double(np) ;
164 debut += n2n3 ; pcos++ ;
165 (*debut) = 0. ;
166 debut += n2n3 ;
167 for (int i=1; i<np/2; i++){
168 *debut = (*pcos)*fac ;
169 debut += n2n3 ; pcos++ ;
170 *debut = -(*psin)*fac ;
171 debut += n2n3 ; psin-- ;
172 }
173 (*debut) = (*pcos)/double(np) ;
174 debut += n2n3 ;
175 (*debut) = 0. ;
176
177 } // fin de la boucle sur r
178 } // fin de la boucle sur theta
179
180
181}
182}
Lorene prototypes.
Definition app_hor.h:64