KADATH
coef_1d.cpp
1 /*
2  Copyright 2017 Philippe Grandclement
3  2018 Ludwig Jens Papenfort
4 
5  This file is part of Kadath.
6 
7  Kadath is free software: you can redistribute it and/or modify
8  it under the terms of the GNU General Public License as published by
9  the Free Software Foundation, either version 3 of the License, or
10  (at your option) any later version.
11 
12  Kadath is distributed in the hope that it will be useful,
13  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  GNU General Public License for more details.
16 
17  You should have received a copy of the GNU General Public License
18  along with Kadath. If not, see <http://www.gnu.org/licenses/>.
19 */
20 
21 #include "base_spectral.hpp"
22 #include "base_fftw.hpp"
23 #include "headcpp.hpp"
24 #include "array.hpp"
25 #include <fftw3.h>
26 #include "math.h"
27 #include <unordered_map>
28 #include <tuple>
29 
30 namespace Kadath {
31 // fftw3 computes optimal algorithm once
32 // then we store it both in a map for later use
33 std::unordered_map<int, fftw_precomp_t> fftw_precomp_map;
34 
35 // get or create buffer and plan
36 fftw_precomp_t& coef_1d_fftw(int n) {
37  auto precomp_it = fftw_precomp_map.find(n);
38  if(precomp_it == fftw_precomp_map.end())
39  precomp_it = fftw_precomp_map.emplace(std::piecewise_construct,
40  std::forward_as_tuple(n),
41  std::forward_as_tuple(n,FFTW_R2HC)).first;
42  return (*precomp_it).second;
43 }
44 
45 void coef_i_1d (int, Array<double>&) ;
46 
47 void coef_1d_pasprevu (Array<double>&) {
48  cout << "Coef_1d not implemented." << endl ;
49  abort() ;
50 }
51 
52 void coef_1d_cheb (Array<double>& tab) {
53 
54  assert (tab.get_ndim()==1) ;
55  int nr = tab.get_size(0) ;
56 
57  double fmoins0 = 0.5 * ( tab(0) - tab(nr-1));
58 
59  auto & fftw_data = coef_1d_fftw(nr-1);
60  double* tab_auxi = fftw_data.buffer;
61  fftw_plan p = fftw_data.plan;
62 
63  for (int i=1 ; i<(nr-1)/2 ; i++) {
64  double fp = 0.5*(tab(i)+tab(nr-1-i)) ;
65  double fms = 0.5*(tab(i)-tab(nr-1-i))*sin(M_PI*i/(nr-1)) ;
66  tab_auxi[i] = fp + fms ;
67  tab_auxi[nr-1-i] = fp-fms ;
68  }
69  tab_auxi[0] = 0.5*(tab(0) + tab(nr-1)) ;
70  tab_auxi[(nr-1)/2] = tab((nr-1)/2) ;
71 
72  fftw_data.execute() ;
73 
74  // Coefficient pairs :
75  tab.set(0) = tab_auxi[0] / (nr-1) ;
76  for (int i=2; i<nr-1; i += 2)
77  tab.set(i) = 2*tab_auxi[(i/2)]/(nr-1) ;
78  tab.set(nr-1) = tab_auxi[(nr-1)/2] / (nr-1) ;
79 
80  // Coefficients impairs :
81  tab.set(1) = 0 ;
82  double som = 0;
83  for (int i = 3; i < nr; i += 2 ) {
84  tab.set(i) = tab(i-2) -4 * tab_auxi[nr-1-(i-1)/2]/(nr-1) ;
85  som += tab(i) ;
86  }
87 
88 // 2. Calcul de c_1 :
89  double c1 = - ( fmoins0 + som ) / ((nr-1)/2) ;
90 
91 // 3. Coef. c_k avec k impair:
92  tab.set(1) = c1 ;
93  for (int i = 3; i < nr; i += 2 )
94  tab.set(i) += c1 ;
95 }
96 
97 void coef_1d_cheb_even (Array<double>& tab) {
98  assert (tab.get_ndim()==1) ;
99  int nr = tab.get_size(0) ;
100 
101  double fmoins0 = - 0.5 * ( tab(0) - tab(nr-1));
102 
103  auto & fftw_data = coef_1d_fftw(nr-1);
104  double* tab_auxi = fftw_data.buffer;
105  fftw_plan p = fftw_data.plan;
106 
107  for (int i=1 ; i<(nr-1)/2 ; i++) {
108  double fp = 0.5*(tab(i)+tab(nr-1-i)) ;
109  double fms = 0.5*(-tab(i)+tab(nr-1-i))*sin(M_PI*i/(nr-1)) ;
110  tab_auxi[i] = fp + fms ;
111  tab_auxi[nr-1-i] = fp-fms ;
112  }
113  tab_auxi[0] = 0.5*(tab(0) + tab(nr-1)) ;
114  tab_auxi[(nr-1)/2] = tab((nr-1)/2) ;
115 
116  fftw_data.execute() ;
117 
118  // Coefficient pairs :
119  tab.set(0) = tab_auxi[0] / (nr-1) ;
120  for (int i=2; i<nr-1; i += 2)
121  tab.set(i) = 2*tab_auxi[(i/2)]/(nr-1) ;
122  tab.set(nr-1) = tab_auxi[(nr-1)/2] / (nr-1) ;
123 
124  // Coefficients impairs :
125  tab.set(1) = 0 ;
126  double som = 0;
127  for (int i = 3; i < nr; i += 2 ) {
128  tab.set(i) = tab(i-2) +4 * tab_auxi[nr-1-i/2]/(nr-1) ;
129  som += tab(i) ;
130  }
131 
132 // 2. Calcul de c_1 :
133  double c1 = ( fmoins0 - som ) / ((nr-1)/2) ;
134 
135 // 3. Coef. c_k avec k impair:
136  tab.set(1) = c1 ;
137  for (int i = 3; i < nr; i += 2 )
138  tab.set(i) += c1 ;
139 }
140 
141 void coef_1d_cheb_odd (Array<double>& tab) {
142  assert (tab.get_ndim()==1) ;
143  int nr = tab.get_size(0) ;
144 
145  double* cf = new double[nr];
146  for (int i=0 ; i<nr ; i++)
147  cf[i] = tab(i) * sin(M_PI/2.*i/(nr-1)) ;
148 
149  double fmoins0 = - 0.5 * ( cf[0] - cf[nr-1] );
150 
151  auto & fftw_data = coef_1d_fftw(nr-1);
152  double* tab_auxi = fftw_data.buffer;
153  fftw_plan p = fftw_data.plan;
154 
155  for (int i=1 ; i<(nr-1)/2 ; i++) {
156  double fp = 0.5*(cf[i]+cf[nr-1-i]) ;
157  double fms = 0.5*(-cf[i]+cf[nr-1-i])*sin(M_PI*i/(nr-1)) ;
158  tab_auxi[i] = fp + fms ;
159  tab_auxi[nr-1-i] = fp-fms ;
160  }
161  tab_auxi[0] = 0.5*(cf[0] + cf[nr-1]) ;
162  tab_auxi[(nr-1)/2] = cf[(nr-1)/2] ;
163 
164  fftw_data.execute() ;
165 
166  // Coefficient pairs :
167  cf[0] = tab_auxi[0] / (nr-1) ;
168  for (int i=2; i<nr-1; i += 2)
169  cf[i] = 2*tab_auxi[(i/2)]/(nr-1) ;
170  cf[nr-1] = tab_auxi[(nr-1)/2] / (nr-1) ;
171 
172  // Coefficients impairs :
173  cf[1] = 0 ;
174  double som = 0;
175  for (int i = 3; i < nr; i += 2 ) {
176  cf[i] = cf[i-2] +4 * tab_auxi[nr-1-i/2]/(nr-1) ;
177  som += cf[i] ;
178  }
179 
180 // 2. Calcul de c_1 :
181  double c1 = ( fmoins0 - som ) / ((nr-1)/2) ;
182 
183 // 3. Coef. c_k avec k impair:
184  cf[1] = c1 ;
185  for (int i = 3; i < nr; i += 2 )
186  cf[i] += c1 ;
187 
188  cf[0] = 2* cf[0] ;
189  for (int i=1; i<nr-1; i++)
190  cf[i] = 2 * cf[i] - cf[i-1] ;
191  cf[nr-1] = 0 ;
192 
193  for (int i=0 ; i<nr ; i++)
194  tab.set(i) = cf[i] ;
195 
196  delete [] cf ;
197 }
198 
199 double coloc_leg(int,int) ;
200 double weight_leg(int,int) ;
201 double gamma_leg(int,int) ;
202 double leg(int,double) ;
203 double coloc_leg_parity(int,int) ;
204 double weight_leg_parity(int,int) ;
205 double gamma_leg_even(int,int) ;
206 double gamma_leg_odd(int,int) ;
207 
208 void coef_1d_leg (Array<double>& tab) {
209 
210  assert (tab.get_ndim()==1) ;
211  int nbr = tab.get_size(0) ;
212 
213  Array<double> res (nbr) ;
214  res = 0 ;
215  for (int i=0 ; i<nbr ; i++)
216  for (int j=0 ; j<nbr ; j++)
217  res.set(i) += tab(j)*leg(i, coloc_leg(j,nbr))*weight_leg(j,nbr)/gamma_leg(i,nbr) ;
218  tab = res ;
219 }
220 
221 void coef_1d_leg_even (Array<double>& tab) {
222 
223  assert (tab.get_ndim()==1) ;
224  int nbr = tab.get_size(0) ;
225 
226  Array<double> res (nbr) ;
227  res = 0 ;
228  for (int i=0 ; i<nbr ; i++)
229  for (int j=0 ; j<nbr ; j++)
230  res.set(i) += tab(j)*leg(i*2, coloc_leg_parity(j,nbr))
231  *weight_leg_parity(j,nbr)/gamma_leg_even(i,nbr) ;
232  tab = res ;
233 }
234 
235 void coef_1d_leg_odd (Array<double>& tab) {
236 
237  assert (tab.get_ndim()==1) ;
238  int nbr = tab.get_size(0) ;
239 
240  Array<double> res (nbr) ;
241  res = 0 ;
242  for (int i=0 ; i<nbr-1 ; i++)
243  for (int j=0 ; j<nbr ; j++)
244  res.set(i) += tab(j)*leg(i*2+1, coloc_leg_parity(j,nbr))
245  *weight_leg_parity(j,nbr)/gamma_leg_odd(i,nbr) ;
246  tab = res ;
247 }
248 
249 void coef_1d_cossin (Array<double>& tab) {
250 
251  assert (tab.get_ndim()==1) ;
252  int nbr = tab.get_size(0) ;
253  if (nbr>3) {
254 
255  auto & fftw_data = coef_1d_fftw(nbr-2);
256  double* cf = fftw_data.buffer;
257  fftw_plan p = fftw_data.plan;
258 
259  for (int i=0 ; i<nbr-2 ; i++)
260  cf[i] = tab(i) ;
261  fftw_data.execute() ;
262 
263  int index = 0 ;
264  double* pcos = cf ;
265  double* psin = cf+nbr-3 ;
266 
267  tab.set(index) = *pcos/double(nbr-2) ;
268  index++ ; pcos ++ ;
269  tab.set(index) = 0 ;
270  index ++ ;
271 
272  for (int i=1 ; i<nbr/2 - 1; i++) {
273  tab.set(index) = 2.*(*pcos)/double(nbr-2) ; // the cosines
274  index++ ; pcos ++ ;
275  tab.set(index) = -2.*(*psin)/double(nbr-2) ; // the sines
276  index ++ ; psin -- ;
277  }
278 
279  tab.set(index) = *pcos/double(nbr-2) ;
280  index ++ ;
281  tab.set(index) = 0 ;
282  }
283  else {
284  tab.set(1) = 0 ;
285  tab.set(2) = 0 ;
286  }
287 }
288 
289 void coef_1d_cos (Array<double>& tab) {
290 
291  assert (tab.get_ndim()==1) ;
292  int nbr = tab.get_size(0) ;
293  // Symetrie taken into account
294  double fmoins0 = 0.5 * ( tab(0) - tab(nbr-1) );
295 
296  auto & fftw_data = coef_1d_fftw(nbr-1);
297  double* tab_auxi = fftw_data.buffer;
298  fftw_plan p = fftw_data.plan;
299 
300  for (int i = 1; i < (nbr-1)/2 ; i++ ) {
301  double fp = 0.5 * ( tab(i) + tab(nbr-1-i) ) ;
302  double fms = 0.5 * ( tab(i) - tab(nbr-1-i) ) * sin(M_PI*i/(nbr-1)) ;
303  tab_auxi[i] = fp + fms ;
304  tab_auxi[nbr-1-i] = fp - fms ;
305  }
306 
307  tab_auxi[0] = 0.5 * ( tab(0) + tab(nbr-1) );
308  tab_auxi[(nbr-1)/2] = tab((nbr-1)/2);
309 
310  fftw_data.execute() ;
311 
312  tab.set(0) = tab_auxi[0] / (nbr-1) ;
313  for (int i=2; i<nbr-1; i += 2 )
314  tab.set(i) = 2*tab_auxi[i/2]/(nbr-1) ;
315  tab.set(nbr-1) = tab_auxi[(nbr-1)/2] / (nbr-1) ;
316 
317  tab.set(1) = 0 ;
318  double som = 0;
319  for (int i = 3; i < nbr; i += 2 ) {
320  tab.set(i) = tab(i-2) + 4 * tab_auxi[nbr-1 - i/2]/(nbr-1) ;
321  som += tab(i) ;
322  }
323 
324 // 2. Calcul de c_1 :
325  double c1 = ( fmoins0 - som ) / ((nbr-1)/2) ;
326 
327 // 3. Coef. c_k avec k impair:
328  tab.set(1) = c1 ;
329  for (int i = 3; i < nbr; i += 2 )
330  tab.set(i) += c1 ;
331 }
332 
333 void coef_1d_sin (Array<double>& tab) {
334 
335  assert (tab.get_ndim()==1) ;
336  int nbr = tab.get_size(0) ;
337  // Symetrie taken into account
338  auto & fftw_data = coef_1d_fftw(nbr-1);
339  double* tab_auxi = fftw_data.buffer;
340  fftw_plan p = fftw_data.plan;
341 
342  for (int i = 1; i < (nbr-1)/2 ; i++ ) {
343  double fp = 0.5 * ( tab(i) + tab(nbr-1-i) ) * sin(M_PI*i/(nbr-1)) ;
344  double fms = 0.5 * ( tab(i) - tab(nbr-1-i) ) ;
345  tab_auxi[i] = fp + fms ;
346  tab_auxi[nbr-1-i] = fp - fms ;
347  }
348 
349  tab_auxi[0] = 0.5 * ( tab(0) + tab(nbr-1) );
350  tab_auxi[(nbr-1)/2] = tab((nbr-1)/2);
351 
352  fftw_data.execute() ;
353 
354  tab.set(0) = 0 ;
355  for (int i=2; i<nbr-1; i += 2 )
356  tab.set(i) = -2*tab_auxi[nbr-1-i/2]/(nbr-1) ;
357  tab.set(nbr-1) = 0 ;
358 
359  tab.set(1) = 2*tab_auxi[0]/(nbr-1) ;
360  for (int i = 3; i < nbr; i += 2 )
361  tab.set(i) = tab(i-2) + 4 * tab_auxi[i/2]/(nbr-1) ;
362 }
363 
364 void coef_1d_cos_even (Array<double>& tab) {
365  assert (tab.get_ndim()==1) ;
366  int nbr = tab.get_size(0) ;
367  if (nbr>3) {
368  // Symetrie taken into account
369  double fmoins0 = 0.5 * ( tab(0) - tab(nbr-1) );
370 
371  auto & fftw_data = coef_1d_fftw(nbr-1);
372  double* tab_auxi = fftw_data.buffer;
373  fftw_plan p = fftw_data.plan;
374 
375  for (int i = 1; i < (nbr-1)/2 ; i++ ) {
376  double fp = 0.5 * ( tab(i) + tab(nbr-1-i) ) ;
377  double fms = 0.5 * ( tab(i) - tab(nbr-1-i) ) * sin(M_PI*i/(nbr-1)) ;
378  tab_auxi[i] = fp + fms ;
379  tab_auxi[nbr-1-i] = fp - fms ;
380  }
381 
382  tab_auxi[0] = 0.5 * ( tab(0) + tab(nbr-1) );
383  tab_auxi[(nbr-1)/2] = tab((nbr-1)/2);
384 
385  fftw_data.execute() ;
386 
387  tab.set(0) = tab_auxi[0] / (nbr-1) ;
388  for (int i=2; i<nbr-1; i += 2 )
389  tab.set(i) = 2*tab_auxi[i/2]/(nbr-1) ;
390  tab.set(nbr-1) = tab_auxi[(nbr-1)/2] / (nbr-1) ;
391 
392  tab.set(1) = 0 ;
393  double som = 0;
394  for (int i = 3; i < nbr; i += 2 ) {
395  tab.set(i) = tab(i-2) + 4 * tab_auxi[nbr-1 - i/2]/(nbr-1) ;
396  som += tab(i) ;
397  }
398 
399 // 2. Calcul de c_1 :
400  double c1 = ( fmoins0 - som ) / ((nbr-1)/2) ;
401 
402 // 3. Coef. c_k avec k impair:
403  tab.set(1) = c1 ;
404  for (int i = 3; i < nbr; i += 2 )
405  tab.set(i) += c1 ;
406  }
407 }
408 
409 void coef_1d_cos_odd (Array<double>& tab) {
410  assert (tab.get_ndim()==1) ;
411  int nbr = tab.get_size(0) ;
412  if (nbr>3) {
413  // Symetrie taken into account
414  double* cf = new double[nbr] ;
415  for (int i=0 ; i<nbr-1 ; i++)
416  cf[i] = tab(i)*sin(M_PI*(nbr-1-i)/2./(nbr-1)) ;
417  cf[nbr-1] = 0 ;
418  double fmoins0 = 0.5 * (cf[0] - cf[nbr-1] );
419 
420  auto & fftw_data = coef_1d_fftw(nbr-1);
421  double* tab_auxi = fftw_data.buffer;
422  fftw_plan p = fftw_data.plan;
423 
424  for (int i = 1; i < (nbr-1)/2 ; i++ ) {
425  double fp = 0.5 * ( cf[i] + cf[nbr-1-i] ) ;
426  double fms = 0.5 * ( cf[i] - cf[nbr-1-i] ) * sin(M_PI*i/(nbr-1)) ;
427  tab_auxi[i] = fp + fms ;
428  tab_auxi[nbr-1-i] = fp - fms ;
429  }
430 
431  tab_auxi[0] = 0.5 * ( cf[0] + cf[nbr-1] );
432  tab_auxi[(nbr-1)/2] = cf[(nbr-1)/2];
433 
434  fftw_data.execute() ;
435 
436  cf[0] = tab_auxi[0] / (nbr-1) ;
437  for (int i=2; i<nbr-1; i += 2 )
438  cf[i] = 2*tab_auxi[i/2]/(nbr-1) ;
439  cf[nbr-1] = tab_auxi[(nbr-1)/2] ;
440 
441  cf[1] = 0 ;
442  double som = 0;
443  for (int i = 3; i < nbr; i += 2 ) {
444  cf[i] = cf[i-2] + 4 * tab_auxi[nbr-1 - i/2]/(nbr-1) ;
445  som += cf[i] ;
446  }
447 
448 // 2. Calcul de c_1 :
449  double c1 = ( fmoins0 - som ) / ((nbr-1)/2) ;
450 
451 // 3. Coef. c_k avec k impair:
452  cf[1] = c1 ;
453  for (int i = 3; i < nbr; i += 2 )
454  cf[i] += c1 ;
455 
456  cf[0] = 2* cf[0] ;
457  for (int i=1; i<nbr-1; i++)
458  cf[i] = 2 * cf[i] - cf[i-1] ;
459  cf[nbr-1] = 0 ;
460  for (int i=0 ; i<nbr ; i++)
461  tab.set(i) = cf[i] ;
462 
463  delete [] cf ;
464  }
465 }
466 
467 void coef_1d_sin_even (Array<double>& tab) {
468  assert (tab.get_ndim()==1) ;
469  int nbr = tab.get_size(0) ;
470  if (nbr>3) {
471  // Symetrie taken into account
472  auto & fftw_data = coef_1d_fftw(nbr-1);
473  double* tab_auxi = fftw_data.buffer;
474  fftw_plan p = fftw_data.plan;
475 
476  for (int i = 1; i < (nbr-1)/2 ; i++ ) {
477  double fp = 0.5 * ( tab(i) + tab(nbr-1-i) ) * sin(M_PI*i/(nbr-1));
478  double fms = 0.5 * ( tab(i) - tab(nbr-1-i) ) ;
479  tab_auxi[i] = fp + fms ;
480  tab_auxi[nbr-1-i] = fp - fms ;
481  }
482 
483  tab_auxi[0] = 0.5 * ( tab(0) -tab(nbr-1) );
484  tab_auxi[(nbr-1)/2] = tab((nbr-1)/2);
485 
486  fftw_data.execute() ;
487 
488  tab.set(0) = 0. ;
489  for (int i=2; i<nbr-1; i += 2 )
490  tab.set(i) = -2*tab_auxi[nbr-1-i/2]/(nbr-1) ;
491  tab.set(nbr-1) = 0 ;
492 
493  tab.set(1) = 2*tab_auxi[0]/(nbr-1) ;
494  for (int i = 3; i < nbr; i += 2 )
495  tab.set(i) = tab(i-2) + 4 * tab_auxi[i/2]/(nbr-1) ;
496  }
497 }
498 
499 void coef_1d_sin_odd (Array<double>& tab) {
500  assert (tab.get_ndim()==1) ;
501  int nbr = tab.get_size(0) ;
502  if (nbr>3) {
503  // Symetrie taken into account
504  double* cf = new double[nbr] ;
505  cf[0] = 0 ;
506  for (int i=1 ; i<nbr ; i++)
507  cf[i] = tab(i)*sin(M_PI*i/2./(nbr-1)) ;
508  double fmoins0 = 0.5 * (cf[0] - cf[nbr-1] );
509 
510  auto & fftw_data = coef_1d_fftw(nbr-1);
511  double* tab_auxi = fftw_data.buffer;
512  fftw_plan p = fftw_data.plan;
513 
514  for (int i = 1; i < (nbr-1)/2 ; i++ ) {
515  double fp = 0.5 * ( cf[i] + cf[nbr-1-i] ) ;
516  double fms = 0.5 * ( cf[i] - cf[nbr-1-i] ) * sin(M_PI*i/(nbr-1)) ;
517  tab_auxi[i] = fp + fms ;
518  tab_auxi[nbr-1-i] = fp - fms ;
519  }
520 
521  tab_auxi[0] = 0.5 * ( cf[0] + cf[nbr-1] );
522  tab_auxi[(nbr-1)/2] = cf[(nbr-1)/2];
523 
524  fftw_data.execute() ;
525 
526  cf[0] = tab_auxi[0] / (nbr-1) ;
527  for (int i=2; i<nbr-1; i += 2 )
528  cf[i] = 2*tab_auxi[i/2]/(nbr-1) ;
529  cf[nbr-1] = tab_auxi[(nbr-1)/2] / (nbr-1) ;
530 
531  cf[1] = 0 ;
532  double som = 0;
533  for (int i = 3; i < nbr; i += 2 ) {
534  cf[i] = cf[i-2] + 4 * tab_auxi[nbr-1 - i/2]/(nbr-1) ;
535  som += cf[i] ;
536  }
537 
538 // 2. Calcul de c_1 :
539  double c1 = ( fmoins0 - som ) / ((nbr-1)/2) ;
540 
541 // 3. Coef. c_k avec k impair:
542  cf[1] = c1 ;
543  for (int i = 3; i < nbr; i += 2 )
544  cf[i] += c1 ;
545 
546  cf[0] = 2* cf[0] ;
547  for (int i=1; i<nbr-1; i++)
548  cf[i] = 2 * cf[i] + cf[i-1] ;
549  cf[nbr-1] = 0 ;
550  for (int i=0 ; i<nbr ; i++)
551  tab.set(i) = cf[i] ;
552 
553  delete [] cf ;
554  }
555 }
556 
557 
558 void coef_1d_cossin_even (Array<double>& tab) {
559 
560  // Copy values in double-size array :
561  assert (tab.get_ndim()==1) ;
562  int nbr = tab.get_size(0) ;
563  Array<double> tab2(2*nbr-2) ;
564  for (int i=0 ; i<nbr-2 ; i++) {
565  tab2.set(i) = tab(i) ;
566  tab2.set(i+nbr-2) = tab(i) ;
567  }
568 
569  coef_1d_cossin (tab2) ;
570 
571  int conte = 0 ;
572  for (int k=0 ; k<nbr-1 ; k+=2) {
573  tab.set(k) = tab2(conte) ;
574  tab.set(k+1) = tab2(conte+1) ;
575  conte += 4 ;
576  }
577 }
578 
579 
580 void coef_1d_cossin_odd (Array<double>& tab) {
581  // Copy values in double-size array :
582  assert (tab.get_ndim()==1) ;
583  int nbr = tab.get_size(0) ;
584  Array<double> tab2(2*nbr-2) ;
585  tab2 = 0 ;
586  for (int i=0 ; i<nbr-2 ; i++) {
587  tab2.set(i) = tab(i) ;
588  tab2.set(i+nbr-2) = -tab(i) ;
589  }
590  coef_1d_cossin (tab2) ;
591 
592  int conte = 2 ;
593  for (int k=0 ; k<nbr-3 ; k+=2) {
594  tab.set(k) = tab2(conte) ;
595  tab.set(k+1) = tab2(conte+1) ;
596  conte += 4 ;
597  }
598  tab.set(nbr-2) = 0. ;
599  tab.set(nbr-1) = 0. ;
600 }
601 
602 void coef_1d (int base, Array<double>& tab) {
603  static void (*coef_1d[NBR_MAX_BASE])(Array<double>&) ;
604  static bool premier_appel = true ;
605 
606  // Premier appel
607  if (premier_appel) {
608  premier_appel = false ;
609 
610  for (int i=0; i<NBR_MAX_BASE; i++)
611  coef_1d[i] = coef_1d_pasprevu ;
612 
613  coef_1d[CHEB] = coef_1d_cheb ;
614  coef_1d[CHEB_EVEN] = coef_1d_cheb_even ;
615  coef_1d[CHEB_ODD] = coef_1d_cheb_odd ;
616  coef_1d[COSSIN] = coef_1d_cossin ;
617  coef_1d[COS_EVEN] = coef_1d_cos_even ;
618  coef_1d[COS_ODD] = coef_1d_cos_odd ;
619  coef_1d[SIN_EVEN] = coef_1d_sin_even ;
620  coef_1d[SIN_ODD] = coef_1d_sin_odd ;
621  coef_1d[COS] = coef_1d_cos ;
622  coef_1d[SIN] = coef_1d_sin ;
623  coef_1d[LEG] = coef_1d_leg ;
624  coef_1d[LEG_EVEN] = coef_1d_leg_even ;
625  coef_1d[LEG_ODD] = coef_1d_leg_odd;
626  coef_1d[COSSIN_EVEN] = coef_1d_cossin_even ;
627  coef_1d[COSSIN_ODD] = coef_1d_cossin_odd ;
628  }
629 
630  coef_1d[base](tab) ;
631 }
632 }
633