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 $" ;
52#include "utilitaires.h"
58 int nwork_colloc = 0 ;
60 double* tab_colloc[nmax] ;
67void poly_leg (
int n,
double& poly,
double& pder,
double& polym1,
double& pderm1,
68 double& polym2,
double& pderm2,
double x) {
87 for (
int i=1 ; i<n ; i++) {
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) ;
99void legendre_collocation_points(
int nr,
double* colloc) {
102 for (
int i=0; ((i<nwork_colloc) && (index<0)); i++)
103 if (nb_colloc[i] == nr) index = i ;
107 index = nwork_colloc ;
109 cout <<
"legendre_collocation_points: " << endl ;
110 cout <<
"too many arrays!" << endl ;
113 double*& t_colloc = tab_colloc[index] ;
114 t_colloc =
new double[nr] ;
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 ;
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) ;
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 ;
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) ;
141 int borne_sup = (nr%2==0) ? nr/2 : (nr+1)/2 ;
143 for (
int j=1 ; j<borne_sup ; j++) {
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 ;
152 for (
int i=0 ; i<j ; i++)
153 sum += 1./(x-t_colloc[nr-i-1]) ;
155 double increm = -poly/(pder-sum*poly) ;
159 if ((fabs(increm) < 1.e-14) || (ite >500))
163 cout <<
"leg_ini: too many iterations..." << endl ;
166 t_colloc[nr-j-1] = x ;
167 double auxi = cs*cd-ss*sd ;
172 t_colloc[(nr-1)/2] = 0 ;
174 for (
int i=0 ; i<borne_sup ; i++)
175 t_colloc[i] = - t_colloc[nr-i-1] ;
176 nb_colloc[index] = nr ;
179 assert((index>=0)&&(index<nmax)) ;
180 for (
int i=0; i<nr; i++)
181 colloc[i] = (tab_colloc[index])[i] ;
191void get_legendre_data(
int np, Tbl*& p_Pni, Tbl*& p_wn) {
194 for (
int i=0; ((i<nwork_leg) && (index<0)); i++)
195 if (nb_leg[i] == np) index = i ;
200 cout <<
"get_legendre_data: " << endl ;
201 cout <<
"too many transformation matrices!" << endl ;
205 tab_pni[index] =
new Tbl(np, np) ;
206 Tbl& Pni = (*tab_pni[index]) ;
208 tab_wn[index] =
new Tbl(np) ;
209 Tbl& wn = (*tab_wn[index]) ;
213 coloc.set_etat_qcq() ;
214 legendre_collocation_points(np, coloc.t) ;
216 for (
int i=0; i<np; i++) {
218 Pni.set(1, i) = coloc(i) ;
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) ;
227 for (
int j=0; j<np; j++)
228 wn.set(j) = 2./( double(np0)*double(np) * Pni(np0,j) * Pni(np0, j) ) ;
233 assert((index>=0)&&(index<nmax)) ;
234 p_Pni = tab_pni[index] ;
235 p_wn = tab_wn[index] ;
Cmp sin(const Cmp &)
Sine.
Cmp cos(const Cmp &)
Cosine.