KADATH
der_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 "base_spectral.hpp"
21 #include "headcpp.hpp"
22 #include "array.hpp"
23 #include "array_math.hpp"
24 namespace Kadath {
25 int der_1d_pasprevu (Array<double>&) {
26  cout << "Der_1d not implemented." << endl ;
27  abort() ;
28 }
29 
30  // CHEB :
31 int der_1d_cheb (Array<double>& tab) {
32 
33  assert (tab.get_ndim()==1) ;
34  int nr = tab.get_size(0) ;
35 
36  Array<double> res (nr) ;
37  res = 0 ;
38  for (int i=0 ;i<nr ; i++)
39  for (int p=i+1 ; p<nr ; p+=2)
40  res.set(i) += p*tab(p) ;
41  res *= 2 ;
42  res.set(0) /= 2 ;
43  tab = res ;
44 
45 
46  // Output base :
47  return CHEB ;
48 }
49  // CHEB EVEN:
50 int der_1d_cheb_even (Array<double>& tab) {
51 
52  assert (tab.get_ndim()==1) ;
53  int nr = tab.get_size(0) ;
54  Array<double> res (nr) ;
55  res = 0 ;
56  for (int i=0 ;i<nr ; i++)
57  for (int p=i+1 ; p<nr ; p++)
58  res.set(i) += 2*p*tab(p) ;
59  res.set(nr-1) = 0 ;
60  res *= 2 ;
61  tab = res ;
62 
63  // Output base :
64  return CHEB_ODD ;
65 }
66 
67  // CHEB ODD:
68 int der_1d_cheb_odd (Array<double>& tab) {
69 
70  assert (tab.get_ndim()==1) ;
71  int nr = tab.get_size(0) ;
72  Array<double> res (nr) ;
73  res = 0 ;
74  for (int i=0 ;i<nr ; i++)
75  for (int p=i ; p<nr ; p++)
76  res.set(i) += (2*p+1)*tab(p) ;
77  res *= 2 ;
78  res.set(0) /= 2 ;
79  tab = res ;
80 
81  // Output base :
82  return CHEB_EVEN ;
83 }
84 
85  // LEG :
86 int der_1d_leg (Array<double>& tab) {
87 
88  assert (tab.get_ndim()==1) ;
89  int nr = tab.get_size(0) ;
90  Array<double> res (nr) ;
91  res = 0 ;
92  for (int i=0 ;i<nr ; i++)
93  for (int p=i+1 ; p<nr ; p+=2)
94  res.set(i) += (2*i+1)*tab(p) ;
95  tab = res ;
96 
97  // Output base :
98  return LEG ;
99 }
100  // LEG EVEN:
101 int der_1d_leg_even (Array<double>& tab) {
102 
103  assert (tab.get_ndim()==1) ;
104  int nr = tab.get_size(0) ;
105  Array<double> res (nr) ;
106  res = 0 ;
107  for (int i=0 ;i<nr ; i++)
108  for (int p=i+1 ; p<nr ; p++)
109  res.set(i) += (4*i+3)*tab(p) ;
110  res.set(nr-1) = 0 ;
111  tab = res ;
112 
113  // Output base :
114  return LEG_ODD ;
115 }
116 
117  // LEG ODD:
118 int der_1d_leg_odd (Array<double>& tab) {
119 
120  assert (tab.get_ndim()==1) ;
121  int nr = tab.get_size(0) ;
122  Array<double> res (nr) ;
123  res = 0 ;
124  for (int i=0 ;i<nr ; i++)
125  for (int p=i ; p<nr ; p++)
126  res.set(i) += (4*i+1)*tab(p) ;
127  tab = res ;
128 
129  // Output base :
130  return LEG_EVEN ;
131 }
132 
133  // COSSIN:
134 int der_1d_cossin (Array<double>& tab) {
135 
136  assert (tab.get_ndim()==1) ;
137  int nr = tab.get_size(0) ;
138  Array<double> res (nr) ;
139  // cosines
140  for (int i=0 ;i<nr-1 ; i+=2)
141  res.set(i+1) = - (i/2)*tab(i) ;
142  // Sines
143  for (int i=1 ;i<nr ; i+=2)
144  res.set(i-1) = ((i-1)/2)*tab(i) ;
145  res.set(nr-1) = 0 ;
146  tab = res ;
147 
148  // Output base :
149  return COSSIN ;
150 }
151 
152  // COSSIN EVEN :
153 int der_1d_cossin_even (Array<double>& tab) {
154 
155  assert (tab.get_ndim()==1) ;
156  int nr = tab.get_size(0) ;
157  Array<double> res (nr) ;
158  // cosines
159  for (int i=0 ;i<nr-1 ; i+=2)
160  res.set(i+1) = - i*tab(i) ;
161  // Sines
162  for (int i=1 ;i<nr ; i+=2)
163  res.set(i-1) = (i-1)*tab(i) ;
164  tab = res ;
165 
166  // Output base :
167  return COSSIN_EVEN ;
168 }
169 
170  // COSSIN ODD :
171 int der_1d_cossin_odd (Array<double>& tab) {
172 
173  assert (tab.get_ndim()==1) ;
174  int nr = tab.get_size(0) ;
175  Array<double> res (nr) ;
176  // cosines
177  for (int i=0 ;i<nr-1 ; i+=2)
178  res.set(i+1) = - (i+1)*tab(i) ;
179  // Sines
180  for (int i=1 ;i<nr ; i+=2)
181  res.set(i-1) = i*tab(i) ;
182  tab = res ;
183 
184  // Output base :
185  return COSSIN_ODD ;
186 }
187 
188  // COS:
189 int der_1d_cos (Array<double>& tab) {
190 
191  assert (tab.get_ndim()==1) ;
192  int nr = tab.get_size(0) ;
193  Array<double> res (nr) ;
194 
195  for (int i=0 ;i<nr ; i++)
196  res.set(i) = - i*tab(i) ;
197  tab = res ;
198 
199  // Output base :
200  return SIN ;
201 }
202 
203 
204  // SIN:
205 int der_1d_sin(Array<double>& tab) {
206 
207  assert (tab.get_ndim()==1) ;
208  int nr = tab.get_size(0) ;
209  Array<double> res (nr) ;
210 
211  for (int i=0 ;i<nr ; i++)
212  res.set(i) = i*tab(i) ;
213  tab = res ;
214 
215  // Output base :
216  return COS ;
217 }
218 
219 
220 
221  // COS EVEN:
222 int der_1d_cos_even (Array<double>& tab) {
223 
224  assert (tab.get_ndim()==1) ;
225  int nr = tab.get_size(0) ;
226  Array<double> res (nr) ;
227 
228  for (int i=0 ;i<nr ; i++)
229  res.set(i) = -2*i*tab(i) ;
230  res.set(nr-1) = 0 ;
231  tab = res ;
232 
233  // Output base :
234  return SIN_EVEN ;
235 }
236 
237 
238  // COS ODD:
239 int der_1d_cos_odd (Array<double>& tab) {
240 
241  assert (tab.get_ndim()==1) ;
242  int nr = tab.get_size(0) ;
243  Array<double> res (nr) ;
244 
245  for (int i=0 ;i<nr ; i++)
246  res.set(i) = -(2*i+1)*tab(i) ;
247  tab = res ;
248 
249  // Output base :
250  return SIN_ODD ;
251 }
252 
253 
254 
255 
256  // SIN EVEN:
257 int der_1d_sin_even (Array<double>& tab) {
258 
259  assert (tab.get_ndim()==1) ;
260  int nr = tab.get_size(0) ;
261  Array<double> res (nr) ;
262 
263  for (int i=0 ;i<nr ; i++)
264  res.set(i) = 2*i*tab(i) ;
265  tab = res ;
266 
267  // Output base :
268  return COS_EVEN ;
269 }
270 
271 
272  // SIN ODD:
273 int der_1d_sin_odd (Array<double>& tab) {
274 
275  assert (tab.get_ndim()==1) ;
276  int nr = tab.get_size(0) ;
277  Array<double> res (nr) ;
278 
279  for (int i=0 ;i<nr ; i++)
280  res.set(i) = (2*i+1)*tab(i) ;
281  res.set(nr-1) = 0 ;
282  tab = res ;
283 
284  // Output base :
285  return COS_ODD ;
286 }
287 
288 
289 
290 int der_1d (int base, Array<double>& tab) {
291  static int (*der_1d[NBR_MAX_BASE])(Array<double>&) ;
292  static bool premier_appel = true ;
293 
294  // Premier appel
295  if (premier_appel) {
296  premier_appel = false ;
297 
298  for (int i=0; i<NBR_MAX_BASE; i++)
299  der_1d[i] = der_1d_pasprevu ;
300 
301  der_1d[CHEB] = der_1d_cheb ;
302  der_1d[CHEB_EVEN] = der_1d_cheb_even ;
303  der_1d[CHEB_ODD] = der_1d_cheb_odd ;
304  der_1d[COSSIN] = der_1d_cossin ;
305  der_1d[COS_EVEN] = der_1d_cos_even ;
306  der_1d[COS_ODD] = der_1d_cos_odd ;
307  der_1d[SIN_EVEN] = der_1d_sin_even ;
308  der_1d[SIN_ODD] = der_1d_sin_odd ;
309  der_1d[COS] = der_1d_cos ;
310  der_1d[SIN] = der_1d_sin ;
311  der_1d[LEG] = der_1d_leg ;
312  der_1d[LEG_EVEN] = der_1d_leg_even ;
313  der_1d[LEG_ODD] = der_1d_leg_odd;
314  der_1d[COSSIN_EVEN] = der_1d_cossin_even ;
315  der_1d[COSSIN_ODD] = der_1d_cossin_odd ;
316  }
317 
318  return der_1d[base](tab) ;
319 }}