KADATH
term_eq_bessel.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 "term_eq.hpp"
21 #include "scalar.hpp"
22 #include "tensor_impl.hpp"
23 #include <gsl/gsl_sf_bessel.h>
24 namespace Kadath {
25 Term_eq bessel_jl (const Term_eq& so, int ll) {
26 
27  int dom = so.get_dom() ;
28 
29  // Double case
30  if (so.type_data == TERM_D) {
31  double val = gsl_sf_bessel_jl (ll, *so.val_d) ;
32  bool doder = (so.der_d == 0x0) ? false : true ;
33  if (doder) {
34  double der = (ll==0) ? (cos (*so.val_d)/(*so.val_d) - sin(*so.val_d)/(*so.val_d)/(*so.val_d))*(*so.der_d) :
35  (gsl_sf_bessel_jl (ll-1,(*so.val_d)) - (ll+1)/(*so.val_d)*val)*(*so.der_d) ;
36  return Term_eq (dom, val, der) ;
37  }
38  else {
39  return Term_eq (dom, val) ;
40  }
41  }
42 
43  // Scalar case
44  if (so.type_data == TERM_T) {
45  if (so.val_t->get_valence()!=0) {
46  cerr << "Bessel function only defined for scalars so far..." << endl ;
47  abort() ;
48  }
49 
50  Val_domain val (so.val_t->get_space().get_domain(dom)) ;
51  val.allocate_conf() ;
52  Val_domain der (so.val_t->get_space().get_domain(dom)) ;
53  der.allocate_conf() ;
54 
55  bool doder = (so.der_t == 0x0) ? false : true ;
56  Index pos (val.get_domain()->get_nbr_points()) ;
57  do {
58  double vv = (*so.val_t)()(dom)(pos) ;
59  val.set(pos) = gsl_sf_bessel_jl (ll, vv) ;
60 
61  if (doder) {
62  double dd = (*so.der_t)()(dom)(pos) ;
63  if (ll==0)
64  der.set(pos) = (cos (vv)/vv - sin(vv)/vv/vv)*dd ;
65  else
66  der.set(pos) = (gsl_sf_bessel_jl (ll-1,vv) - (ll+1)/vv*val(pos))*dd ;
67  }
68  }
69  while (pos.inc()) ;
70 
71  val.std_base() ;
72  if (doder)
73  der.std_base() ;
74 
75  Scalar valres (so.val_t->get_space()) ;
76  valres.set_domain(dom) = val ;
77  if (doder) {
78  Scalar derres (so.val_t->get_space()) ;
79  derres.set_domain(dom) = der ;
80  return Term_eq (dom, valres, derres) ;
81  }
82  else {
83  return Term_eq (dom, valres) ;
84  }
85  }
86 
87  // Unknown case
88  cerr << "Unknown type of Term_eq in bessel_jl" << endl ;
89  abort() ;
90 }
91 
92 
93 Term_eq bessel_yl (const Term_eq& so, int ll) {
94 
95  int dom = so.get_dom() ;
96 
97  // Double case
98  if (so.type_data == TERM_D) {
99  double val = gsl_sf_bessel_yl (ll, *so.val_d) ;
100  bool doder = (so.der_d == 0x0) ? false : true ;
101  if (doder) {
102  double der = (ll==0) ? (sin (*so.val_d)/(*so.val_d) + cos(*so.val_d)/(*so.val_d)/(*so.val_d))*(*so.der_d) :
103  (gsl_sf_bessel_yl (ll-1,(*so.val_d)) - (ll+1)/(*so.val_d)*val)*(*so.der_d) ;
104  return Term_eq (dom, val, der) ;
105  }
106  else {
107  return Term_eq (dom, val) ;
108  }
109  }
110 
111  // Scalar case
112  if (so.type_data == TERM_T) {
113  if (so.val_t->get_valence()!=0) {
114  cerr << "Bessel function only defined for scalars so far..." << endl ;
115  abort() ;
116  }
117 
118  Val_domain val (so.val_t->get_space().get_domain(dom)) ;
119  val.allocate_conf() ;
120  Val_domain der (so.val_t->get_space().get_domain(dom)) ;
121  der.allocate_conf() ;
122 
123  bool doder = (so.der_t == 0x0) ? false : true ;
124  Index pos (val.get_domain()->get_nbr_points()) ;
125  do {
126  double vv = (*so.val_t)()(dom)(pos) ;
127  val.set(pos) = gsl_sf_bessel_yl (ll, vv) ;
128  if (doder) {
129  double dd = (*so.der_t)()(dom)(pos) ;
130  if (ll==0)
131  der.set(pos) = (sin (vv)/vv + cos(vv)/vv/vv)*dd ;
132  else
133  der.set(pos) = (gsl_sf_bessel_yl (ll-1,vv) - (ll+1)/vv*val(pos))*dd ;
134  }
135  }
136  while (pos.inc()) ;
137 
138  val.std_base() ;
139  if (doder)
140  der.std_base() ;
141 
142  Scalar valres (so.val_t->get_space()) ;
143  valres.set_domain(dom) = val ;
144  if (doder) {
145  Scalar derres (so.val_t->get_space()) ;
146  derres.set_domain(dom) = der ;
147  return Term_eq (dom, valres, derres) ;
148  }
149  else {
150  return Term_eq (dom, valres) ;
151  }
152  }
153 
154  // Unknown case
155  cerr << "Unknown type of Term_eq in bessel_yl" << endl ;
156  abort() ;
157 }
158 
159 Term_eq bessel_djl (const Term_eq& so, int ll) {
160 
161  int dom = so.get_dom() ;
162 
163  // Double case
164  if (so.type_data == TERM_D) {
165  double val ;
166 
167  double jl = gsl_sf_bessel_jl (ll,(*so.val_d)) ;
168  double jlm1 = (ll==0) ? 0 : gsl_sf_bessel_jl (ll-1 ,(*so.val_d)) ;
169  double jlm2 = (ll<2) ? 0 : gsl_sf_bessel_jl (ll-2 ,(*so.val_d)) ;
170 
171  if (ll==0) {
172  val = cos(*so.val_d)/(*so.val_d) - sin(*so.val_d)/(*so.val_d)/(*so.val_d) ;
173  }
174  else {
175  val = jlm1 - (ll+1)/(*so.val_d)*jl ;
176  }
177  bool doder = (so.der_d == 0x0) ? false : true ;
178  if (doder) {
179  double der ;
180 
181  if (ll==0) {
182  der = (-sin(*so.val_d)/(*so.val_d) -2*cos(*so.val_d)/(*so.val_d)/(*so.val_d) + 2*sin(*so.val_d)/(*so.val_d)/(*so.val_d)/(*so.val_d))*(*so.der_d) ;
183  }
184  else {
185  if (ll==1) {
186  der = (cos(*so.val_d)/(*so.val_d) -3*sin(*so.val_d)/(*so.val_d)/(*so.val_d) -6*cos(*so.val_d)/(*so.val_d)/(*so.val_d)/(*so.val_d) + 6*sin(*so.val_d) /(*so.val_d)/(*so.val_d)/(*so.val_d)/(*so.val_d))*(*so.der_d) ;
187  }
188  else {
189  der = (jlm2 - double(2*ll+1)*jlm1/(*so.val_d) + double((ll+1)*(ll+2))*jl/(*so.val_d)/(*so.val_d))*(*so.der_d) ;
190  }
191  }
192  return Term_eq (dom, val, der) ;
193  }
194  else {
195  return Term_eq (dom, val) ;
196  }
197  }
198 
199  // Scalar case
200  if (so.type_data == TERM_T) {
201  if (so.val_t->get_valence()!=0) {
202  cerr << "Bessel function only defined for scalars so far..." << endl ;
203  abort() ;
204  }
205 
206  Val_domain val (so.val_t->get_space().get_domain(dom)) ;
207  val.allocate_conf() ;
208  Val_domain der (so.val_t->get_space().get_domain(dom)) ;
209  der.allocate_conf() ;
210 
211  bool doder = (so.der_t == 0x0) ? false : true ;
212  Index pos (val.get_domain()->get_nbr_points()) ;
213  do {
214  double vv = (*so.val_t)()(dom)(pos) ;
215 
216  double jl = gsl_sf_bessel_jl (ll, vv) ;
217  double jlm1 = (ll==0) ? 0 : gsl_sf_bessel_jl (ll-1 , vv) ;
218  double jlm2 = (ll<2) ? 0 : gsl_sf_bessel_jl (ll-2 , vv) ;
219 
220  if (ll==0) {
221  val.set(pos) = cos(vv)/vv - sin(vv)/vv/vv ;
222  }
223  else {
224  val.set(pos) = jlm1 - (ll+1)/vv*jl ;
225  }
226 
227  if (doder) {
228  double dd = (*so.der_t)()(dom)(pos) ;
229  if (ll==0) {
230  der.set(pos) = (-sin(vv)/vv -2*cos(vv)/vv/vv + 2*sin(vv)/vv/vv/vv)*dd ;
231  }
232  else {
233  if (ll==1) {
234  der.set(pos) = (cos(vv)/vv -3*sin(vv)/vv/vv -6*cos(vv)/vv/vv/vv + 6*sin(vv) /vv/vv/vv/vv)*dd ;
235  }
236  else {
237  der = (jlm2 - double(2*ll+1)*jlm1/vv + double((ll+1)*(ll+2))*jl/vv/vv)*dd ;
238  }
239 
240  }
241  }
242  }
243  while (pos.inc()) ;
244 
245  val.std_base() ;
246  if (doder)
247  der.std_base() ;
248 
249  Scalar valres (so.val_t->get_space()) ;
250  valres.set_domain(dom) = val ;
251  if (doder) {
252  Scalar derres (so.val_t->get_space()) ;
253  derres.set_domain(dom) = der ;
254  return Term_eq (dom, valres, derres) ;
255  }
256  else {
257  return Term_eq (dom, valres) ;
258  }
259  }
260 
261  // Unknown case
262  cerr << "Unknown type of Term_eq in bessel_djl" << endl ;
263  abort() ;
264 }
265 
266 
267 Term_eq bessel_dyl (const Term_eq& so, int ll) {
268 
269  int dom = so.get_dom() ;
270 
271  // Double case
272  if (so.type_data == TERM_D) {
273  double val ;
274 
275  double yl = gsl_sf_bessel_yl (ll,(*so.val_d)) ;
276  double ylm1 = (ll==0) ? 0 : gsl_sf_bessel_yl (ll-1 ,(*so.val_d)) ;
277  double ylm2 = (ll<2) ? 0 : gsl_sf_bessel_yl (ll-2 ,(*so.val_d)) ;
278 
279  if (ll==0) {
280  val = sin(*so.val_d)/(*so.val_d) + cos(*so.val_d)/(*so.val_d)/(*so.val_d) ;
281  }
282  else {
283  val = ylm1 - (ll+1)/(*so.val_d)*yl ;
284  }
285  bool doder = (so.der_d == 0x0) ? false : true ;
286  if (doder) {
287  double der ;
288 
289  if (ll==0) {
290  der = (cos(*so.val_d)/(*so.val_d) -2*sin(*so.val_d)/(*so.val_d)/(*so.val_d) - 2*cos(*so.val_d)/(*so.val_d)/(*so.val_d)/(*so.val_d))*(*so.der_d) ;
291  }
292  else {
293  if (ll==1) {
294  der = (sin(*so.val_d)/(*so.val_d) +3*cos(*so.val_d)/(*so.val_d)/(*so.val_d) -6*sin(*so.val_d)/(*so.val_d)/(*so.val_d)/(*so.val_d) -6*cos(*so.val_d) /(*so.val_d)/(*so.val_d)/(*so.val_d)/(*so.val_d))*(*so.der_d) ;
295  }
296  else {
297  der = (ylm2 - double(2*ll+1)*ylm1/(*so.val_d) + double((ll+1)*(ll+2))*yl/(*so.val_d)/(*so.val_d))*(*so.der_d) ;
298  }
299  }
300  return Term_eq (dom, val, der) ;
301  }
302  else {
303  return Term_eq (dom, val) ;
304  }
305  }
306 
307  // Scalar case
308  if (so.type_data == TERM_T) {
309  if (so.val_t->get_valence()!=0) {
310  cerr << "Bessel function only defined for scalars so far..." << endl ;
311  abort() ;
312  }
313 
314  Val_domain val (so.val_t->get_space().get_domain(dom)) ;
315  val.allocate_conf() ;
316  Val_domain der (so.val_t->get_space().get_domain(dom)) ;
317  der.allocate_conf() ;
318 
319  bool doder = (so.der_t == 0x0) ? false : true ;
320  Index pos (val.get_domain()->get_nbr_points()) ;
321  do {
322  double vv = (*so.val_t)()(dom)(pos) ;
323 
324  double yl = gsl_sf_bessel_yl (ll, vv) ;
325  double ylm1 = (ll==0) ? 0 : gsl_sf_bessel_yl (ll-1 , vv) ;
326  double ylm2 = (ll<2) ? 0 : gsl_sf_bessel_yl (ll-2 , vv) ;
327 
328  if (ll==0) {
329  val.set(pos) = sin(vv)/vv + cos(vv)/vv/vv ;
330  }
331  else {
332  val.set(pos) = ylm1 - (ll+1)/vv*yl ;
333  }
334 
335  if (doder) {
336  double dd = (*so.der_t)()(dom)(pos) ;
337  if (ll==0) {
338  der.set(pos) = (cos(vv)/vv -2*sin(vv)/vv/vv - 2*cos(vv)/vv/vv/vv)*dd ;
339  }
340  else {
341  if (ll==1) {
342  der.set(pos) = (sin(vv)/vv +3*cos(vv)/vv/vv -6*sin(vv)/vv/vv/vv - 6*cos(vv) /vv/vv/vv/vv)*dd ;
343  }
344  else {
345  der = (ylm2 - double(2*ll+1)*ylm1/vv + double((ll+1)*(ll+2))*yl/vv/vv)*dd ;
346  }
347 
348  }
349  }
350  }
351  while (pos.inc()) ;
352 
353  val.std_base() ;
354  if (doder)
355  der.std_base() ;
356 
357  Scalar valres (so.val_t->get_space()) ;
358  valres.set_domain(dom) = val ;
359  if (doder) {
360  Scalar derres (so.val_t->get_space()) ;
361  derres.set_domain(dom) = der ;
362  return Term_eq (dom, valres, derres) ;
363  }
364  else {
365  return Term_eq (dom, valres) ;
366  }
367  }
368 
369  // Unknown case
370  cerr << "Unknown type of Term_eq in bessel_dyl" << endl ;
371  abort() ;
372 }}
Kadath::Domain::get_nbr_points
Dim_array const & get_nbr_points() const
Returns the number of points.
Definition: space.hpp:92
Kadath::Index
Class that gives the position inside a multi-dimensional Array.
Definition: index.hpp:38
Kadath::Index::inc
bool inc(int increm, int var=0)
Increments the position of the Index.
Definition: index.hpp:99
Kadath::Scalar
The class Scalar does not really implements scalars in the mathematical sense but rather tensorial co...
Definition: scalar.hpp:67
Kadath::Scalar::set_domain
Val_domain & set_domain(int)
Read/write of a particular Val_domain.
Definition: scalar.hpp:555
Kadath::Space::get_domain
const Domain * get_domain(int i) const
returns a pointer on the domain.
Definition: space.hpp:1385
Kadath::Tensor::get_valence
int get_valence() const
Returns the valence.
Definition: tensor.hpp:509
Kadath::Tensor::get_space
const Space & get_space() const
Returns the Space.
Definition: tensor.hpp:499
Kadath::Term_eq
This class is intended to describe the manage objects appearing in the equations.
Definition: term_eq.hpp:62
Kadath::Term_eq::der_t
Tensor * der_t
Pointer on the variation, if the Term_eq is a Tensor.
Definition: term_eq.hpp:69
Kadath::Term_eq::type_data
const int type_data
Flag describing the type of data :
Definition: term_eq.hpp:75
Kadath::Term_eq::val_d
double * val_d
Pointer on the value, if the Term_eq is a double.
Definition: term_eq.hpp:66
Kadath::Term_eq::get_dom
int get_dom() const
Definition: term_eq.hpp:135
Kadath::Term_eq::der_d
double * der_d
Pointer on the variation if the Term_eq is a double.
Definition: term_eq.hpp:67
Kadath::Term_eq::val_t
Tensor * val_t
Pointer on the value, if the Term_eq is a Tensor.
Definition: term_eq.hpp:68
Kadath::Val_domain
Class for storing the basis of decompositions of a field and its values on both the configuration and...
Definition: val_domain.hpp:69
Kadath::Val_domain::set
double & set(const Index &pos)
Read/write the value of the field in the configuration space.
Definition: val_domain.cpp:171
Kadath::Val_domain::std_base
void std_base()
Sets the standard basis of decomposition.
Definition: val_domain.cpp:246
Kadath::Val_domain::allocate_conf
void allocate_conf()
Allocates the values in the configuration space and destroys the values in the coefficients space.
Definition: val_domain.cpp:209
Kadath::Val_domain::get_domain
const Domain * get_domain() const
Definition: val_domain.hpp:111