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