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
23char 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
99namespace Lorene {
100//******************************************************************************
101
102double* legendre_norm(int m, int nt) {
103
104int 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}
Cmp sqrt(const Cmp &)
Square root.
Definition cmp_math.C:220
Lorene prototypes.
Definition app_hor.h:64