LORENE
FFT991/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/FFT991/cfpcossin.C,v 1.4 2014/10/15 12:48:19 j_novak Exp $" ;
24
25/*
26 * Transformation de Fourier sur le premier indice d'un tableau 3-D
27 * par le biais de la routine FFT Fortran FFT991
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 3^q 5^r
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.4 2014/10/15 12:48:19 j_novak Exp $
72 * $Log: cfpcossin.C,v $
73 * Revision 1.4 2014/10/15 12:48:19 j_novak
74 * Corrected namespace declaration.
75 *
76 * Revision 1.3 2014/10/13 08:53:15 j_novak
77 * Lorene classes and functions now belong to the namespace Lorene.
78 *
79 * Revision 1.2 2014/10/06 15:18:44 j_novak
80 * Modified #include directives to use c++ syntax.
81 *
82 * Revision 1.1 2004/12/21 17:06:01 j_novak
83 * Added all files for using fftw3.
84 *
85 * Revision 1.4 2003/01/31 10:31:23 e_gourgoulhon
86 * Suppressed the directive #include <malloc.h> for malloc is defined
87 * in <stdlib.h>
88 *
89 * Revision 1.3 2002/10/16 14:36:43 j_novak
90 * Reorganization of #include instructions of standard C++, in order to
91 * use experimental version 3 of gcc.
92 *
93 * Revision 1.2 2002/09/09 13:00:39 e_gourgoulhon
94 * Modification of declaration of Fortran 77 prototypes for
95 * a better portability (in particular on IBM AIX systems):
96 * All Fortran subroutine names are now written F77_* and are
97 * defined in the new file C++/Include/proto_f77.h.
98 *
99 * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
100 * LORENE
101 *
102 * Revision 2.0 1999/02/22 15:48:58 hyc
103 * *** empty log message ***
104 *
105 *
106 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/FFT991/cfpcossin.C,v 1.4 2014/10/15 12:48:19 j_novak Exp $
107 *
108 */
109// headers du C
110#include <cstdlib>
111#include <cassert>
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 ) ;
121//*****************************************************************************
122
123void cfpcossin(const int* deg, const int* dim, double* cf)
124
125{
126
127int i, j, k, index ;
128
129// Dimensions du tableau cf :
130 int n1 = dim[0] ;
131 int n2 = dim[1] ;
132 int n3 = dim[2] ;
133
134// Nombres de degres de liberte en phi :
135 int np = deg[0] ;
136
137// Tests de dimension:
138 if (np+2 > n1) {
139 cout << "cfpcossin: np+2 > n1 : np+2 = " << np+2 << " , n1 = "
140 << n1 << endl ;
141 abort () ;
142 exit(-1) ;
143 }
144
145// Recherche des tables
146
147 int* facto = facto_ini(np) ;
148 double* trigo = trigo_ini(np) ;
149
150// Tableau de travail
151 double* t1 = (double*)( malloc( (np+2)*sizeof(double) ) ) ;
152
153// Parametres pour la routine FFT991
154 int jump = 1 ;
155 int inc = n2 * n3 ;
156 int lot = 1 ;
157 int isign = - 1 ;
158
159// boucle sur theta et r
160 for (j=0; j<n2; j++) {
161 for (k=0; k<n3; k++) {
162 index = n3 * j + k ;
163// FFT
164 double* debut = cf + index ;
165
166 F77_fft991( debut, t1, trigo, facto, &inc, &jump, &np,
167 &lot, &isign) ;
168 } // fin de la boucle sur r
169 } // fin de la boucle sur theta
170
171// Normalisation des cosinus
172 int n2n3 = n2 * n3 ;
173 for (i=2; i<np; i += 2 ) {
174 for (j=0; j<n2; j++) {
175 for (k=0; k<n3; k++) {
176 index = n2n3 * i + n3 * j + k ;
177 cf[index] *= 2. ;
178 }
179 }
180 }
181
182// Normalisation des sinus
183 for (i=1; i<np+1; i += 2 ) {
184 for (j=0; j<n2; j++) {
185 for (k=0; k<n3; k++) {
186 index = n2n3 * i + n3 * j + k ;
187 cf[index] *= -2. ;
188 }
189 }
190 }
191
192 // Menage
193 free(t1) ;
194
195}
196}
Lorene prototypes.
Definition app_hor.h:64