LORENE
solh.C
1 /*
2  * Copyright (c) 1999-2001 Philippe Grandclement
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 
24 
25 /*
26  * $Id: solh.C,v 1.11 2016/12/05 16:18:10 j_novak Exp $
27  * $Log: solh.C,v $
28  * Revision 1.11 2016/12/05 16:18:10 j_novak
29  * Suppression of some global variables (file names, loch, ...) to prevent redefinitions
30  *
31  * Revision 1.10 2014/10/13 08:53:30 j_novak
32  * Lorene classes and functions now belong to the namespace Lorene.
33  *
34  * Revision 1.9 2014/10/06 15:16:10 j_novak
35  * Modified #include directives to use c++ syntax.
36  *
37  * Revision 1.8 2008/02/18 13:53:43 j_novak
38  * Removal of special indentation instructions.
39  *
40  * Revision 1.7 2007/12/20 09:11:09 jl_cornou
41  * Correction of an error in op_sxpun about Jacobi(0,2) polynomials
42  *
43  * Revision 1.6 2007/12/13 15:48:46 jl_cornou
44  * *** empty log message ***
45  *
46  * Revision 1.5 2007/12/12 12:30:49 jl_cornou
47  * *** empty log message ***
48  *
49  * Revision 1.4 2004/10/05 15:44:21 j_novak
50  * Minor speed enhancements.
51  *
52  * Revision 1.3 2003/01/31 08:49:58 e_gourgoulhon
53  * Increased the number nmax of stored matrices from 100 to 200.
54  *
55  * Revision 1.2 2002/10/16 14:37:12 j_novak
56  * Reorganization of #include instructions of standard C++, in order to
57  * use experimental version 3 of gcc.
58  *
59  * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
60  * LORENE
61  *
62  * Revision 2.13 2000/01/18 14:15:50 phil
63  * enleve assert sur nobre de points min en r
64  *
65  * Revision 2.12 2000/01/04 18:59:39 phil
66  * Double nmax
67  *
68  * Revision 2.11 1999/10/11 14:29:37 phil
69  * &-> &&
70  *
71  * Revision 2.10 1999/09/30 09:23:25 phil
72  * remplacement des && en &
73  *
74  * Revision 2.9 1999/09/17 15:26:25 phil
75  * correction definition de NMAX
76  *
77  * Revision 2.8 1999/06/23 12:35:00 phil
78  * *** empty log message ***
79  *
80  * Revision 2.7 1999/04/28 10:49:19 phil
81  * augmentation de nmax a 50
82  *
83  * Revision 2.6 1999/04/27 13:12:34 phil
84  * *** empty log message ***
85  *
86  * Revision 2.5 1999/04/19 14:07:02 phil
87  * *** empty log message ***
88  *
89  * Revision 2.4 1999/04/16 13:19:07 phil
90  * *** empty log message ***
91  *
92  * Revision 2.3 1999/04/16 13:17:02 phil
93  * *** empty log message ***
94  *
95  * Revision 2.2 1999/04/14 13:57:05 phil
96  * Sauvegarde des solutions deja calculees
97  *
98  * Revision 2.1 1999/04/07 14:39:23 phil
99  * esthetique
100  *
101  * Revision 2.0 1999/04/07 14:11:18 phil
102  * *** empty log message ***
103  *
104  *
105  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/PDE/solh.C,v 1.11 2016/12/05 16:18:10 j_novak Exp $
106  *
107  */
108 
109 //fichiers includes
110 #include <cstdio>
111 #include <cstdlib>
112 #include <cmath>
113 
114 #include "proto.h"
115 #include "matrice.h"
116 #include "type_parite.h"
117 
118 
119 /*
120  *
121  * Renvoie une ou 2 solutions homogenes
122  * Si base_r = R_CHEB deux solutions (x+echelle)^l dans (0, *) et
123  * 1/(x+echelle)^(l+1) dans (1, *)
124  * Si base_r = R_CHEBU 1 solution (x-1)^l+1 dans (*)
125  * Si base_r = R_CHEBP ou R_CHEBI x^l dans (*)
126  *
127  * Entree :
128  * n : nbre de points en r
129  * l : nbre quantique associe
130  * echelle : cf ci-dessus, utile que dans le cas R_CHEB
131  * base_r : base de decomposition
132  *
133  * Sortie :
134  * Tbl contenant les coefficient de la ou des solution homogenes
135  *
136  */
137 
138  //------------------------------------
139  // Routine pour les cas non prevus --
140  //------------------------------------
141 namespace Lorene {
142 Tbl _solh_pas_prevu (int n, int l, double echelle) {
143 
144  cout << " Solution homogene pas prevue ..... : "<< endl ;
145  cout << " N : " << n << endl ;
146  cout << " l : " << l << endl ;
147  cout << " echelle : " << echelle << endl ;
148  abort() ;
149  exit(-1) ;
150  Tbl res(1) ;
151  return res;
152 }
153 
154 
155  //-------------------
156  //-- R_CHEB ------
157  //-------------------
158 
159 Tbl _solh_r_cheb (int n, int l, double echelle) {
160 
161  const int nmax = 200 ; // Nombre de Tbl stockes
162  static Tbl* tab[nmax] ; // les Tbl calcules
163  static int nb_dejafait = 0 ; // nbre de Tbl calcules
164  static int l_dejafait[nmax] ;
165  static int nr_dejafait[nmax] ;
166  static double vieux_echelle = 0;
167 
168  // Si on a change l'echelle : on detruit tout :
169  if (vieux_echelle != echelle) {
170  for (int i=0 ; i<nb_dejafait ; i++) {
171  l_dejafait[i] = -1 ;
172  nr_dejafait[i] = -1 ;
173  delete tab[i] ;
174  }
175  nb_dejafait = 0 ;
176  vieux_echelle = echelle ;
177  }
178 
179  int indice = -1 ;
180 
181  // On determine si la matrice a deja ete calculee :
182  for (int conte=0 ; conte<nb_dejafait ; conte ++)
183  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
184  indice = conte ;
185 
186  // Calcul a faire :
187  if (indice == -1) {
188  if (nb_dejafait >= nmax) {
189  cout << "_solh_r_cheb : trop de Tbl" << endl ;
190  abort() ;
191  exit (-1) ;
192  }
193 
194 
195  l_dejafait[nb_dejafait] = l ;
196  nr_dejafait[nb_dejafait] = n ;
197 
198  // assert (l < n) ;
199 
200  tab[nb_dejafait] = new Tbl(2, n) ;
201  Tbl* pres = tab[nb_dejafait] ;
202  pres->set_etat_qcq() ;
203  double* coloc = new double[n] ;
204 
205  int * deg = new int[3] ;
206  deg[0] = 1 ;
207  deg[1] = 1 ;
208  deg[2] = n ;
209 
210  //Construction de la premiere solution homogene :
211  // cad celle polynomiale.
212 
213  if (l==0) {
214  pres->set(0, 0) = 1 ;
215  for (int i=1 ; i<n ; i++)
216  pres->set(0, i) = 0 ;
217  }
218  else {
219  for (int i=0 ; i<n ; i++)
220  coloc[i] = pow(echelle-cos(M_PI*i/(n-1)), double(l)) ;
221 
222  cfrcheb(deg, deg, coloc, deg, coloc) ;
223  for (int i=0 ; i<n ;i++)
224  pres->set(0, i) = coloc[i] ;
225  }
226 
227 
228  // construction de la seconde solution homogene :
229  // cad celle fractionnelle.
230  for (int i=0 ; i<n ; i++)
231  coloc[i] = 1/pow(echelle-cos(M_PI*i/(n-1)), double(l+1)) ;
232 
233  cfrcheb(deg, deg, coloc, deg, coloc) ;
234  for (int i=0 ; i<n ;i++)
235  pres->set(1, i) = coloc[i] ;
236 
237 
238  delete [] coloc ;
239  delete [] deg ;
240  indice = nb_dejafait ;
241  nb_dejafait ++ ;
242  }
243 
244  return *tab[indice] ;
245 }
246 
247 
248  //-------------------
249  //-- R_JACO02 ------
250  //-------------------
251 
252 Tbl _solh_r_jaco02 (int n, int l, double echelle) {
253 
254  const int nmax = 200 ; // Nombre de Tbl stockes
255  static Tbl* tab[nmax] ; // les Tbl calcules
256  static int nb_dejafait = 0 ; // nbre de Tbl calcules
257  static int l_dejafait[nmax] ;
258  static int nr_dejafait[nmax] ;
259  static double vieux_echelle = 0;
260 
261  // Si on a change l'echelle : on detruit tout :
262  if (vieux_echelle != echelle) {
263  for (int i=0 ; i<nb_dejafait ; i++) {
264  l_dejafait[i] = -1 ;
265  nr_dejafait[i] = -1 ;
266  delete tab[i] ;
267  }
268  nb_dejafait = 0 ;
269  vieux_echelle = echelle ;
270  }
271 
272  int indice = -1 ;
273 
274  // On determine si la matrice a deja ete calculee :
275  for (int conte=0 ; conte<nb_dejafait ; conte ++)
276  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
277  indice = conte ;
278 
279  // Calcul a faire :
280  if (indice == -1) {
281  if (nb_dejafait >= nmax) {
282  cout << "_solh_r_jaco02 : trop de Tbl" << endl ;
283  abort() ;
284  exit (-1) ;
285  }
286 
287 
288  l_dejafait[nb_dejafait] = l ;
289  nr_dejafait[nb_dejafait] = n ;
290 
291  // assert (l < n) ;
292 
293  tab[nb_dejafait] = new Tbl(2, n) ;
294  Tbl* pres = tab[nb_dejafait] ;
295  pres->set_etat_qcq() ;
296  double* coloc = new double[n] ;
297 
298  int * deg = new int[3] ;
299  deg[0] = 1 ;
300  deg[1] = 1 ;
301  deg[2] = n ;
302 
303 
304  double* zeta = pointsgausslobatto(n-1) ;
305  //Construction de la premiere solution homogene :
306  // cad celle polynomiale.
307 
308  if (l==0) {
309  pres->set(0, 0) = 1 ;
310  for (int i=1 ; i<n ; i++)
311  pres->set(0, i) = 0 ;
312  }
313  else {
314  for (int i=0 ; i<n ; i++)
315  coloc[i] = pow((echelle + zeta[i]), double(l)) ;
316 
317  cfrjaco02(deg, deg, coloc, deg, coloc) ;
318  for (int i=0 ; i<n ;i++)
319  pres->set(0, i) = coloc[i] ;
320  }
321 
322 
323  // construction de la seconde solution homogene :
324  // cad celle fractionnelle.
325  for (int i=0 ; i<n ; i++)
326  coloc[i] = 1/pow((echelle + zeta[i]), double(l+1)) ;
327 
328  cfrjaco02(deg, deg, coloc, deg, coloc) ;
329  for (int i=0 ; i<n ;i++)
330  pres->set(1, i) = coloc[i] ;
331 
332 
333  delete [] coloc ;
334  delete [] deg ;
335  indice = nb_dejafait ;
336  nb_dejafait ++ ;
337  }
338 
339  return *tab[indice] ;
340 }
341 
342 
343 
344  //-------------------
345  //-- R_CHEBP ------
346  //-------------------
347 
348 Tbl _solh_r_chebp (int n, int l, double) {
349 
350  const int nmax = 200 ; // Nombre de Tbl stockes
351  static Tbl* tab[nmax] ; // les Tbl calcules
352  static int nb_dejafait = 0 ; // nbre de Tbl calcules
353  static int l_dejafait[nmax] ;
354  static int nr_dejafait[nmax] ;
355 
356  int indice = -1 ;
357 
358  // On determine si la matrice a deja ete calculee :
359  for (int conte=0 ; conte<nb_dejafait ; conte ++)
360  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
361  indice = conte ;
362 
363  // Calcul a faire :
364  if (indice == -1) {
365  if (nb_dejafait >= nmax) {
366  cout << "_solh_r_chebp : trop de Tbl" << endl ;
367  abort() ;
368  exit (-1) ;
369  }
370 
371 
372  l_dejafait[nb_dejafait] = l ;
373  nr_dejafait[nb_dejafait] = n ;
374 
375  assert (div(l, 2).rem ==0) ;
376 
377  // assert (l < 2*n-1) ;
378 
379  tab[nb_dejafait] = new Tbl(n) ;
380  Tbl* pres = tab[nb_dejafait] ;
381  pres->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  if (l==0) {
390  pres->set(0) = 1 ;
391  for (int i=1 ; i<n ; i++)
392  pres->set(i) = 0 ;
393  }
394  else {
395  for (int i=0 ; i<n ; i++)
396  coloc[i] = pow(sin(M_PI*i/2/(n-1)), double(l)) ;
397 
398  cfrchebp(deg, deg, coloc, deg, coloc) ;
399  for (int i=0 ; i<n ;i++)
400  pres->set(i) = coloc[i] ;
401  }
402 
403  delete [] coloc ;
404  delete [] deg ;
405  indice = nb_dejafait ;
406  nb_dejafait ++ ;
407  }
408 
409  return *tab[indice] ;
410 }
411 
412 
413  //-------------------
414  //-- R_CHEBI -----
415  //-------------------
416 
417 Tbl _solh_r_chebi (int n, int l, double) {
418 
419  const int nmax = 200 ; // Nombre de Tbl stockes
420  static Tbl* tab[nmax] ; // les Tbl calcules
421  static int nb_dejafait = 0 ; // nbre de Tbl calcules
422  static int l_dejafait[nmax] ;
423  static int nr_dejafait[nmax] ;
424 
425  int indice = -1 ;
426 
427  // On determine si la matrice a deja ete calculee :
428  for (int conte=0 ; conte<nb_dejafait ; conte ++)
429  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
430  indice = conte ;
431 
432  // Calcul a faire :
433  if (indice == -1) {
434  if (nb_dejafait >= nmax) {
435  cout << "_solh_r_chebi : trop de Tbl" << endl ;
436  abort() ;
437  exit (-1) ;
438  }
439 
440 
441  l_dejafait[nb_dejafait] = l ;
442  nr_dejafait[nb_dejafait] = n ;
443 
444 
445  assert (div(l, 2).rem == 1) ;
446  // assert (l < 2*n) ;
447 
448  tab[nb_dejafait] = new Tbl(n) ;
449  Tbl* pres = tab[nb_dejafait] ;
450  pres->set_etat_qcq() ;
451  double* coloc = new double[n] ;
452 
453  int * deg = new int[3] ;
454  deg[0] = 1 ;
455  deg[1] = 1 ;
456  deg[2] = n ;
457 
458  if (l==1) {
459  pres->set(0) = 1 ;
460  for (int i=1 ; i<n ; i++)
461  pres->set(i) = 0 ;
462  }
463  else {
464  for (int i=0 ; i<n ; i++)
465  coloc[i] = pow(sin(M_PI*i/2/(n-1)), double(l)) ;
466 
467  cfrchebi(deg, deg, coloc, deg, coloc) ;
468  for (int i=0 ; i<n ;i++)
469  pres->set(i) = coloc[i] ;
470  }
471 
472  delete [] coloc ;
473  delete [] deg ;
474  indice = nb_dejafait ;
475  nb_dejafait ++ ;
476  }
477 
478  return *tab[indice] ;
479 }
480 
481 
482 
483  //-------------------
484  //-- R_CHEBU -----
485  //-------------------
486 
487 Tbl _solh_r_chebu (int n, int l, double) {
488 
489  const int nmax = 200 ; // Nombre de Tbl stockes
490  static Tbl* tab[nmax] ; // les Tbl calcules
491  static int nb_dejafait = 0 ; // nbre de Tbl calcules
492  static int l_dejafait[nmax] ;
493  static int nr_dejafait[nmax] ;
494 
495  int indice = -1 ;
496 
497  // On determine si la matrice a deja ete calculee :
498  for (int conte=0 ; conte<nb_dejafait ; conte ++)
499  if ((l_dejafait[conte] == l) && (nr_dejafait[conte] == n))
500  indice = conte ;
501 
502  // Calcul a faire :
503  if (indice == -1) {
504  if (nb_dejafait >= nmax) {
505  cout << "_solh_r_chebu : trop de Tbl" << endl ;
506  abort() ;
507  exit (-1) ;
508  }
509 
510 
511  l_dejafait[nb_dejafait] = l ;
512  nr_dejafait[nb_dejafait] = n ;
513 
514 
515  // assert (l < n-1) ;
516  tab[nb_dejafait] = new Tbl(n) ;
517  Tbl* pres = tab[nb_dejafait] ;
518  pres->set_etat_qcq() ;
519  double* coloc = new double[n] ;
520 
521  int * deg = new int[3] ;
522  deg[0] = 1 ;
523  deg[1] = 1 ;
524  deg[2] = n ;
525 
526  for (int i=0 ; i<n ; i++)
527  coloc[i] = pow(-1-cos(M_PI*i/(n-1)), double(l+1)) ;
528 
529  cfrcheb(deg, deg, coloc, deg, coloc) ;
530  for (int i=0 ; i<n ;i++)
531  pres->set(i) = coloc[i] ;
532 
533  delete [] coloc ;
534  delete [] deg ;
535  indice = nb_dejafait ;
536  nb_dejafait ++ ;
537  }
538 
539  return *tab[indice] ;
540 }
541 
542 
543 
544 
545  //-------------------
546  //-- Fonction ----
547  //-------------------
548 
549 
550 Tbl solh(int n, int l, double echelle, int base_r) {
551 
552  // Routines de derivation
553  static Tbl (*solh[MAX_BASE])(int, int, double) ;
554  static int nap = 0 ;
555 
556  // Premier appel
557  if (nap==0) {
558  nap = 1 ;
559  for (int i=0 ; i<MAX_BASE ; i++) {
560  solh[i] = _solh_pas_prevu ;
561  }
562  // Les routines existantes
563  solh[R_CHEB >> TRA_R] = _solh_r_cheb ;
564  solh[R_CHEBU >> TRA_R] = _solh_r_chebu ;
565  solh[R_CHEBP >> TRA_R] = _solh_r_chebp ;
566  solh[R_CHEBI >> TRA_R] = _solh_r_chebi ;
567  solh[R_JACO02 >> TRA_R] = _solh_r_jaco02 ;
568  }
569 
570  return solh[base_r](n, l, echelle) ;
571 }
572 }
Lorene prototypes.
Definition: app_hor.h:67
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
Definition: type_parite.h:188
Cmp cos(const Cmp &)
Cosine.
Definition: cmp_math.C:97
#define TRA_R
Translation en R, used for a bitwise shift (in hex)
Definition: type_parite.h:158
#define R_CHEBI
base de Cheb. impaire (rare) seulement
Definition: type_parite.h:170
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition: type_parite.h:168
Cmp pow(const Cmp &, int)
Power .
Definition: cmp_math.C:351
Cmp sin(const Cmp &)
Sine.
Definition: cmp_math.C:72
#define R_CHEBU
base de Chebychev ordinaire (fin), dev. en 1/r
Definition: type_parite.h:180
#define MAX_BASE
Nombre max. de bases differentes.
Definition: type_parite.h:144
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition: type_parite.h:166