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 
28 char 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 
88 namespace Lorene {
89 Tbl _sh_pvect_r_pas_prevu (int, int, double) ;
90 Tbl _sh_pvect_r_cheb (int, int, double) ;
91 Tbl _sh_pvect_r_chebp (int, int, double) ;
92 Tbl _sh_pvect_r_chebi (int, int, double) ;
93 Tbl _sh_pvect_r_chebu (int, int, double) ;
94 
95  //------------------------------------
96  // Routine pour les cas non prevus --
97  //------------------------------------
98 Tbl _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 
115 Tbl _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 
215 Tbl _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 
281 Tbl _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 
351 Tbl _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 
414 Tbl 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 }
R_CHEB
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:348
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
R_CHEBI
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
TRA_R
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
Definition: type_parite.h:158
Lorene::cos
Cmp cos(const Cmp &)
Cosine.
Definition: cmp_math.C:94
MAX_BASE
#define MAX_BASE
Nombre max. de bases differentes.
Definition: type_parite.h:144
R_CHEBU
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
Definition: type_parite.h:180
Lorene::sin
Cmp sin(const Cmp &)
Sine.
Definition: cmp_math.C:69