LORENE
legendre.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_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/legendre.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)) / (2m-1)!!
28 * aux points
29 * theta_j = pi/2 j/(nt-1) 0 <= j <= nt-1
30 * qui echantillonnent uniformement l'intervalle [0, pi/2].
31 *
32 *
33 * Entree:
34 * -------
35 * int m : ordre de la fonction de Legendre associee P_l^m
36 * int nt : nombre de points en theta
37 *
38 * Sortie (valeur de retour) :
39 * -------------------------
40 * double* legendre : ensemble des (nt-m)*nt valeurs
41 * P_l^m(cos(theta))/(2m-1)!!
42 * stokees comme suit:
43 *
44 * legendre[nt* (l-m) + j] = P_l^m( cos(theta_j) ) / (2m-1)!!
45 *
46 * avec m <= l <= nt-1.
47 *
48 * NB: Cette routine effectue le calcul a chaque appel et ne renvoie pas
49 * un pointeur sur des valeurs precedemment calculees.
50 */
51
52
53/*
54 * $Id: legendre.C,v 1.6 2014/10/13 08:53:13 j_novak Exp $
55 * $Log: legendre.C,v $
56 * Revision 1.6 2014/10/13 08:53:13 j_novak
57 * Lorene classes and functions now belong to the namespace Lorene.
58 *
59 * Revision 1.5 2014/10/06 15:16:02 j_novak
60 * Modified #include directives to use c++ syntax.
61 *
62 * Revision 1.4 2005/02/18 13:14:13 j_novak
63 * Changing of malloc/free to new/delete + suppression of some unused variables
64 * (trying to avoid compilation warnings).
65 *
66 * Revision 1.3 2003/01/31 10:31:24 e_gourgoulhon
67 * Suppressed the directive #include <malloc.h> for malloc is defined
68 * in <stdlib.h>
69 *
70 * Revision 1.2 2002/10/16 14:36:54 j_novak
71 * Reorganization of #include instructions of standard C++, in order to
72 * use experimental version 3 of gcc.
73 *
74 * Revision 1.1.1.1 2001/11/20 15:19:28 e_gourgoulhon
75 * LORENE
76 *
77 * Revision 2.0 1999/02/22 15:37:13 hyc
78 * *** empty log message ***
79 *
80 *
81 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/legendre.C,v 1.6 2014/10/13 08:53:13 j_novak Exp $
82 *
83 */
84
85// headers du C
86#include <cstdlib>
87#include <cassert>
88#include <cmath>
89
90#include "headcpp.h"
91
92namespace Lorene {
93//******************************************************************************
94
95double* legendre(int m, int nt) {
96
97int i, j, l ;
98
99 int lmax = nt - 1 ;
100 assert(m >= 0) ;
101 assert(m <= lmax) ;
102
103 double dt = M_PI / double(2*(nt-1)) ;
104
105// Allocation memoire pour le tableau resultat
106//--------------------------------------------
107
108 double* resu = new double[(lmax-m+1)*nt] ; //(double *)(malloc( (lmax-m+1)*nt * sizeof(double) )) ;
109
110 // Tableau de travail
111 double* cost = new double[nt] ; //(double*)( malloc( nt*sizeof(double) ) ) ;
112
113//-----------------------
114// 1/ Calcul de P_m^m
115//-----------------------
116
117 if (m==0) {
118 for (j=0; j<nt; j++) {
119 resu[j] = 1. ; // P_0^0(x) = 1.
120 }
121 }
122 else {
123
124//... P_m^m(x) = (-1)^m (1-x^2)^{m/2} <--- cette formule donne un P_m^m
125// plus petit par un facteur
126// (2m-1)!! que celui de la litterature
127
128 for (j=0; j<nt; j++) {
129 double y = 1. ;
130 double s = sin(j*dt) ;
131 for (i=1 ; i<2*m; i+=2) {
132 y *= - s ;
133// NB: Pour obtenir le P_m^m de la litterature, il faudrait remplacer la ligne
134// ci-dessus par : y *= - i*s ;
135 }
136 resu[j] = y ;
137//## resu[j] = pow(-s, double(m)) ;
138 }
139 } // fin du cas m != 0
140
141 if (lmax==m) {
142 delete [] cost ;
143 return resu ;
144 }
145 else {
146
147//-----------------------
148// 2/ Calcul de P_{m+1}^m
149//-----------------------
150
151//... Calcul des cos( theta_j ) :
152 for (j=0; j<nt; j++) {
153 cost[j] = cos(j*dt) ;
154 }
155
156 for (j=0; j<nt; j++) {
157 resu[nt+j] = cost[j] * (2.*m+1) * resu[j] ;
158 }
159
160//-----------------------
161// 3/ Calcul de P_l^m pour m+2 <= l <= lmax
162//-----------------------
163
164 for (l=m+2; l < lmax+1 ; l++) {
165 int i_l = nt*(l-m) ;
166 int i_lm1 = nt*(l-1-m) ;
167 int i_lm2 = nt*(l-2-m) ;
168 int a = 2*l - 1 ;
169 int b = l + m - 1 ;
170 int c = l - m ;
171
172 for (j=0; j<nt; j++) {
173 resu[i_l+j] = ( cost[j] * a * resu[i_lm1+j]
174 - b * resu[i_lm2+j] ) / c ;
175 }
176 }
177
178 delete [] cost ; //free (cost) ;
179 return resu ;
180
181 } // fin du cas lmax > m
182
183}
184
185
186
187}
Cmp sin(const Cmp &)
Sine.
Definition cmp_math.C:69
Cmp cos(const Cmp &)
Cosine.
Definition cmp_math.C:94
Lorene prototypes.
Definition app_hor.h:64