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 
23 char 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:
118 namespace Lorene {
119 int* facto_ini(int ) ;
120 double* trigo_ini(int ) ;
121 //*****************************************************************************
122 
123 void cfpcossin(const int* deg, const int* dim, double* cf)
124 
125 {
126 
127 int 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
Lorene prototypes.
Definition: app_hor.h:64