21 #include "headcpp.hpp"
23 #define NBR_LEG_MAX 20
26 void legendre (
int n,
double& poly,
double& pder,
double& polym1,
double& pderm1,
27 double& polym2,
double& pderm2,
double x) {
46 for (
int i=1 ; i<n ; i++) {
51 poly = ((2*i+1)*x*polym1 - i*polym2)/(i+1) ;
52 pder = ((2*i+1)*polym1+(2*i+1)*x*pderm1-i*pderm2)/(i+1) ;
57 double leg (
int n,
double x) {
59 double p, dp, p1, dp1, p2, dp2 ;
60 legendre (n, p, dp, p1, dp1, p2, dp2, x) ;
64 void gauss_lobato_legendre (
int n, Array<double>& coloc, Array<double>& weight,
65 double prec = PRECISION,
int itemax = 500) {
67 assert (coloc.get_ndim()==1) ;
68 assert (coloc.get_size(0)==n) ;
69 assert (weight.get_ndim()==1) ;
70 assert (weight.get_size(0)==n) ;
74 double p_plus, dp_plus, p_plus_m1, dp_plus_m1, p_plus_m2, dp_plus_m2 ;
75 double p_moins, dp_moins, p_moins_m1, dp_moins_m1, p_moins_m2, dp_moins_m2 ;
76 double p, dp, p_m1, dp_m1, p_m2, dp_m2 ;
78 legendre (n, p_plus, dp_plus, p_plus_m1, dp_plus_m1, p_plus_m2, dp_plus_m2, x_plus) ;
79 legendre (n, p_moins, dp_moins, p_moins_m1, dp_moins_m1, p_moins_m2, dp_moins_m2, x_moins) ;
81 double det = p_plus_m1*p_moins_m2 - p_moins_m1*p_plus_m2 ;
82 double r_plus = -p_plus ;
83 double r_moins = -p_moins ;
84 double a = (r_plus*p_moins_m2 - r_moins*p_plus_m2)/det ;
85 double b = (r_moins*p_plus_m1 - r_plus*p_moins_m1)/det ;
88 double dth = M_PI/(2*n+1) ;
89 double cd = cos (2*dth) ;
90 double sd = sin (2*dth) ;
91 double cs = cos(dth) ;
92 double ss = sin(dth) ;
94 int borne_sup = (n%2==0) ?
97 for (
int j=1 ; j<borne_sup ; j++) {
102 legendre (n, p, dp, p_m1, dp_m1, p_m2, dp_m2, x) ;
103 double poly = p + a*p_m1 + b*p_m2 ;
104 double pder = dp + a * dp_m1 + b*dp_m2 ;
106 for (
int i=0 ; i<j ; i++)
107 sum += 1./(x-coloc(n-i-1)) ;
109 double increm = -poly/(pder-sum*poly) ;
113 if ((fabs(increm) < prec) || (ite >itemax))
117 cout <<
"Too many iterations..." << endl ;
120 coloc.set(n-j-1) = x ;
121 double auxi = cs*cd-ss*sd ;
126 coloc.set((n-1)/2) = 0 ;
128 for (
int i=0 ; i<borne_sup ; i++)
129 coloc.set(i) = - coloc(n-i-1) ;
131 for (
int i=0 ; i<n ; i++) {
132 legendre (n-1, p, dp, p_m1, dp_m1, p_m2, dp_m2, coloc(i)) ;
133 weight.set(i) = 2./(n-1)/n/p/p ;
138 double get_things_leg (
int indice,
int n,
int thing) {
140 static bool first = true ;
141 static int nbr_done = 0 ;
142 static int n_done[NBR_LEG_MAX] ;
143 static Array<double>** coloc =
new Array<double>* [NBR_LEG_MAX] ;
144 static Array<double>** weight =
new Array<double>* [NBR_LEG_MAX] ;
145 static Array<double>** gamma =
new Array<double>* [NBR_LEG_MAX] ;
149 for (
int i=0 ; i<NBR_LEG_MAX ; i++) {
160 for (
int i=0 ; i<nbr_done ; i++)
166 if (nbr_done>=NBR_LEG_MAX) {
167 cout <<
"Too many different number of points for Legendre..." << endl ;
171 Array<double> col (n) ;
172 Array<double> poids(n) ;
174 gauss_lobato_legendre (n, col, poids) ;
175 coloc[nbr_done] =
new Array<double>(col) ;
176 weight[nbr_done] =
new Array<double>(poids) ;
177 gamma[nbr_done] =
new Array<double>(n) ;
178 (*gamma[nbr_done]) = 0 ;
179 for (
int i=0 ; i<n ; i++)
180 for (
int j=0 ; j<n ; j++)
181 gamma[nbr_done]->set(i) += pow(leg(i,(*coloc[nbr_done])(j)),2) *
182 (*weight[nbr_done])(j) ;
183 n_done[nbr_done] = n ;
185 index = nbr_done - 1 ;
190 return (*coloc[index])(indice) ;
193 return (*weight[index])(indice) ;
196 return (*gamma[index])(indice) ;
199 cout <<
"Strange error in get_things_legendre..." << endl ;
205 double coloc_leg (
int i,
int nbr) {
206 assert ((i>=0) && (i<nbr)) ;
207 return get_things_leg(i, nbr, 0) ;
210 double weight_leg (
int i,
int nbr) {
211 assert ((i>=0) && (i<nbr)) ;
212 return get_things_leg(i, nbr, 1) ;
215 double gamma_leg (
int i,
int nbr) {
216 assert ((i>=0) && (i<nbr)) ;
217 return get_things_leg(i, nbr, 2) ;
222 double get_things_leg_parity (
int indice,
int n,
int thing) {
224 static bool first = true ;
225 static int nbr_done = 0 ;
226 static int n_done[NBR_LEG_MAX] ;
227 static Array<double>** coloc =
new Array<double>* [NBR_LEG_MAX] ;
228 static Array<double>** weight =
new Array<double>* [NBR_LEG_MAX] ;
232 for (
int i=0 ; i<NBR_LEG_MAX ; i++) {
242 for (
int i=0 ; i<nbr_done ; i++)
248 if (nbr_done>=NBR_LEG_MAX) {
249 cout <<
"Too many number of points for Legendre..." << endl ;
253 Array<double> col (2*n-1) ;
254 Array<double> poids(2*n-1) ;
256 gauss_lobato_legendre (2*n-1, col, poids) ;
257 coloc[nbr_done] =
new Array<double>(n) ;
258 for (
int i=0 ; i<n ; i++)
259 coloc[nbr_done]->set(i) = col(i+n-1) ;
260 weight[nbr_done] =
new Array<double>(n) ;
261 weight[nbr_done]->set(0) = poids(n-1)/2. ;
262 for (
int i=1 ; i<n ; i++)
263 weight[nbr_done]->set(i) = poids(i+n-1) ;
264 n_done[nbr_done] = n ;
266 index = nbr_done - 1 ;
271 return (*coloc[index])(indice) ;
274 return (*weight[index])(indice) ;
277 cout <<
"Strange error in get_things_legendre..." << endl ;
283 double coloc_leg_parity (
int i,
int nbr) {
284 assert ((i>=0) && (i<nbr)) ;
285 return get_things_leg_parity(i, nbr, 0) ;
288 double weight_leg_parity (
int i,
int nbr) {
289 assert ((i>=0) && (i<nbr)) ;
290 return get_things_leg_parity(i, nbr, 1) ;
294 double gamma_leg_even (
int nn,
int nbr) {
295 assert ((nn>=0) && (nn<nbr)) ;
296 static bool first = true ;
297 static int nbr_done = 0 ;
298 static int n_done[NBR_LEG_MAX] ;
299 static Array<double>** gamma =
new Array<double>* [NBR_LEG_MAX] ;
302 for (
int i=0 ; i<NBR_LEG_MAX ; i++) {
311 for (
int i=0 ; i<nbr_done ; i++)
312 if (n_done[i] == nbr)
317 if (nbr_done>=NBR_LEG_MAX) {
318 cout <<
"Too many number of points for Legendre..." << endl ;
321 gamma[nbr_done] =
new Array<double> (nbr) ;
322 *gamma[nbr_done] = 0 ;
323 for (
int i=0 ; i<nbr ; i++)
324 for (
int j=0 ; j<nbr ; j++)
325 gamma[nbr_done]->set(i) += pow(leg(i*2,coloc_leg_parity(j,nbr)),2) * weight_leg_parity(j, nbr) ;
326 n_done[nbr_done] = nbr ;
328 index = nbr_done - 1 ;
331 return (*gamma[index])(nn) ;
335 double gamma_leg_odd (
int nn,
int nbr) {
336 assert ((nn>=0) && (nn<nbr)) ;
337 static bool first = true ;
338 static int nbr_done = 0 ;
339 static int n_done[NBR_LEG_MAX] ;
340 static Array<double>** gamma =
new Array<double>* [NBR_LEG_MAX] ;
343 for (
int i=0 ; i<NBR_LEG_MAX ; i++) {
352 for (
int i=0 ; i<nbr_done ; i++)
353 if (n_done[i] == nbr)
358 if (nbr_done>=NBR_LEG_MAX) {
359 cout <<
"Too many number of points for Legendre..." << endl ;
362 gamma[nbr_done] =
new Array<double> (nbr) ;
363 *gamma[nbr_done] = 0 ;
364 for (
int i=0 ; i<nbr ; i++)
365 for (
int j=0 ; j<nbr ; j++)
366 gamma[nbr_done]->set(i) += pow(leg(i*2+1,coloc_leg_parity(j,nbr)),2) * weight_leg_parity(j, nbr) ;
367 n_done[nbr_done] = nbr ;
369 index = nbr_done - 1 ;
372 return (*gamma[index])(nn) ;