LORENE
legendre_norm.C
1 /*
2  * Copyright (c) 1999-2001 Eric Gourgoulhon
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 char legendre_norm_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/legendre_norm.C,v 1.6 2014/10/13 08:53:13 j_novak Exp $" ;
24 
25 /*
26  * Calcule les valeurs des fonctions de Legendre associees
27  * P_l^m(cos(theta)) normalisees de facon a ce que
28  *
29  * int_0^pi [ P_l^m(cos(theta)) ]^2 sin(theta) dtheta = 1
30  *
31  * NB: Cette normalisation est differente de celle de la litterature
32  *
33  * Le calcul est effectue aux 2*nt-1 points
34  * theta_j = pi/2 j/(2*nt-2) 0 <= j <= 2*nt-2
35  * qui echantillonnent uniformement l'intervalle [0, pi/2].
36  *
37  *
38  * Entree:
39  * -------
40  * int m : ordre de la fonction de Legendre associee P_l^m
41  * int nt : nombre de points en theta
42  *
43  * Sortie (valeur de retour) :
44  * -------------------------
45  * double* legendre_norm : ensemble des (2*nt-1-m)*(2*nt-1) valeurs
46  * P_l^m(cos(theta))
47  * stokees comme suit:
48  *
49  * legendre_norm[(2*nt-1)* (l-m) + j] = P_l^m( cos(theta_j) )
50  *
51  * avec m <= l <= 2*nt-2.
52  *
53  * NB: Cette routine effectue le calcul a chaque appel et ne renvoie pas
54  * un pointeur sur des valeurs precedemment calculees.
55  */
56 
57 
58 /*
59  * $Id: legendre_norm.C,v 1.6 2014/10/13 08:53:13 j_novak Exp $
60  * $Log: legendre_norm.C,v $
61  * Revision 1.6 2014/10/13 08:53:13 j_novak
62  * Lorene classes and functions now belong to the namespace Lorene.
63  *
64  * Revision 1.5 2014/10/06 15:16:02 j_novak
65  * Modified #include directives to use c++ syntax.
66  *
67  * Revision 1.4 2005/02/18 13:14:13 j_novak
68  * Changing of malloc/free to new/delete + suppression of some unused variables
69  * (trying to avoid compilation warnings).
70  *
71  * Revision 1.3 2003/01/31 10:31:24 e_gourgoulhon
72  * Suppressed the directive #include <malloc.h> for malloc is defined
73  * in <stdlib.h>
74  *
75  * Revision 1.2 2002/10/16 14:36:54 j_novak
76  * Reorganization of #include instructions of standard C++, in order to
77  * use experimental version 3 of gcc.
78  *
79  * Revision 1.1.1.1 2001/11/20 15:19:29 e_gourgoulhon
80  * LORENE
81  *
82  * Revision 2.0 1999/02/22 15:37:00 hyc
83  * *** empty log message ***
84  *
85  *
86  * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/legendre_norm.C,v 1.6 2014/10/13 08:53:13 j_novak Exp $
87  *
88  */
89 
90 // headers du C
91 #include <cstdlib>
92 #include <cassert>
93 #include <cmath>
94 
95 // Prototypage
96 #include "headcpp.h"
97 #include "proto.h"
98 
99 namespace Lorene {
100 //******************************************************************************
101 
102 double* legendre_norm(int m, int nt) {
103 
104 int l, j ;
105 
106  int lmax = 2*nt - 2 ;
107 
108 // Sur-echantillonnage pour calculer les carres sans aliasing:
109  int nt2 = 2*nt - 1 ;
110  int nt2m1 = nt2 - 1 ;
111 
112  int deg[3] ;
113  deg[0] = 1 ;
114  deg[1] = 1 ;
115  deg[2] = nt2 ;
116 
117 // Tableau de travail
118  double* yy = new double[nt2] ; //(double*)( malloc( nt2*sizeof(double) ) ) ;
119 
120 // Recherche des fonctions de legendre associees non normalisees
121 // -------------------------------------------------------------
122 // (NB: elles different de celles de la litterature par un facteur (2m-1)!!) :
123 
124  double* leg = legendre(m, nt2) ;
125 
126 // Normalisation
127 // -------------
128  for (l=m; l<lmax+1; l++) {
129 
130  int ml = (m+1)*(l+1) ;
131 
132 // On divise les fonctions de Legendre par (m+1)*(l+1)
133 // pour obtenir des nombres pas trop grands:
134 
135  for (j=0; j<nt2; j++) {
136  leg[nt2*(l-m)+j] /= ml ;
137  }
138 
139 // Carre :
140  for (j=0; j<nt2; j++) {
141  double w = leg[nt2*(l-m)+j] ;
142  yy[nt2m1-j] = w * w ; // le rangement est celui qui convient
143  // a cfrchebp
144  }
145 
146 // Developpement en polynomes de Tchebyshev pairs (x=cos(theta)) :
147 
148  cfrchebp(deg, deg, yy, deg, yy) ;
149 
150 // Integrale sur [0,Pi] = 2 fois l'integrale sur [0,1] pour x = cos(theta) :
151  double integ = 2.*int1d_chebp(nt2, yy) ;
152 
153 // Facteur de normalisation
154  double fact = 1. / sqrt(integ) ;
155 
156 /* Test: Comparaison avec le resultat analytique
157  *
158  * double fact_test = ml* factorielle2(2*m-1) * sqrt( double(2*l+1)/2.
159  * * factorielle(l-m) / factorielle(l+m) ) ;
160  * double diff = (fact - fact_test) / fact_test ;
161  *
162  * cout << "m, l : "<< m << " " << l << " : " << fact << " " << fact_test
163  * << " " << diff << endl ;
164  */
165 
166  for (j=0; j<nt2; j++) {
167  leg[nt2*(l-m)+j] *= fact ;
168  }
169 
170  } // fin de la boucle sur l
171 
172 // Liberation espace memoire :
173  delete [] yy ;
174 
175  return leg ;
176 
177 }
178 
179 
180 
181 }
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::sqrt
Cmp sqrt(const Cmp &)
Square root.
Definition: cmp_math.C:220