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