1 // GetDP - Copyright (C) 1997-2021 P. Dular and C. Geuzaine, University of Liege
2 //
3 // See the LICENSE.txt file for license information. Please report all
4 // issues on https://gitlab.onelab.info/getdp/getdp/issues.
5
6 // g++ -std=c++11 on mingw does not define bessel functions
7 #if defined (WIN32) && !defined(__CYGWIN__)
8 #undef __STRICT_ANSI__
9 #endif
10
11 #include <math.h>
12 #include <complex>
13 #include "ProData.h"
14 #include "F.h"
15 #include "MallocUtils.h"
16 #include "Message.h"
17 #include "Bessel.h"
18
19 extern struct CurrentData Current ;
20
21 /* ------------------------------------------------------------------------ */
22 /* C math functions (scalar, 1 argument, imaginary part set to zero) */
23 /* ------------------------------------------------------------------------ */
24
25 #define scalar_real_1_arg(func, string) \
26 int k; \
27 \
28 if(A->Type != SCALAR) \
29 Message::Error("Non scalar argument for function '" string "'"); \
30 if(Current.NbrHar > 1 && A->Val[MAX_DIM] != 0.) \
31 Message::Error("Function '" string "' only accepts real arguments");\
32 V->Val[0] = func(A->Val[0]) ; \
33 if (Current.NbrHar > 1){ \
34 V->Val[MAX_DIM] = 0. ; \
35 for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) \
36 V->Val[MAX_DIM*k] = V->Val[MAX_DIM*(k+1)] = 0. ; \
37 } \
38 V->Type = SCALAR;
39
F_Floor(F_ARG)40 void F_Floor (F_ARG) { scalar_real_1_arg (floor,"Floor") }
F_Ceil(F_ARG)41 void F_Ceil (F_ARG) { scalar_real_1_arg (ceil, "Ceil") }
F_Fabs(F_ARG)42 void F_Fabs (F_ARG) { scalar_real_1_arg (fabs, "Fabs") }
F_Asin(F_ARG)43 void F_Asin (F_ARG) { scalar_real_1_arg (asin, "Asin") }
F_Acos(F_ARG)44 void F_Acos (F_ARG) { scalar_real_1_arg (acos, "Acos") }
F_Atan(F_ARG)45 void F_Atan (F_ARG) { scalar_real_1_arg (atan, "Atan") }
F_Atanh(F_ARG)46 void F_Atanh (F_ARG) { scalar_real_1_arg (atanh, "Atanh")}
47
48 #undef scalar_real_1_arg
49
50 /* ------------------------------------------------------------------------ */
51 /* C math functions (complex scalar, 1 argument) */
52 /* ------------------------------------------------------------------------ */
53
54 #define scalar_cmplx_1_arg(func, string) \
55 int k; \
56 std::complex<double> tmp; \
57 \
58 if(A->Type != SCALAR) \
59 Message::Error("Non scalar argument for function '" string "'"); \
60 \
61 switch(Current.NbrHar){ \
62 case 1: \
63 V->Val[0] = func(A->Val[0]) ; \
64 V->Val[MAX_DIM] = 0. ; \
65 break ; \
66 case 2: \
67 tmp = std::complex<double>(A->Val[0], A->Val[MAX_DIM]); \
68 tmp = func(tmp); \
69 V->Val[0] = std::real(tmp); \
70 V->Val[MAX_DIM] = std::imag(tmp); \
71 break; \
72 default: \
73 tmp = std::complex<double>(A->Val[0], A->Val[MAX_DIM]); \
74 tmp = func(tmp); \
75 V->Val[0] = std::real(tmp); \
76 V->Val[MAX_DIM] = std::imag(tmp); \
77 for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) { \
78 V->Val[MAX_DIM*k] = std::real(tmp); \
79 V->Val[MAX_DIM*(k+1)] = std::imag(tmp); \
80 } \
81 } \
82 V->Type = SCALAR; \
83
F_Exp(F_ARG)84 void F_Exp (F_ARG) { scalar_cmplx_1_arg (exp, "Exp") }
F_Log(F_ARG)85 void F_Log (F_ARG) { scalar_cmplx_1_arg (log, "Log") }
F_Log10(F_ARG)86 void F_Log10 (F_ARG) { scalar_cmplx_1_arg (log10,"Log10") }
F_Sqrt(F_ARG)87 void F_Sqrt (F_ARG) { scalar_cmplx_1_arg (sqrt, "Sqrt") }
F_Sin(F_ARG)88 void F_Sin (F_ARG) { scalar_cmplx_1_arg (sin, "Sin") }
F_Cos(F_ARG)89 void F_Cos (F_ARG) { scalar_cmplx_1_arg (cos, "Cos" ) }
F_Tan(F_ARG)90 void F_Tan (F_ARG) { scalar_cmplx_1_arg (tan, "Tan") }
F_Sinh(F_ARG)91 void F_Sinh (F_ARG) { scalar_cmplx_1_arg (sinh, "Sinh") }
F_Cosh(F_ARG)92 void F_Cosh (F_ARG) { scalar_cmplx_1_arg (cosh, "Cosh") }
F_Tanh(F_ARG)93 void F_Tanh (F_ARG) { scalar_cmplx_1_arg (tanh, "Tanh") }
F_Abs(F_ARG)94 void F_Abs (F_ARG) { scalar_cmplx_1_arg (std::abs, "Abs") }
95
96 #undef scalar_cmplx_1_arg
97
98 /* ------------------------------------------------------------------------ */
99 /* C math functions (scalar, 2 arguments, imaginary part set to zero) */
100 /* ------------------------------------------------------------------------ */
101
102 #define scalar_real_2_arg(func, string) \
103 int k; \
104 \
105 if(A->Type != SCALAR || (A+1)->Type != SCALAR) \
106 Message::Error("Non scalar argument(s) for function '" string "'"); \
107 if(Current.NbrHar > 1 && (A->Val[MAX_DIM] != 0. || \
108 (A+1)->Val[MAX_DIM] != 0.)) \
109 Message::Error("Function '" string "' only accepts real arguments");\
110 \
111 V->Val[0] = func(A->Val[0], (A+1)->Val[0]) ; \
112 if (Current.NbrHar > 1){ \
113 V->Val[MAX_DIM] = 0. ; \
114 for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) { \
115 V->Val[MAX_DIM*k] = func(A->Val[0], (A+1)->Val[0]) ; \
116 V->Val[MAX_DIM*(k+1)] = 0. ; \
117 } \
118 } \
119 V->Type = SCALAR;
120
F_Atan2(F_ARG)121 void F_Atan2 (F_ARG) { scalar_real_2_arg (atan2, "Atan2") }
F_Fmod(F_ARG)122 void F_Fmod (F_ARG) { scalar_real_2_arg (fmod, "Fmod") }
123
124 #undef scalar_real_2_arg
125
126 /* ------------------------------------------------------------------------ */
127 /* Sign function */
128 /* ------------------------------------------------------------------------ */
129
F_Sign(F_ARG)130 void F_Sign(F_ARG)
131 {
132 int k;
133 double x;
134
135 if(A->Type != SCALAR)
136 Message::Error("Non scalar argument for function 'Sign'");
137 x = A->Val[0];
138
139 if(x >= 0.)
140 V->Val[0] = 1.;
141 else if(x < 0.)
142 V->Val[0] = -1.;
143 else
144 V->Val[0] = 0.;
145
146 if (Current.NbrHar > 1){
147 V->Val[MAX_DIM] = 0. ;
148 for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2)
149 V->Val[MAX_DIM*k] = V->Val[MAX_DIM*(k+1)] = 0. ;
150 }
151 V->Type = SCALAR;
152 }
153
154 /* ------------------------------------------------------------------------ */
155 /* Min, Max */
156 /* ------------------------------------------------------------------------ */
157
F_Min(F_ARG)158 void F_Min(F_ARG)
159 {
160 int k;
161
162 if(A->Type != SCALAR || (A+1)->Type != SCALAR)
163 Message::Error("Non scalar argument(s) for function Min");
164
165 V->Val[0] = (A->Val[0] < (A+1)->Val[0]) ? A->Val[0] : (A+1)->Val[0];
166 if (Current.NbrHar > 1){
167 V->Val[MAX_DIM] = 0. ;
168 for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2)
169 V->Val[MAX_DIM*k] = V->Val[MAX_DIM*(k+1)] = 0. ;
170 }
171 V->Type = SCALAR;
172 }
173
F_Max(F_ARG)174 void F_Max(F_ARG)
175 {
176 int k;
177
178 if(A->Type != SCALAR || (A+1)->Type != SCALAR)
179 Message::Error("Non scalar argument(s) for function Max");
180
181 V->Val[0] = (A->Val[0] > (A+1)->Val[0]) ? A->Val[0] : (A+1)->Val[0];
182 if (Current.NbrHar > 1){
183 V->Val[MAX_DIM] = 0. ;
184 for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2)
185 V->Val[MAX_DIM*k] = V->Val[MAX_DIM*(k+1)] = 0. ;
186 }
187 V->Type = SCALAR;
188 }
189
190 /* ------------------------------------------------------------------------ */
191 /* Bessel functions jn, yn and their derivatives */
192 /* ------------------------------------------------------------------------ */
193
F_Jn(F_ARG)194 void F_Jn(F_ARG)
195 {
196 int k, n;
197 double x;
198
199 if(A->Type != SCALAR || (A+1)->Type != SCALAR)
200 Message::Error("Non scalar argument(s) for Bessel function of the first kind 'Jn'");
201 n = (int)A->Val[0];
202 x = (A+1)->Val[0];
203
204 V->Val[0] = jn(n, x);
205
206 if (Current.NbrHar > 1){
207 V->Val[MAX_DIM] = 0. ;
208 for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2)
209 V->Val[MAX_DIM*k] = V->Val[MAX_DIM*(k+1)] = 0. ;
210 }
211 V->Type = SCALAR;
212 }
213
F_JnComplex(F_ARG)214 void F_JnComplex(F_ARG)
215 {
216 if(A->Type != SCALAR || (A+1)->Type != SCALAR)
217 Message::Error("Non scalar argument(s) for Bessel function of the first kind 'JnComplex'");
218 int n = (int)A->Val[0];
219 double xr = (A+1)->Val[0];
220 double xi = (A+1)->Val[MAX_DIM];
221 double valr, vali;
222
223 BesselJnComplex(n, 1, xr, xi, &valr, &vali);
224
225 V->Val[0] = valr;
226 V->Val[MAX_DIM] = vali;
227
228 if (Current.NbrHar > 1){
229 for (int k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2)
230 V->Val[MAX_DIM*k] = V->Val[MAX_DIM*(k+1)] = 0. ;
231 }
232 V->Type = SCALAR;
233 }
234
F_KnComplex(F_ARG)235 void F_KnComplex(F_ARG)
236 {
237 if(A->Type != SCALAR || (A+1)->Type != SCALAR)
238 Message::Error("Non scalar argument(s) for modified Bessel function of the second kind 'KnComplex'");
239 int n = (int)A->Val[0];
240 double xr = (A+1)->Val[0];
241 double xi = (A+1)->Val[MAX_DIM];
242 double valr, vali;
243
244 BesselKnComplex(n, 1, xr, xi, &valr, &vali);
245
246 V->Val[0] = valr;
247 V->Val[MAX_DIM] = vali;
248
249 if (Current.NbrHar > 1){
250 for (int k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2)
251 V->Val[MAX_DIM*k] = V->Val[MAX_DIM*(k+1)] = 0. ;
252 }
253 V->Type = SCALAR;
254 }
255
F_Yn(F_ARG)256 void F_Yn(F_ARG)
257 {
258 int k, n;
259 double x;
260
261 if(A->Type != SCALAR || (A+1)->Type != SCALAR)
262 Message::Error("Non scalar argument(s) for Bessel function of the second 'Yn'");
263 n = (int)A->Val[0];
264 x = (A+1)->Val[0];
265
266 V->Val[0] = yn(n, x);
267
268 if (Current.NbrHar > 1){
269 V->Val[MAX_DIM] = 0. ;
270 for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2)
271 V->Val[MAX_DIM*k] = V->Val[MAX_DIM*(k+1)] = 0. ;
272 }
273 V->Type = SCALAR;
274 }
275
dBessel(double * tab,int n,double x)276 double dBessel(double *tab, int n, double x)
277 {
278 if(n == 0){
279 return - tab[1];
280 }
281 else{
282 return tab[n-1] - (double)n/x * tab[n];
283 }
284 }
285
F_dJn(F_ARG)286 void F_dJn(F_ARG)
287 {
288
289 int k, n;
290 double x, *jntab;
291
292 if(A->Type != SCALAR || (A+1)->Type != SCALAR)
293 Message::Error("Non scalar argument(s) for the derivative of the Bessel "
294 "function of the first kind 'dJn'");
295 n = (int)A->Val[0];
296 x = (A+1)->Val[0];
297
298 jntab = (double*)Malloc((n + 2) * sizeof(double));
299 for(k = 0; k < n + 2; k++){
300 jntab[k] = jn(k, x);
301 }
302 V->Val[0] = dBessel(jntab, n, x);
303 Free(jntab);
304
305 if (Current.NbrHar > 1){
306 V->Val[MAX_DIM] = 0. ;
307 for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2)
308 V->Val[MAX_DIM*k] = V->Val[MAX_DIM*(k+1)] = 0. ;
309 }
310 V->Type = SCALAR;
311 }
312
F_dYn(F_ARG)313 void F_dYn(F_ARG)
314 {
315 int k, n;
316 double x, *yntab;
317
318 if(A->Type != SCALAR || (A+1)->Type != SCALAR)
319 Message::Error("Non scalar argument(s) for the derivative of the Bessel "
320 "function of the second kind 'dYn'");
321 n = (int)A->Val[0];
322 x = (A+1)->Val[0];
323
324 yntab = (double*)Malloc((n + 2) * sizeof(double));
325 for(k = 0; k < n + 2; k++){
326 yntab[k] = yn(k, x);
327 }
328 V->Val[0] = dBessel(yntab, n, x);
329 Free(yntab);
330
331 if (Current.NbrHar > 1){
332 V->Val[MAX_DIM] = 0. ;
333 for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2)
334 V->Val[MAX_DIM*k] = V->Val[MAX_DIM*(k+1)] = 0. ;
335 }
336 V->Type = SCALAR;
337 }
338
339 /* ------------------------------------------------------------------------ */
340 /* Spherical Bessel functions jn, yn and their derivatives */
341 /* ------------------------------------------------------------------------ */
342
F_JnSph(F_ARG)343 void F_JnSph(F_ARG)
344 {
345 if(A->Type != SCALAR || (A+1)->Type != SCALAR)
346 Message::Error("Non scalar argument(s) for function 'JnSph' (spherical Bessel function)");
347
348 int n = (int)A->Val[0];
349 double x = (A+1)->Val[0];
350
351 V->Type = SCALAR;
352 V->Val[0] = Spherical_j_n(n, x);
353 }
354
F_YnSph(F_ARG)355 void F_YnSph(F_ARG)
356 {
357 if(A->Type != SCALAR || (A+1)->Type != SCALAR)
358 Message::Error("Non scalar argument(s) for function 'YnSph' (spherical Bessel function)");
359
360 int n = (int)A->Val[0];
361 double x = (A+1)->Val[0];
362
363 V->Type = SCALAR;
364 V->Val[0] = Spherical_y_n(n, x);
365 }
366
F_dJnSph(F_ARG)367 void F_dJnSph(F_ARG)
368 {
369 if(A->Type != SCALAR || (A+1)->Type != SCALAR)
370 Message::Error("Non scalar argument(s) for function 'dJnSph' (derivative of spherical Bessel function)");
371
372 int n = (int)A->Val[0];
373 double x = (A+1)->Val[0];
374
375 V->Type = SCALAR;
376 V->Val[0] = (n/x) * Spherical_j_n(n, x) - Spherical_j_n(n+1, x);
377 }
378
F_dYnSph(F_ARG)379 void F_dYnSph(F_ARG)
380 {
381 if(A->Type != SCALAR || (A+1)->Type != SCALAR)
382 Message::Error("Non scalar argument(s) for function 'dYnSph' (derivative of spherical Bessel function)");
383
384 int n = (int)A->Val[0];
385 double x = (A+1)->Val[0];
386
387 V->Type = SCALAR;
388 V->Val[0] = (n/x) * Spherical_y_n(n, x) - Spherical_y_n(n+1, x);
389 }
390