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 //------------------------------------
141namespace Lorene {
142Tbl _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
159Tbl _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
252Tbl _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
348Tbl _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
417Tbl _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
487Tbl _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
550Tbl 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::Tbl
Basic array class.
Definition tbl.h:161
Lorene::Tbl::set_etat_qcq
void set_etat_qcq()
Sets the logical state to ETATQCQ (ordinary state).
Definition tbl.C:364
Lorene::sin
Cmp sin(const Cmp &)
Sine.
Definition cmp_math.C:72
Lorene::pow
Cmp pow(const Cmp &, int)
Power .
Definition cmp_math.C:351
Lorene::cos
Cmp cos(const Cmp &)
Cosine.
Definition cmp_math.C:97
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
R_JACO02
#define R_JACO02
base de Jacobi(0,2) ordinaire (finjac)
Definition type_parite.h:188
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
R_CHEB
#define R_CHEB
base de Chebychev ordinaire (fin)
Definition type_parite.h:166
R_CHEBP
#define R_CHEBP
base de Cheb. paire (rare) seulement
Definition type_parite.h:168
Lorene
Lorene prototypes.
Definition app_hor.h:67