LORENE
sh_pvect_r.C
1/*
2 * Methods for computing the homogeneous solutions for Ope_pois_vect_r.
3 *
4 * (see file Ope_elementary for documentation).
5 *
6 */
7
8/*
9 * Copyright (c) 2004 Jerome Novak
10 *
11 * This file is part of LORENE.
12 *
13 * LORENE is free software; you can redistribute it and/or modify
14 * it under the terms of the GNU General Public License version 2
15 * as published by the Free Software Foundation.
16 *
17 * LORENE is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU General Public License for more details.
21 *
22 * You should have received a copy of the GNU General Public License
23 * along with LORENE; if not, write to the Free Software
24 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25 *
26 */
27
28char sh_pvect_r_C[] = "$Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_pois_vect_r/sh_pvect_r.C,v 1.5 2014/10/13 08:53:34 j_novak Exp $" ;
29
30/*
31 * $Id: sh_pvect_r.C,v 1.5 2014/10/13 08:53:34 j_novak Exp $
32 * $Log: sh_pvect_r.C,v $
33 * Revision 1.5 2014/10/13 08:53:34 j_novak
34 * Lorene classes and functions now belong to the namespace Lorene.
35 *
36 * Revision 1.4 2014/10/06 15:16:13 j_novak
37 * Modified #include directives to use c++ syntax.
38 *
39 * Revision 1.3 2004/12/23 16:00:10 j_novak
40 * Modif. comments.
41 *
42 * Revision 1.2 2004/05/10 15:36:42 j_novak
43 * Corrected a missing #include
44 *
45 * Revision 1.1 2004/05/10 15:28:22 j_novak
46 * First version of functions for the solution of the r-component of the
47 * vector Poisson equation.
48 *
49 *
50 * $Header: /cvsroot/Lorene/C++/Source/Ope_elementary/Ope_pois_vect_r/sh_pvect_r.C,v 1.5 2014/10/13 08:53:34 j_novak Exp $
51 *
52 */
53
54//fichiers includes
55#include <cstdlib>
56#include <cmath>
57
58#include "matrice.h"
59#include "type_parite.h"
60#include "proto.h"
61
62/*
63 *
64 * Renvoie une ou 2 solution homogene
65 * Si l <> 0 :
66 * Si base_r = R_CHEB deux solutions (x+echelle)^(l-1) dans (0, *) et
67 * 1/(x+echelle)^(l+2) dans (1, *)
68 * Si base_r = R_CHEBU 1 solution (x-1)^(l+2) dans (*)
69 * Si base_r = R_CHEBP ou R_CHEBI x^(l-1) dans (*)
70 *
71 * Si l = 0:
72 * Si base_r = R_CHEB deux solutions (x+echelle) dans (0, *) et
73 * 1/(x+echelle)^(l+2) dans (1, *)
74 * Si base_r = R_CHEBU 1 solution (x-1)^(l+2) dans (*)
75 * Si base_r = R_CHEBP ou R_CHEBI x dans (*)
76 *
77 * Entree :
78 * n : nbre de points en r
79 * l : nbre quantique associe
80 * echelle : cf ci-dessus, utile que dans le cas R_CHEB
81 * base_r : base de decomposition
82 *
83 * Sortie :
84 * Tbl contenant les coefficient de la ou des solution homogenes
85 *
86 */
87
88namespace Lorene {
89Tbl _sh_pvect_r_pas_prevu (int, int, double) ;
90Tbl _sh_pvect_r_cheb (int, int, double) ;
91Tbl _sh_pvect_r_chebp (int, int, double) ;
92Tbl _sh_pvect_r_chebi (int, int, double) ;
93Tbl _sh_pvect_r_chebu (int, int, double) ;
94
95 //------------------------------------
96 // Routine pour les cas non prevus --
97 //------------------------------------
98Tbl _sh_pvect_r_pas_prevu (int n, int l, double echelle) {
99
100 cout << " Solution homogene pas prevue ..... : "<< endl ;
101 cout << " N : " << n << endl ;
102 cout << " l : " << l << endl ;
103 cout << " echelle : " << echelle << endl ;
104 abort() ;
105 exit(-1) ;
106 Tbl res(1) ;
107 return res;
108}
109
110
111 //-------------------
112 //-- R_CHEB ------
113 //-------------------
114
115Tbl _sh_pvect_r_cheb (int n, int l, double echelle) {
116
117 const int nmax = 200 ; // Nombre de Tbl stockes
118 static Tbl* tab[nmax] ; // les Tbl calcules
119 static int nb_dejafait = 0 ; // nbre de Tbl calcules
120 static int l_dejafait[nmax] ;
121 static int nr_dejafait[nmax] ;
122 static double vieux_echelle = 0;
123
124 // Si on a change l'echelle : on detruit tout :
125 if (vieux_echelle != echelle) {
126 for (int i=0 ; i<nb_dejafait ; i++) {
127 l_dejafait[i] = -1 ;
128 nr_dejafait[i] = -1 ;
129 delete tab[i] ;
130 }
131 nb_dejafait = 0 ;
132 vieux_echelle = echelle ;
133 }
134
135 int indice = -1 ;
136
137 // On determine si la matrice a deja ete calculee :
138 for (int conte=0 ; conte<nb_dejafait ; conte ++)
139 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
140 indice = conte ;
141
142 // Calcul a faire :
143 if (indice == -1) {
144 if (nb_dejafait >= nmax) {
145 cout << "_sh_pvect_r_cheb : trop de Tbl" << endl ;
146 abort() ;
147 exit (-1) ;
148 }
149
150
151 l_dejafait[nb_dejafait] = l ;
152 nr_dejafait[nb_dejafait] = n ;
153
154 // assert (l < n) ;
155
156 Tbl res(2, n) ;
157 res.set_etat_qcq() ;
158 double* coloc = new double[n] ;
159
160 int * deg = new int[3] ;
161 deg[0] = 1 ;
162 deg[1] = 1 ;
163 deg[2] = n ;
164
165 //Construction de la premiere solution homogene :
166 // cad celle polynomiale.
167
168 if (l==0) {
169 for (int i=0 ; i<n ; i++)
170 coloc[i] = echelle - cos(M_PI*i/(n-1)) ;
171
172 cfrcheb(deg, deg, coloc, deg, coloc) ;
173 for (int i=0 ; i<n ;i++)
174 res.set(0, i) = coloc[i] ;
175 }
176 else if (l==1) {
177 res.set(0, 0) = 1 ;
178 for (int i=1 ; i<n ; i++)
179 res.set(0, i) = 0 ;
180 }
181 else {
182
183 for (int i=0 ; i<n ; i++)
184 coloc[i] = pow(echelle-cos(M_PI*i/(n-1)), double(l-1)) ;
185
186 cfrcheb(deg, deg, coloc, deg, coloc) ;
187 for (int i=0 ; i<n ;i++)
188 res.set(0, i) = coloc[i] ;
189 }
190
191
192 // construction de la seconde solution homogene :
193 // cad celle fractionnelle.
194 for (int i=0 ; i<n ; i++)
195 coloc[i] = 1/pow(echelle-cos(M_PI*i/(n-1)), double(l+2)) ;
196
197 cfrcheb(deg, deg, coloc, deg, coloc) ;
198 for (int i=0 ; i<n ;i++)
199 res.set(1, i) = coloc[i] ;
200
201 delete [] coloc ;
202 delete [] deg ;
203 tab[nb_dejafait] = new Tbl(res) ;
204 nb_dejafait ++ ;
205 return res ;
206 }
207
208 else return *tab[indice] ;
209}
210
211 //-------------------
212 //-- R_CHEBP ------
213 //-------------------
214
215Tbl _sh_pvect_r_chebp (int n, int l, double) {
216
217 const int nmax = 200 ; // Nombre de Tbl stockes
218 static Tbl* tab[nmax] ; // les Tbl calcules
219 static int nb_dejafait = 0 ; // nbre de Tbl calcules
220 static int l_dejafait[nmax] ;
221 static int nr_dejafait[nmax] ;
222
223 int indice = -1 ;
224
225 // On determine si la matrice a deja ete calculee :
226 for (int conte=0 ; conte<nb_dejafait ; conte ++)
227 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
228 indice = conte ;
229
230 // Calcul a faire :
231 if (indice == -1) {
232 if (nb_dejafait >= nmax) {
233 cout << "_sh_pvect_r_chebp : trop de Tbl" << endl ;
234 abort() ;
235 exit (-1) ;
236 }
237
238
239 l_dejafait[nb_dejafait] = l ;
240 nr_dejafait[nb_dejafait] = n ;
241
242 Tbl res(n) ;
243 res.set_etat_qcq() ;
244 double* coloc = new double[n] ;
245
246 int * deg = new int[3] ;
247 deg[0] = 1 ;
248 deg[1] = 1 ;
249 deg[2] = n ;
250
251 assert (div(l, 2).rem == 1) ;
252 if (l==1) {
253 res.set(0) = 1 ;
254 for (int i=1 ; i<n ; i++)
255 res.set(i) = 0 ;
256 }
257 else {
258 for (int i=0 ; i<n ; i++)
259 coloc[i] = pow(sin(M_PI*i/2/(n-1)), double(l-1)) ;
260
261 cfrchebp(deg, deg, coloc, deg, coloc) ;
262 for (int i=0 ; i<n ;i++)
263 res.set(i) = coloc[i] ;
264 }
265
266 delete [] coloc ;
267 delete [] deg ;
268 tab[nb_dejafait] = new Tbl(res) ;
269 nb_dejafait ++ ;
270 return res ;
271 }
272
273 else return *tab[indice] ;
274}
275
276
277 //-------------------
278 //-- R_CHEBI -----
279 //-------------------
280
281Tbl _sh_pvect_r_chebi (int n, int l, double) {
282
283 const int nmax = 200 ; // Nombre de Tbl stockes
284 static Tbl* tab[nmax] ; // les Tbl calcules
285 static int nb_dejafait = 0 ; // nbre de Tbl calcules
286 static int l_dejafait[nmax] ;
287 static int nr_dejafait[nmax] ;
288
289 int indice = -1 ;
290
291 // On determine si la matrice a deja ete calculee :
292 for (int conte=0 ; conte<nb_dejafait ; conte ++)
293 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
294 indice = conte ;
295
296 // Calcul a faire :
297 if (indice == -1) {
298 if (nb_dejafait >= nmax) {
299 cout << "_sh_pvect_r_chebi : trop de Tbl" << endl ;
300 abort() ;
301 exit (-1) ;
302 }
303
304
305 l_dejafait[nb_dejafait] = l ;
306 nr_dejafait[nb_dejafait] = n ;
307
308
309 assert (div(l, 2).rem == 0) ;
310 // assert (l < 2*n) ;
311
312 Tbl res(n) ;
313 res.set_etat_qcq() ;
314 double* coloc = new double[n] ;
315
316 int * deg = new int[3] ;
317 deg[0] = 1 ;
318 deg[1] = 1 ;
319 deg[2] = n ;
320
321 if (l==0) {
322 res.set(0) = 1 ;
323 for (int i=1 ; i<n ; i++)
324 res.set(i) = 0 ;
325 }
326 else {
327 for (int i=0 ; i<n ; i++)
328 coloc[i] = pow(sin(M_PI*i/2/(n-1)), double(l-1)) ;
329
330 cfrchebi(deg, deg, coloc, deg, coloc) ;
331 for (int i=0 ; i<n ;i++)
332 res.set(i) = coloc[i] ;
333 }
334
335 delete [] coloc ;
336 delete [] deg ;
337 tab[nb_dejafait] = new Tbl(res) ;
338 nb_dejafait ++ ;
339 return res ;
340 }
341
342 else return *tab[indice] ;
343}
344
345
346
347 //-------------------
348 //-- R_CHEBU -----
349 //-------------------
350
351Tbl _sh_pvect_r_chebu (int n, int l, double) {
352
353 const int nmax = 200 ; // Nombre de Tbl stockes
354 static Tbl* tab[nmax] ; // les Tbl calcules
355 static int nb_dejafait = 0 ; // nbre de Tbl calcules
356 static int l_dejafait[nmax] ;
357 static int nr_dejafait[nmax] ;
358
359 int indice = -1 ;
360
361 // On determine si la matrice a deja ete calculee :
362 for (int conte=0 ; conte<nb_dejafait ; conte ++)
363 if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
364 indice = conte ;
365
366 // Calcul a faire :
367 if (indice == -1) {
368 if (nb_dejafait >= nmax) {
369 cout << "_sh_pvect_r_chebu : trop de Tbl" << endl ;
370 abort() ;
371 exit (-1) ;
372 }
373
374
375 l_dejafait[nb_dejafait] = l ;
376 nr_dejafait[nb_dejafait] = n ;
377
378
379 // assert (l < n-1) ;
380 Tbl res(n) ;
381 res.set_etat_qcq() ;
382 double* coloc = new double[n] ;
383
384 int * deg = new int[3] ;
385 deg[0] = 1 ;
386 deg[1] = 1 ;
387 deg[2] = n ;
388
389 for (int i=0 ; i<n ; i++)
390 coloc[i] = pow(-1-cos(M_PI*i/(n-1)), double(l+2)) ;
391
392 cfrcheb(deg, deg, coloc, deg, coloc) ;
393 for (int i=0 ; i<n ;i++)
394 res.set(i) = coloc[i] ;
395
396 delete [] coloc ;
397 delete [] deg ;
398 tab[nb_dejafait] = new Tbl(res) ;
399 nb_dejafait ++ ;
400 return res ;
401 }
402
403 else return *tab[indice] ;
404}
405
406
407
408
409 //-------------------
410 //-- Fonction ----
411 //-------------------
412
413
414Tbl sh_pvect_r(int n, int l, double echelle, int base_r) {
415
416 // Routines de derivation
417 static Tbl (*sh_pvect_r[MAX_BASE])(int, int, double) ;
418 static int nap = 0 ;
419
420 // Premier appel
421 if (nap==0) {
422 nap = 1 ;
423 for (int i=0 ; i<MAX_BASE ; i++) {
424 sh_pvect_r[i] = _sh_pvect_r_pas_prevu ;
425 }
426 // Les routines existantes
427 sh_pvect_r[R_CHEB >> TRA_R] = _sh_pvect_r_cheb ;
428 sh_pvect_r[R_CHEBU >> TRA_R] = _sh_pvect_r_chebu ;
429 sh_pvect_r[R_CHEBP >> TRA_R] = _sh_pvect_r_chebp ;
430 sh_pvect_r[R_CHEBI >> TRA_R] = _sh_pvect_r_chebi ;
431 }
432
433 Tbl res(sh_pvect_r[base_r](n, l, echelle)) ;
434 return res ;
435}
436}
Cmp sin(const Cmp &)
Sine.
Definition cmp_math.C:69
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:348
Cmp cos(const Cmp &)
Cosine.
Definition cmp_math.C:94
#define MAX_BASE
Nombre max. de bases differentes.
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
#define R_CHEBI
base de Cheb. impaire (rare) seulement
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
#define R_CHEB
base de Chebychev ordinaire (fin)
#define R_CHEBP
base de Cheb. paire (rare) seulement
Lorene prototypes.
Definition app_hor.h:64