LORENE
FFTW3/cfpcossini.C
1/*
2 * Copyright (c) 1999-2001 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 cfpcossini_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/FFTW3/cfpcossini.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 * Cas d'une fonction antisymetrique par la transformation phi --> phi + Pi
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 = Pi k/np 0 <= k <= np-1
44 * La convention de stokage utilisee est la suivante
45 * cf[ dim[2]*dim[1]*k + dim[2]*j + i ] = f(phi_k) 0 <= k <= np-1,
46 * les indices j et i correspondant respectivement
47 * a theta et a r.
48 * L'espace memoire correspondant au
49 * pointeur cf doit etre dim[0]*dim[1]*dim[2] et doit
50 * etre alloue avant l'appel a la routine.
51 *
52 * sortie: tableau des coefficients de la fonction suivant
53 * la convention de stokage
54 * cf[ dim[2]*dim[1]*k + dim[2]*j + i ] = c_k 0<= k < np,
55 * ou les indices j et i correspondent respectivement
56 * a theta et a r et ou les c_k sont les coefficients
57 * du developpement de f en series de Fourier:
58 *
59 * f(phi) = c_0 cos(phi) + c_2 sin(phi)
60 * + som_{l=1}^{np/2-1} c_{2l+1} cos( (2l+1) phi )
61 * + c_{2l+2} sin( (2l+1) phi ),
62 *
63 * En particulier: c_1 = 0 et c_{np+1} = 0
64 */
65
66/*
67 * $Id: cfpcossini.C,v 1.3 2014/10/13 08:53:18 j_novak Exp $
68 * $Log: cfpcossini.C,v $
69 * Revision 1.3 2014/10/13 08:53:18 j_novak
70 * Lorene classes and functions now belong to the namespace Lorene.
71 *
72 * Revision 1.2 2014/10/06 15:18:47 j_novak
73 * Modified #include directives to use c++ syntax.
74 *
75 * Revision 1.1 2004/12/21 17:06:02 j_novak
76 * Added all files for using fftw3.
77 *
78 * Revision 1.2 2002/10/16 14:36:43 j_novak
79 * Reorganization of #include instructions of standard C++, in order to
80 * use experimental version 3 of gcc.
81 *
82 * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
83 * LORENE
84 *
85 * Revision 2.1 2000/09/08 15:55:04 eric
86 * Premiere version testee.
87 *
88 * Revision 2.0 2000/09/07 15:15:06 eric
89 * *** empty log message ***
90 *
91 *
92 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/FFTW3/cfpcossini.C,v 1.3 2014/10/13 08:53:18 j_novak Exp $
93 *
94 */
95
96// Headers C
97#include <cstdlib>
98
99// Headers Lorene
100#include "headcpp.h"
101#include "proto.h"
102
103namespace Lorene {
104//*****************************************************************************
105
106void cfpcossini(const int* deg, const int* dim, double* cf) {
107
108 // Dimensions du tableau cf :
109 int n1 = dim[0] ;
110 int nt = dim[1] ;
111 int nr = dim[2] ;
112 int ntnr = nt * nr ;
113
114 // Nombres de degres de liberte en phi :
115 int np = deg[0] ;
116
117 // Tests de dimension:
118 if (np+2 > n1) {
119 cout << "cfpcossini: np+2 > n1 : np+2 = " << np+2 << " , n1 = "
120 << n1 << endl ;
121 abort() ;
122 }
123
124 // Tableau contenant les points de 0 a 2 Pi (et non seulement de 0 a Pi)
125 int np2 = 2*np ;
126 int deg2[] = {np2, nt, nr} ;
127 int dim2[] = {np2+2, nt, nr} ;
128
129 double* cf2 = new double[(np2+2)*nt*nr] ;
130
131 // Recopie des valeurs pour phi dans [0, Pi[ :
132 for (int k=0; k<np; k++) {
133 for (int ij = 0; ij <ntnr; ij++) {
134 cf2[k*ntnr + ij] = cf[k*ntnr + ij] ;
135 }
136 }
137
138 // Valeurs pour phi dans [Pi, 2Pi[ obtenues par antisymetrie:
139 int npntnr = np * ntnr ;
140 for (int k=0; k<np; k++) {
141 for (int ij = 0; ij <ntnr; ij++) {
142 cf2[npntnr + k*ntnr + ij] = - cf[k*ntnr + ij] ;
143 }
144 }
145
146 // Transformation de Fourier sur cf2 :
147 cfpcossin(deg2, dim2, cf2) ;
148
149 // Recopie des coefficients dans cf :
150
151 // Terme k=0 cos(phi)
152 for (int ij = 0; ij <ntnr; ij++) {
153 cf[ij] = cf2[2*ntnr + ij] ;
154 }
155
156 // Terme k=1
157 for (int ij = 0; ij <ntnr; ij++) {
158 cf[ntnr + ij] = 0 ;
159 }
160
161 // Terme k=2 sin(phi)
162 for (int ij = 0; ij <ntnr; ij++) {
163 cf[2*ntnr + ij] = cf2[3*ntnr + ij] ;
164 }
165
166 // Termes suivants:
167 for (int k=3; k<np; k+=2) {
168 int k2 = 2*(k-1) + 2 ;
169 // Terme en cosinus
170 for (int ij = 0; ij <ntnr; ij++) {
171 cf[k*ntnr + ij] = cf2[k2*ntnr + ij] ;
172 }
173 // Terme en sinus
174 k2++ ;
175 int k1 = k+1 ;
176 for (int ij = 0; ij <ntnr; ij++) {
177 cf[k1*ntnr + ij] = cf2[k2*ntnr + ij] ;
178 }
179 }
180
181 // Terme k=np+1
182 for (int ij = 0; ij <ntnr; ij++) {
183 cf[(np+1)*ntnr + ij] = 0 ;
184 }
185
186
187
188 // Menage
189 delete [] cf2 ;
190
191}
192}
Lorene prototypes.
Definition app_hor.h:64