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 
23 char 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:
114 namespace Lorene {
115 fftw_plan prepare_fft(int, Tbl*&) ;
116 //*****************************************************************************
117 
118 void 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