LORENE
leg_ini.C
1/*
2 * Copyright (c) 2013 Jerome Novak
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 leg_ini_C[] = "$Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/leg_ini.C,v 1.3 2014/10/13 08:53:13 j_novak Exp $" ;
24
25
26/*
27 * $Id: leg_ini.C,v 1.3 2014/10/13 08:53:13 j_novak Exp $
28 * $Log: leg_ini.C,v $
29 * Revision 1.3 2014/10/13 08:53:13 j_novak
30 * Lorene classes and functions now belong to the namespace Lorene.
31 *
32 * Revision 1.2 2013/06/06 15:31:33 j_novak
33 * Functions to compute Legendre coefficients (not fully tested yet).
34 *
35 * Revision 1.1 2013/06/05 15:08:13 j_novak
36 * Initial revision. Not ready yet...
37 *
38 *
39 *
40 * $Header: /cvsroot/Lorene/C++/Source/Non_class_members/Coef/leg_ini.C,v 1.3 2014/10/13 08:53:13 j_novak Exp $
41 *
42 */
43
44
45// headers du C
46#include <cstdlib>
47#include <cassert>
48#include <cmath>
49
50//Lorene prototypes
51#include "tbl.h"
52#include "utilitaires.h"
53
54namespace Lorene {
55
56namespace {
57 const int nmax = 50 ; //Maximal number of Legendre transforms sizes
58 int nwork_colloc = 0 ;
59 int nwork_leg = 0 ;
60 double* tab_colloc[nmax] ;
61 int nb_colloc[nmax] ;
62 Tbl* tab_pni[nmax] ;
63 Tbl* tab_wn[nmax] ;
64 int nb_leg[nmax] ;
65}
66
67void poly_leg (int n, double& poly, double& pder, double& polym1, double& pderm1,
68 double& polym2, double& pderm2, double x) {
69
70
71 if (n==0) {
72 poly = 1 ;
73 pder = 0 ;
74 }
75 else
76 if (n==1) {
77 polym1 = 1 ;
78 pderm1 = 0 ;
79 poly = x ;
80 pder = 1 ;
81 }
82 else {
83 polym1 = 1 ;
84 pderm1 = 0 ;
85 poly = x ;
86 pder = 1 ;
87 for (int i=1 ; i<n ; i++) {
88 polym2 = polym1 ;
89 pderm2 = pderm1 ;
90 polym1 = poly ;
91 pderm1 = pder ;
92 poly = ((2*i+1)*x*polym1 - i*polym2)/(i+1) ;
93 pder = ((2*i+1)*polym1+(2*i+1)*x*pderm1-i*pderm2)/(i+1) ;
94 }
95 }
96}
97
98/************************************************************************/
99void legendre_collocation_points(int nr, double* colloc) {
100
101 int index = -1 ;
102 for (int i=0; ((i<nwork_colloc) && (index<0)); i++)
103 if (nb_colloc[i] == nr) index = i ; //Have the collocation points already been
104 // computed?
105
106 if (index <0) { //New array needed
107 index = nwork_colloc ;
108 if (index >= nmax) {
109 cout << "legendre_collocation_points: " << endl ;
110 cout << "too many arrays!" << endl ;
111 abort() ;
112 }
113 double*& t_colloc = tab_colloc[index] ;
114 t_colloc = new double[nr] ;
115 int nr0 = nr - 1 ;
116
117 double x_plus = 1 ;
118 double x_moins = -1 ;
119 double p_plus, dp_plus, p_plus_m1, dp_plus_m1, p_plus_m2, dp_plus_m2 ;
120 double p_moins, dp_moins, p_moins_m1, dp_moins_m1, p_moins_m2, dp_moins_m2 ;
121 double p, dp, p_m1, dp_m1, p_m2, dp_m2 ;
122
123 poly_leg (nr, p_plus, dp_plus, p_plus_m1, dp_plus_m1, p_plus_m2,
124 dp_plus_m2, x_plus) ;
125 poly_leg (nr, p_moins, dp_moins, p_moins_m1, dp_moins_m1, p_moins_m2,
126 dp_moins_m2, x_moins) ;
127
128 double det = p_plus_m1*p_moins_m2 - p_moins_m1*p_plus_m2 ;
129 double r_plus = -p_plus ;
130 double r_moins = -p_moins ;
131 double a = (r_plus*p_moins_m2 - r_moins*p_plus_m2)/det ;
132 double b = (r_moins*p_plus_m1 - r_plus*p_moins_m1)/det ;
133
134 t_colloc[nr0] = 1 ;
135 double dth = M_PI/(2*nr+1) ;
136 double cd = cos (2*dth) ;
137 double sd = sin (2*dth) ;
138 double cs = cos(dth) ;
139 double ss = sin(dth) ;
140
141 int borne_sup = (nr%2==0) ? nr/2 : (nr+1)/2 ;
142
143 for (int j=1 ; j<borne_sup ; j++) {
144 double x = cs ;
145 bool loop = true ;
146 int ite = 0 ;
147 while (loop) {
148 poly_leg (nr, p, dp, p_m1, dp_m1, p_m2, dp_m2, x) ;
149 double poly = p + a*p_m1 + b*p_m2 ;
150 double pder = dp + a * dp_m1 + b*dp_m2 ;
151 double sum = 0 ;
152 for (int i=0 ; i<j ; i++)
153 sum += 1./(x-t_colloc[nr-i-1]) ;
154
155 double increm = -poly/(pder-sum*poly) ;
156
157 x += increm ;
158 ite ++ ;
159 if ((fabs(increm) < 1.e-14) || (ite >500))
160 loop = false ;
161 }
162 if (ite > 500) {
163 cout << "leg_ini: too many iterations..." << endl ;
164 abort() ;
165 }
166 t_colloc[nr-j-1] = x ;
167 double auxi = cs*cd-ss*sd ;
168 ss = cs*sd+ss*cd ;
169 cs = auxi ;
170 }
171 if (nr%2==1)
172 t_colloc[(nr-1)/2] = 0 ;
173 // Copy of the symetric ones :
174 for (int i=0 ; i<borne_sup ; i++)
175 t_colloc[i] = - t_colloc[nr-i-1] ;
176 nb_colloc[index] = nr ;
177 nwork_colloc++ ;
178 }
179 assert((index>=0)&&(index<nmax)) ;
180 for (int i=0; i<nr; i++)
181 colloc[i] = (tab_colloc[index])[i] ;
182
183 return ;
184
185}
186
187
188
189/************************************************************************/
190
191void get_legendre_data(int np, Tbl*& p_Pni, Tbl*& p_wn) {
192
193 int index = -1 ;
194 for (int i=0; ((i<nwork_leg) && (index<0)); i++)
195 if (nb_leg[i] == np) index = i ; //Has the plan already been estimated?
196
197 if (index <0) { //New plan needed
198 index = nwork_leg ;
199 if (index >= nmax) {
200 cout << "get_legendre_data: " << endl ;
201 cout << "too many transformation matrices!" << endl ;
202 abort() ;
203 }
204 int np0 = np - 1 ;
205 tab_pni[index] = new Tbl(np, np) ;
206 Tbl& Pni = (*tab_pni[index]) ;
207 Pni.set_etat_qcq() ;
208 tab_wn[index] = new Tbl(np) ;
209 Tbl& wn = (*tab_wn[index]) ;
210 wn.set_etat_qcq() ;
211
212 Tbl coloc(np) ;
213 coloc.set_etat_qcq() ;
214 legendre_collocation_points(np, coloc.t) ;
215
216 for (int i=0; i<np; i++) {
217 Pni.set(0, i) = 1 ;
218 Pni.set(1, i) = coloc(i) ;
219 }
220 for (int n=2; n<np; n++) {
221 for (int i=0; i<np; i++) {
222 Pni.set(n,i) = (2. - 1./double(n))*coloc(i)*Pni(n-1, i)
223 - (1. - 1./double(n))*Pni(n-2, i) ;
224 }
225 }
226
227 for (int j=0; j<np; j++)
228 wn.set(j) = 2./( double(np0)*double(np) * Pni(np0,j) * Pni(np0, j) ) ;
229
230 nb_leg[index] = np ;
231 nwork_leg++ ;
232 }
233 assert((index>=0)&&(index<nmax)) ;
234 p_Pni = tab_pni[index] ;
235 p_wn = tab_wn[index] ;
236
237 return ;
238}
239
240}
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