LORENE
chb_sin_legmi.C
1/*
2 * Copyright (c) 1999-2001 Eric Gourgoulhon
3 * Copyright (c) 2009 Jerome Novak
4 *
5 * This file is part of LORENE.
6 *
7 * LORENE is free software; you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation; either version 2 of the License, or
10 * (at your option) any later version.
11 *
12 * LORENE is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with LORENE; if not, write to the Free Software
19 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 *
21 */
22
23
24char chb_sin_legmi_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/chb_sin_legmi.C,v 1.3 2014/10/13 08:53:11 j_novak Exp $" ;
25
26/*
27 * Calcule les coefficients du developpement (suivant theta) en fonctions
28 * associees de Legendre P_l^m(cos(theta)) a partir des coefficients du
29 * developpement en sin(j*theta)
30 * representant une fonction 3-D antisymetrique par le retournement
31 * (x, y, z) --> (-x, -y, z).
32 *
33 * Entree:
34 * -------
35 * const int* deg : tableau du nombre effectif de degres de liberte dans chacune
36 * des 3 dimensions:
37 * deg[0] = np : nombre de points de collocation en phi
38 * deg[1] = nt : nombre de points de collocation en theta
39 * deg[2] = nr : nombre de points de collocation en r
40 *
41 * const double* cfi : tableau des coefficients c_j du develop. en cos defini
42 * comme suit (a r et phi fixes)
43 *
44 * f(theta) = som_{j=0}^{nt-1} c_j sin( j theta )
45 *
46 * L'espace memoire correspondant au pointeur cfi doit etre
47 * nr*nt*(np+2) et doit avoir ete alloue avant
48 * l'appel a la routine.
49 * Le coefficient c_j (0 <= j <= nt-2) doit etre stoke dans le
50 * tableau cfi comme suit
51 * c_j = cfi[ nr*nt* k + i + nr* j ]
52 * ou k et i sont les indices correspondant a
53 * phi et r respectivement.
54 *
55 * Sortie:
56 * -------
57 * double* cfo : tableau des coefficients a_j du develop. en fonctions de
58 * Legendre associees P_l^m (m impair)
59 *
60 * f(theta) =
61 * som_{l=m}^{nt-1} a_j P_j^m( cos(theta) )
62 *
63 * avec m impair.
64 *
65 * P_l^m(x) represente la fonction de Legendre associee
66 * de degre l et d'ordre m normalisee de facon a ce que
67 *
68 * int_0^pi [ P_l^m(cos(theta)) ]^2 sin(theta) dtheta = 1
69 *
70 * L'espace memoire correspondant au pointeur cfo doit etre
71 * nr*nt*(np+2) et doit avoir ete alloue avant
72 * l'appel a la routine.
73 * Le coefficient a_j (0 <= j <= nt-1) est stoke dans le
74 * tableau cfo comme suit
75 * a_j = cfo[ nr*nt* k + i + nr* j ]
76 * ou k et i sont les indices correspondant a phi et r
77 * respectivement: m = 2( k/2 ).
78 * NB: pour j< m, a_j = 0
79 *
80 * NB:
81 * ---
82 * Il n'est pas possible d'avoir le pointeur cfo egal a cfi.
83 */
84
85/*
86 * $Id: chb_sin_legmi.C,v 1.3 2014/10/13 08:53:11 j_novak Exp $
87 * $Log: chb_sin_legmi.C,v $
88 * Revision 1.3 2014/10/13 08:53:11 j_novak
89 * Lorene classes and functions now belong to the namespace Lorene.
90 *
91 * Revision 1.2 2014/10/06 15:16:01 j_novak
92 * Modified #include directives to use c++ syntax.
93 *
94 * Revision 1.1 2009/10/23 12:54:47 j_novak
95 * New base T_LEG_MI
96 *
97 *
98 *
99 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/chb_sin_legmi.C,v 1.3 2014/10/13 08:53:11 j_novak Exp $
100 *
101 */
102
103// headers du C
104#include <cassert>
105#include <cstdlib>
106
107// Prototypage
108#include "headcpp.h"
109#include "proto.h"
110
111namespace Lorene {
112//******************************************************************************
113
114void chb_sin_legmi(const int* deg , const double* cfi, double* cfo) {
115
116int k2, l, jmin, j, i, m ;
117
118 // Nombres de degres de liberte en phi et theta :
119 int np = deg[0] ;
120 int nt = deg[1] ;
121 int nr = deg[2] ;
122
123 assert(np < 4*nt) ;
124 assert( cfi != cfo ) ;
125
126 // Tableau de travail
127 double* som = new double[nr] ;
128
129 // Recherche de la matrice de passage sin --> Legendre
130 double* aa = mat_sin_legmi(np, nt) ;
131
132 // Increment en m pour la matrice aa :
133 int maa = nt * nt ;
134
135 // Pointeurs de travail :
136 double* resu = cfo ;
137 const double* cc = cfi ;
138
139 // Increment en phi :
140 int ntnr = nt * nr ;
141
142 // Indice courant en phi :
143 int k = 0 ;
144
145 // Cas k=0 (m=1 : cos(phi))
146 // ------------------------
147
148 // Cas l=0 : a_l = 0
149 for (i=0; i<nr; i++) {
150 *resu = 0. ;
151 resu++ ;
152 }
153
154 // ... produit matriciel
155 for (l=1; l<nt-1; l++) {
156 for (i=0; i<nr; i++) {
157 som[i] = 0 ;
158 }
159
160 jmin = l ; // pour m=1, aa_lj = 0 pour j<l
161
162 for (j=jmin; j<nt-1; j++) {
163 double amlj = aa[nt*l + j] ;
164 for (i=0; i<nr; i++) {
165 som[i] += amlj * cc[nr*j + i] ;
166 }
167 }
168
169 for (i=0; i<nr; i++) {
170 *resu = som[i] ;
171 resu++ ;
172 }
173
174 } // fin de la boucle sur l
175
176 // Dernier coef en l=nt-1 mis a zero pour le cas m impair :
177 for (i=0; i<nr; i++) {
178 *resu = 0 ;
179 resu++ ;
180 }
181
182 // Special case np=1 (axisymmetry)
183 // -------------------------------
184 if (np==1) {
185 for (i=0; i<2*ntnr; i++) {
186 *resu = 0 ;
187 resu++ ;
188 }
189 delete [] som ;
190 return ;
191 }
192
193
194 // On passe au phi suivant :
195 cc = cc + ntnr ;
196 k++ ;
197
198 // Cas k=1 : tout est mis a zero
199 // -----------------------------
200
201 for (l=0; l<nt; l++) {
202 for (i=0; i<nr; i++) {
203 *resu = 0 ;
204 resu++ ;
205 }
206 }
207
208 // On passe au phi suivant :
209 cc = cc + ntnr ;
210 k++ ;
211
212 // Cas k=2 (m=1 : sin(phi))
213 // ------------------------
214
215 // Cas l=0 : a_l = 0
216 for (i=0; i<nr; i++) {
217 *resu = 0. ;
218 resu++ ;
219 }
220
221 // ... produit matriciel
222 for (l=1; l<nt-1; l++) {
223 for (i=0; i<nr; i++) {
224 som[i] = 0 ;
225 }
226
227 jmin = l ; // pour m=1, aa_lj = 0 pour j<l
228
229 for (j=jmin; j<nt-1; j++) {
230 double amlj = aa[nt*l + j] ;
231 for (i=0; i<nr; i++) {
232 som[i] += amlj * cc[nr*j + i] ;
233 }
234 }
235
236 for (i=0; i<nr; i++) {
237 *resu = som[i] ;
238 resu++ ;
239 }
240
241 } // fin de la boucle sur l
242
243 // Dernier coef en l=nt-1 mis a zero pour le cas m impair :
244 for (i=0; i<nr; i++) {
245 *resu = 0 ;
246 resu++ ;
247 }
248
249 // On passe au phi suivant :
250 cc = cc + ntnr ;
251 k++ ;
252
253 // On passe au m suivant
254 aa += maa ; // pointeur sur la nouvelle matrice de passage
255
256 // Cas k >= 3
257 // ----------
258
259 for (m=3; m < np ; m+=2) {
260
261 for (k2=0; k2 < 2; k2++) { // k2=0 : cos(m phi) ; k2=1 : sin(m phi)
262 int lmax = (m<nt-1 ? m : nt -1) ;
263 for (l=0; l<lmax; l++) {
264 for (i=0; i<nr; i++) {
265 *resu = 0 ;
266 resu++ ;
267 }
268 }
269
270 // ... produit matriciel
271 for (l=m; l<nt-1; l++) {
272 for (i=0; i<nr; i++) {
273 som[i] = 0 ;
274 }
275
276 jmin = 1 ;
277
278 for (j=jmin; j<nt-1; j++) {
279 double amlj = aa[nt*l + j] ;
280 for (i=0; i<nr; i++) {
281 som[i] += amlj * cc[nr*j + i] ;
282 }
283 }
284
285 for (i=0; i<nr; i++) {
286 *resu = som[i] ;
287 resu++ ;
288 }
289
290 } // fin de la boucle sur l
291
292 // Dernier coef en l=nt-1 mis a zero pour le cas m impair :
293 for (i=0; i<nr; i++) {
294 *resu = 0 ;
295 resu++ ;
296 }
297
298 // On passe au phi suivant :
299 cc = cc + ntnr ;
300 k++ ;
301
302 } // fin de la boucle sur k2
303
304 // On passe a l'harmonique en phi suivante :
305
306 aa += maa ; // pointeur sur la nouvelle matrice de passage
307
308 } // fin de la boucle (m) sur phi
309
310 // Cas k=np+1 : tout est mis a zero
311 // --------------------------------
312
313 for (l=0; l<nt; l++) {
314 for (i=0; i<nr; i++) {
315 *resu = 0 ;
316 resu++ ;
317 }
318 }
319
320
321//## verif :
322 assert(resu == cfo + (np+2)*ntnr) ;
323
324 // Menage
325 delete [] som ;
326
327}
328}
Lorene prototypes.
Definition app_hor.h:64