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 
23 char 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 
54 namespace Lorene {
55 
56 namespace {
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 
67 void 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 /************************************************************************/
99 void 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 
191 void 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 }
Lorene
Lorene prototypes.
Definition: app_hor.h:64
Lorene::cos
Cmp cos(const Cmp &)
Cosine.
Definition: cmp_math.C:94
Lorene::sin
Cmp sin(const Cmp &)
Sine.
Definition: cmp_math.C:69