KADATH
integral_1d.cpp
1 /*
2  Copyright 2017 Philippe Grandclement
3 
4  This file is part of Kadath.
5 
6  Kadath 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 3 of the License, or
9  (at your option) any later version.
10 
11  Kadath 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 Kadath. If not, see <http://www.gnu.org/licenses/>.
18 */
19 
20 #include "point.hpp"
21 #include "base_spectral.hpp"
22 #include "array.hpp"
23 namespace Kadath {
24 double integral_1d_pasprevu ( const Array<double>&) {
25  cout << "Integral_1d not implemented." << endl ;
26  abort() ;
27 }
28 
29 double integral_1d_cos_even (const Array<double>& so) {
30  return M_PI*so(0) ;
31 }
32 
33 double integral_1d_sin_odd (const Array<double>& so) {
34  int nn = so.get_size(0) ;
35  double res = 0. ;
36  for (int i=0 ; i<nn ; i++)
37  res += so(i)*2./(2.*i+1) ;
38  return res ;
39 }
40 
41 double integral_1d_sin_even (const Array<double>&) {
42  return 0. ;
43 }
44 
45 double integral_1d_cos_odd (const Array<double>&) {
46  return 0. ;
47 }
48 
49 double integral_1d_cheb_even (const Array<double>& so) {
50  int nn = so.get_size(0) ;
51  double res = 0. ;
52  for (int i=0 ; i<nn ; i++)
53  res += -so(i)/(4.*i*i-1) ;
54  return res ;
55 }
56 double integral_1d_cheb (const Array<double>& so) {
57  int nn = so.get_size(0) ;
58  double res = 0. ;
59  for (int i=0 ; i<nn ; i+=2)
60  res += -2*so(i)/(i*i-1.) ;
61  return res ;
62 }
63 
64 double integral_1d (int base, const Array<double>& tab) {
65 
66  static double (*integral_1d[NBR_MAX_BASE])(const Array<double>&) ;
67  static bool premier_appel = true ;
68  // Premier appel
69  if (premier_appel) {
70  premier_appel = false ;
71 
72  for (int i=0; i<NBR_MAX_BASE; i++)
73  integral_1d[i] = integral_1d_pasprevu ;
74 
75  integral_1d[COS_EVEN] = integral_1d_cos_even ;
76  integral_1d[COS_ODD] = integral_1d_cos_odd ;
77  integral_1d[SIN_ODD] = integral_1d_sin_odd ;
78  integral_1d[SIN_EVEN] = integral_1d_sin_even ;
79  integral_1d[CHEB_EVEN] = integral_1d_cheb_even ;
80  integral_1d[CHEB] = integral_1d_cheb ;
81  }
82 
83  return integral_1d[base](tab) ;
84 }}