1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 1997-2021 The Octave Project Developers
4 //
5 // See the file COPYRIGHT.md in the top-level directory of this
6 // distribution or <https://octave.org/copyright/>.
7 //
8 // This file is part of Octave.
9 //
10 // Octave is free software: you can redistribute it and/or modify it
11 // under the terms of the GNU General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // Octave is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU General Public License for more details.
19 //
20 // You should have received a copy of the GNU General Public License
21 // along with Octave; see the file COPYING. If not, see
22 // <https://www.gnu.org/licenses/>.
23 //
24 ////////////////////////////////////////////////////////////////////////
25
26 #if defined (HAVE_CONFIG_H)
27 # include "config.h"
28 #endif
29
30 #include "lo-specfun.h"
31 #include "quit.h"
32
33 #include "defun.h"
34 #include "error.h"
35 #include "ovl.h"
36 #include "utils.h"
37
38 enum bessel_type
39 {
40 BESSEL_J,
41 BESSEL_Y,
42 BESSEL_I,
43 BESSEL_K,
44 BESSEL_H1,
45 BESSEL_H2
46 };
47
48 #define DO_BESSEL(type, alpha, x, scaled, ierr, result) \
49 do \
50 { \
51 switch (type) \
52 { \
53 case BESSEL_J: \
54 result = octave::math::besselj (alpha, x, scaled, ierr); \
55 break; \
56 \
57 case BESSEL_Y: \
58 result = octave::math::bessely (alpha, x, scaled, ierr); \
59 break; \
60 \
61 case BESSEL_I: \
62 result = octave::math::besseli (alpha, x, scaled, ierr); \
63 break; \
64 \
65 case BESSEL_K: \
66 result = octave::math::besselk (alpha, x, scaled, ierr); \
67 break; \
68 \
69 case BESSEL_H1: \
70 result = octave::math::besselh1 (alpha, x, scaled, ierr); \
71 break; \
72 \
73 case BESSEL_H2: \
74 result = octave::math::besselh2 (alpha, x, scaled, ierr); \
75 break; \
76 \
77 default: \
78 break; \
79 } \
80 } \
81 while (0)
82
83 octave_value_list
do_bessel(enum bessel_type type,const char * fn,const octave_value_list & args,int nargout)84 do_bessel (enum bessel_type type, const char *fn,
85 const octave_value_list& args, int nargout)
86 {
87 int nargin = args.length ();
88
89 if (nargin < 2 || nargin > 3)
90 print_usage ();
91
92 octave_value_list retval (nargout > 1 ? 2 : 1);
93
94 bool scaled = false;
95 if (nargin == 3)
96 {
97 octave_value opt_arg = args(2);
98 bool rpt_error = false;
99
100 if (! opt_arg.is_scalar_type ())
101 rpt_error = true;
102 else if (opt_arg.isnumeric ())
103 {
104 double opt_val = opt_arg.double_value ();
105 if (opt_val != 0.0 && opt_val != 1.0)
106 rpt_error = true;
107 scaled = (opt_val == 1.0);
108 }
109 else if (opt_arg.islogical ())
110 scaled = opt_arg.bool_value ();
111
112 if (rpt_error)
113 error ("%s: OPT must be 0 (or false) or 1 (or true)", fn);
114 }
115
116 octave_value alpha_arg = args(0);
117 octave_value x_arg = args(1);
118
119 if (alpha_arg.is_single_type () || x_arg.is_single_type ())
120 {
121 if (alpha_arg.is_scalar_type ())
122 {
123 float alpha = args(0).xfloat_value ("%s: ALPHA must be a scalar or matrix", fn);
124
125 if (x_arg.is_scalar_type ())
126 {
127 FloatComplex x = x_arg.xfloat_complex_value ("%s: X must be a scalar or matrix", fn);
128
129 octave_idx_type ierr;
130 octave_value result;
131
132 DO_BESSEL (type, alpha, x, scaled, ierr, result);
133
134 retval(0) = result;
135 if (nargout > 1)
136 retval(1) = static_cast<float> (ierr);
137 }
138 else
139 {
140 FloatComplexNDArray x = x_arg.xfloat_complex_array_value ("%s: X must be a scalar or matrix", fn);
141
142 Array<octave_idx_type> ierr;
143 octave_value result;
144
145 DO_BESSEL (type, alpha, x, scaled, ierr, result);
146
147 retval(0) = result;
148 if (nargout > 1)
149 retval(1) = NDArray (ierr);
150 }
151 }
152 else
153 {
154 dim_vector dv0 = args(0).dims ();
155 dim_vector dv1 = args(1).dims ();
156
157 bool args0_is_row_vector = (dv0(1) == dv0.numel ());
158 bool args1_is_col_vector = (dv1(0) == dv1.numel ());
159
160 if (args0_is_row_vector && args1_is_col_vector)
161 {
162 FloatRowVector ralpha = args(0).xfloat_row_vector_value ("%s: ALPHA must be a scalar or matrix", fn);
163
164 FloatComplexColumnVector cx
165 = x_arg.xfloat_complex_column_vector_value ("%s: X must be a scalar or matrix", fn);
166
167 Array<octave_idx_type> ierr;
168 octave_value result;
169
170 DO_BESSEL (type, ralpha, cx, scaled, ierr, result);
171
172 retval(0) = result;
173 if (nargout > 1)
174 retval(1) = NDArray (ierr);
175 }
176 else
177 {
178 FloatNDArray alpha = args(0).xfloat_array_value ("%s: ALPHA must be a scalar or matrix", fn);
179
180 if (x_arg.is_scalar_type ())
181 {
182 FloatComplex x = x_arg.xfloat_complex_value ("%s: X must be a scalar or matrix", fn);
183
184 Array<octave_idx_type> ierr;
185 octave_value result;
186
187 DO_BESSEL (type, alpha, x, scaled, ierr, result);
188
189 retval(0) = result;
190 if (nargout > 1)
191 retval(1) = NDArray (ierr);
192 }
193 else
194 {
195 FloatComplexNDArray x = x_arg.xfloat_complex_array_value ("%s: X must be a scalar or matrix", fn);
196
197 Array<octave_idx_type> ierr;
198 octave_value result;
199
200 DO_BESSEL (type, alpha, x, scaled, ierr, result);
201
202 retval(0) = result;
203 if (nargout > 1)
204 retval(1) = NDArray (ierr);
205 }
206 }
207 }
208 }
209 else
210 {
211 if (alpha_arg.is_scalar_type ())
212 {
213 double alpha = args(0).xdouble_value ("%s: ALPHA must be a scalar or matrix", fn);
214
215 if (x_arg.is_scalar_type ())
216 {
217 Complex x = x_arg.xcomplex_value ("%s: X must be a scalar or matrix", fn);
218
219 octave_idx_type ierr;
220 octave_value result;
221
222 DO_BESSEL (type, alpha, x, scaled, ierr, result);
223
224 retval(0) = result;
225 if (nargout > 1)
226 retval(1) = static_cast<double> (ierr);
227 }
228 else
229 {
230 ComplexNDArray x = x_arg.xcomplex_array_value ("%s: X must be a scalar or matrix", fn);
231
232 Array<octave_idx_type> ierr;
233 octave_value result;
234
235 DO_BESSEL (type, alpha, x, scaled, ierr, result);
236
237 retval(0) = result;
238 if (nargout > 1)
239 retval(1) = NDArray (ierr);
240 }
241 }
242 else
243 {
244 dim_vector dv0 = args(0).dims ();
245 dim_vector dv1 = args(1).dims ();
246
247 bool args0_is_row_vector = (dv0(1) == dv0.numel ());
248 bool args1_is_col_vector = (dv1(0) == dv1.numel ());
249
250 if (args0_is_row_vector && args1_is_col_vector)
251 {
252 RowVector ralpha = args(0).xrow_vector_value ("%s: ALPHA must be a scalar or matrix", fn);
253
254 ComplexColumnVector cx
255 = x_arg.xcomplex_column_vector_value ("%s: X must be a scalar or matrix", fn);
256
257 Array<octave_idx_type> ierr;
258 octave_value result;
259
260 DO_BESSEL (type, ralpha, cx, scaled, ierr, result);
261
262 retval(0) = result;
263 if (nargout > 1)
264 retval(1) = NDArray (ierr);
265 }
266 else
267 {
268 NDArray alpha = args(0).xarray_value ("%s: ALPHA must be a scalar or matrix", fn);
269
270 if (x_arg.is_scalar_type ())
271 {
272 Complex x = x_arg.xcomplex_value ("%s: X must be a scalar or matrix", fn);
273
274 Array<octave_idx_type> ierr;
275 octave_value result;
276
277 DO_BESSEL (type, alpha, x, scaled, ierr, result);
278
279 retval(0) = result;
280 if (nargout > 1)
281 retval(1) = NDArray (ierr);
282 }
283 else
284 {
285 ComplexNDArray x = x_arg.xcomplex_array_value ("%s: X must be a scalar or matrix", fn);
286
287 Array<octave_idx_type> ierr;
288 octave_value result;
289
290 DO_BESSEL (type, alpha, x, scaled, ierr, result);
291
292 retval(0) = result;
293 if (nargout > 1)
294 retval(1) = NDArray (ierr);
295 }
296 }
297 }
298 }
299
300 return retval;
301 }
302
303 DEFUN (besselj, args, nargout,
304 doc: /* -*- texinfo -*-
305 @deftypefn {} {@var{J} =} besselj (@var{alpha}, @var{x})
306 @deftypefnx {} {@var{J} =} besselj (@var{alpha}, @var{x}, @var{opt})
307 @deftypefnx {} {[@var{J}, @var{ierr}] =} besselj (@dots{})
308 Compute Bessel functions of the first kind.
309
310 The order of the Bessel function @var{alpha} must be real. The points for
311 evaluation @var{x} may be complex.
312
313 If the optional argument @var{opt} is 1 or true, the result @var{J} is
314 multiplied by @w{@code{exp (-abs (imag (@var{x})))}}.
315
316 If @var{alpha} is a scalar, the result is the same size as @var{x}. If @var{x}
317 is a scalar, the result is the same size as @var{alpha}. If @var{alpha} is a
318 row vector and @var{x} is a column vector, the result is a matrix with
319 @code{length (@var{x})} rows and @code{length (@var{alpha})} columns.
320 Otherwise, @var{alpha} and @var{x} must conform and the result will be the same
321 size.
322
323 If requested, @var{ierr} contains the following status information and is the
324 same size as the result.
325
326 @enumerate 0
327 @item
328 Normal return.
329
330 @item
331 Input error, return @code{NaN}.
332
333 @item
334 Overflow, return @code{Inf}.
335
336 @item
337 Loss of significance by argument reduction results in less than half of machine
338 accuracy.
339
340 @item
341 Loss of significance by argument reduction, output may be inaccurate.
342
343 @item
344 Error---no computation, algorithm termination condition not met, return
345 @code{NaN}.
346 @end enumerate
347
348 @seealso{bessely, besseli, besselk, besselh}
349 @end deftypefn */)
350 {
351 return do_bessel (BESSEL_J, "besselj", args, nargout);
352 }
353
354 /*
355 %!# Function besselj is tested along with other bessels at the end of this file
356 */
357
358 DEFUN (bessely, args, nargout,
359 doc: /* -*- texinfo -*-
360 @deftypefn {} {@var{Y} =} bessely (@var{alpha}, @var{x})
361 @deftypefnx {} {@var{Y} =} bessely (@var{alpha}, @var{x}, @var{opt})
362 @deftypefnx {} {[@var{Y}, @var{ierr}] =} bessely (@dots{})
363 Compute Bessel functions of the second kind.
364
365 The order of the Bessel function @var{alpha} must be real. The points for
366 evaluation @var{x} may be complex.
367
368 If the optional argument @var{opt} is 1 or true, the result @var{Y} is
369 multiplied by @w{@code{exp (-abs (imag (@var{x})))}}.
370
371 If @var{alpha} is a scalar, the result is the same size as @var{x}. If @var{x}
372 is a scalar, the result is the same size as @var{alpha}. If @var{alpha} is a
373 row vector and @var{x} is a column vector, the result is a matrix with
374 @code{length (@var{x})} rows and @code{length (@var{alpha})} columns.
375 Otherwise, @var{alpha} and @var{x} must conform and the result will be the same
376 size.
377
378 If requested, @var{ierr} contains the following status information and is the
379 same size as the result.
380
381 @enumerate 0
382 @item
383 Normal return.
384
385 @item
386 Input error, return @code{NaN}.
387
388 @item
389 Overflow, return @code{Inf}.
390
391 @item
392 Loss of significance by argument reduction results in less than half of machine
393 accuracy.
394
395 @item
396 Complete loss of significance by argument reduction, return @code{NaN}.
397
398 @item
399 Error---no computation, algorithm termination condition not met, return
400 @code{NaN}.
401 @end enumerate
402
403 @seealso{besselj, besseli, besselk, besselh}
404 @end deftypefn */)
405 {
406 return do_bessel (BESSEL_Y, "bessely", args, nargout);
407 }
408
409 /*
410 %!# Function bessely is tested along with other bessels at the end of this file
411 */
412
413 DEFUN (besseli, args, nargout,
414 doc: /* -*- texinfo -*-
415 @deftypefn {} {@var{I} =} besseli (@var{alpha}, @var{x})
416 @deftypefnx {} {@var{I} =} besseli (@var{alpha}, @var{x}, @var{opt})
417 @deftypefnx {} {[@var{I}, @var{ierr}] =} besseli (@dots{})
418 Compute modified Bessel functions of the first kind.
419
420 The order of the Bessel function @var{alpha} must be real. The points for
421 evaluation @var{x} may be complex.
422
423 If the optional argument @var{opt} is 1 or true, the result @var{I} is
424 multiplied by @w{@code{exp (-abs (real (@var{x})))}}.
425
426 If @var{alpha} is a scalar, the result is the same size as @var{x}. If @var{x}
427 is a scalar, the result is the same size as @var{alpha}. If @var{alpha} is a
428 row vector and @var{x} is a column vector, the result is a matrix with
429 @code{length (@var{x})} rows and @code{length (@var{alpha})} columns.
430 Otherwise, @var{alpha} and @var{x} must conform and the result will be the same
431 size.
432
433 If requested, @var{ierr} contains the following status information and is the
434 same size as the result.
435
436 @enumerate 0
437 @item
438 Normal return.
439
440 @item
441 Input error, return @code{NaN}.
442
443 @item
444 Overflow, return @code{Inf}.
445
446 @item
447 Loss of significance by argument reduction results in less than half of machine
448 accuracy.
449
450 @item
451 Complete loss of significance by argument reduction, return @code{NaN}.
452
453 @item
454 Error---no computation, algorithm termination condition not met, return
455 @code{NaN}.
456 @end enumerate
457
458 @seealso{besselk, besselj, bessely, besselh}
459 @end deftypefn */)
460
461 {
462 return do_bessel (BESSEL_I, "besseli", args, nargout);
463 }
464
465 /*
466 %!# Function besseli is tested along with other bessels at the end of this file
467 */
468
469 DEFUN (besselk, args, nargout,
470 doc: /* -*- texinfo -*-
471 @deftypefn {} {@var{K} =} besselk (@var{alpha}, @var{x})
472 @deftypefnx {} {@var{K} =} besselk (@var{alpha}, @var{x}, @var{opt})
473 @deftypefnx {} {[@var{K}, @var{ierr}] =} besselk (@dots{})
474
475 Compute modified Bessel functions of the second kind.
476
477 The order of the Bessel function @var{alpha} must be real. The points for
478 evaluation @var{x} may be complex.
479
480 If the optional argument @var{opt} is 1 or true, the result @var{K} is
481 multiplied by @w{@code{exp (@var{x})}}.
482
483 If @var{alpha} is a scalar, the result is the same size as @var{x}. If @var{x}
484 is a scalar, the result is the same size as @var{alpha}. If @var{alpha} is a
485 row vector and @var{x} is a column vector, the result is a matrix with
486 @code{length (@var{x})} rows and @code{length (@var{alpha})} columns.
487 Otherwise, @var{alpha} and @var{x} must conform and the result will be the same
488 size.
489
490 If requested, @var{ierr} contains the following status information and is the
491 same size as the result.
492
493 @enumerate 0
494 @item
495 Normal return.
496
497 @item
498 Input error, return @code{NaN}.
499
500 @item
501 Overflow, return @code{Inf}.
502
503 @item
504 Loss of significance by argument reduction results in less than half of machine
505 accuracy.
506
507 @item
508 Complete loss of significance by argument reduction, return @code{NaN}.
509
510 @item
511 Error---no computation, algorithm termination condition not met, return
512 @code{NaN}.
513 @end enumerate
514
515 @seealso{besseli, besselj, bessely, besselh}
516 @end deftypefn */)
517 {
518 return do_bessel (BESSEL_K, "besselk", args, nargout);
519 }
520
521 /*
522 %!# Function besselk is tested along with other bessels at the end of this file
523 */
524
525 DEFUN (besselh, args, nargout,
526 doc: /* -*- texinfo -*-
527 @deftypefn {} {@var{H} =} besselh (@var{alpha}, @var{x})
528 @deftypefnx {} {@var{H} =} besselh (@var{alpha}, @var{k}, @var{x})
529 @deftypefnx {} {@var{H} =} besselh (@var{alpha}, @var{k}, @var{x}, @var{opt})
530 @deftypefnx {} {[@var{H}, @var{ierr}] =} besselh (@dots{})
531 Compute Bessel functions of the third kind (Hankel functions).
532
533 The order of the Bessel function @var{alpha} must be real. The kind of Hankel
534 function is specified by @var{k} and may be either first (@var{k} = 1) or
535 second (@var{k} = 2). The default is Hankel functions of the first kind. The
536 points for evaluation @var{x} may be complex.
537
538 If the optional argument @var{opt} is 1 or true, the result is multiplied
539 by @code{exp (-I*@var{x})} for @var{k} = 1 or @code{exp (I*@var{x})} for
540 @var{k} = 2.
541
542 If @var{alpha} is a scalar, the result is the same size as @var{x}. If @var{x}
543 is a scalar, the result is the same size as @var{alpha}. If @var{alpha} is a
544 row vector and @var{x} is a column vector, the result is a matrix with
545 @code{length (@var{x})} rows and @code{length (@var{alpha})} columns.
546 Otherwise, @var{alpha} and @var{x} must conform and the result will be the same
547 size.
548
549 If requested, @var{ierr} contains the following status information and is the
550 same size as the result.
551
552 @enumerate 0
553 @item
554 Normal return.
555
556 @item
557 Input error, return @code{NaN}.
558
559 @item
560 Overflow, return @code{Inf}.
561
562 @item
563 Loss of significance by argument reduction results in less than half of machine
564 accuracy.
565
566 @item
567 Complete loss of significance by argument reduction, return @code{NaN}.
568
569 @item
570 Error---no computation, algorithm termination condition not met, return
571 @code{NaN}.
572 @end enumerate
573
574 @seealso{besselj, bessely, besseli, besselk}
575 @end deftypefn */)
576 {
577 int nargin = args.length ();
578
579 if (nargin < 2 || nargin > 4)
580 print_usage ();
581
582 octave_value_list retval;
583
584 if (nargin == 2)
585 {
586 retval = do_bessel (BESSEL_H1, "besselh", args, nargout);
587 }
588 else
589 {
590 octave_idx_type kind = args(1).xint_value ("besselh: invalid value of K");
591
592 octave_value_list tmp_args;
593
594 if (nargin == 4)
595 tmp_args(2) = args(3);
596
597 tmp_args(1) = args(2);
598 tmp_args(0) = args(0);
599
600 if (kind == 1)
601 retval = do_bessel (BESSEL_H1, "besselh", tmp_args, nargout);
602 else if (kind == 2)
603 retval = do_bessel (BESSEL_H2, "besselh", tmp_args, nargout);
604 else
605 error ("besselh: K must be 1 or 2");
606 }
607
608 return retval;
609 }
610
611 /*
612 %!# Function besselh is tested along with other bessels at the end of this file
613 */
614
615 DEFUN (airy, args, nargout,
616 doc: /* -*- texinfo -*-
617 @deftypefn {} {[@var{a}, @var{ierr}] =} airy (@var{k}, @var{z}, @var{opt})
618 Compute Airy functions of the first and second kind, and their derivatives.
619
620 @example
621 @group
622 K Function Scale factor (if "opt" is supplied)
623 --- -------- ---------------------------------------
624 0 Ai (Z) exp ((2/3) * Z * sqrt (Z))
625 1 dAi(Z)/dZ exp ((2/3) * Z * sqrt (Z))
626 2 Bi (Z) exp (-abs (real ((2/3) * Z * sqrt (Z))))
627 3 dBi(Z)/dZ exp (-abs (real ((2/3) * Z * sqrt (Z))))
628 @end group
629 @end example
630
631 The function call @code{airy (@var{z})} is equivalent to
632 @code{airy (0, @var{z})}.
633
634 The result is the same size as @var{z}.
635
636 If requested, @var{ierr} contains the following status information and
637 is the same size as the result.
638
639 @enumerate 0
640 @item
641 Normal return.
642
643 @item
644 Input error, return @code{NaN}.
645
646 @item
647 Overflow, return @code{Inf}.
648
649 @item
650 Loss of significance by argument reduction results in less than half
651 of machine accuracy.
652
653 @item
654 Loss of significance by argument reduction, output may be inaccurate.
655
656 @item
657 Error---no computation, algorithm termination condition not met,
658 return @code{NaN}.
659 @end enumerate
660 @end deftypefn */)
661 {
662 int nargin = args.length ();
663
664 if (nargin < 1 || nargin > 3)
665 print_usage ();
666
667 octave_value_list retval (nargout > 1 ? 2 : 1);
668
669 int kind = 0;
670 if (nargin > 1)
671 {
672 kind = args(0).xint_value ("airy: K must be an integer value");
673
674 if (kind < 0 || kind > 3)
675 error ("airy: K must be 0, 1, 2, or 3");
676 }
677
678 bool scale = (nargin == 3);
679
680 int idx = (nargin == 1 ? 0 : 1);
681
682 if (args(idx).is_single_type ())
683 {
684 FloatComplexNDArray z = args(idx).xfloat_complex_array_value ("airy: Z must be a complex matrix");
685
686 Array<octave_idx_type> ierr;
687 octave_value result;
688
689 if (kind > 1)
690 result = octave::math::biry (z, kind == 3, scale, ierr);
691 else
692 result = octave::math::airy (z, kind == 1, scale, ierr);
693
694 retval(0) = result;
695 if (nargout > 1)
696 retval(1) = NDArray (ierr);
697 }
698 else
699 {
700 ComplexNDArray z = args(idx).xcomplex_array_value ("airy: Z must be a complex matrix");
701
702 Array<octave_idx_type> ierr;
703 octave_value result;
704
705 if (kind > 1)
706 result = octave::math::biry (z, kind == 3, scale, ierr);
707 else
708 result = octave::math::airy (z, kind == 1, scale, ierr);
709
710 retval(0) = result;
711 if (nargout > 1)
712 retval(1) = NDArray (ierr);
713 }
714
715 return retval;
716 }
717
718 /*
719 FIXME: Function airy does not yet have BIST tests
720 */
721
722 /*
723 ## Test values computed with GP/PARI version 2.3.3
724 %!shared alpha, x, jx, yx, ix, kx, nix
725 %!
726 %! ## Bessel functions, even order, positive and negative x
727 %! alpha = 2; x = 1.25;
728 %! jx = 0.1710911312405234823613091417;
729 %! yx = -1.193199310178553861283790424;
730 %! ix = 0.2220184483766341752692212604;
731 %! kx = 0.9410016167388185767085460540;
732 %!
733 %!assert (besselj (alpha,x), jx, 100*eps)
734 %!assert (bessely (alpha,x), yx, 100*eps)
735 %!assert (besseli (alpha,x), ix, 100*eps)
736 %!assert (besselk (alpha,x), kx, 100*eps)
737 %!assert (besselh (alpha,1,x), jx + I*yx, 100*eps)
738 %!assert (besselh (alpha,2,x), jx - I*yx, 100*eps)
739 %!
740 %!assert (besselj (alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
741 %!assert (bessely (alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
742 %!assert (besseli (alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
743 %!assert (besselk (alpha,x,1), kx*exp(x), 100*eps)
744 %!assert (besselh (alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
745 %!assert (besselh (alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
746 %!
747 %!assert (besselj (-alpha,x), jx, 100*eps)
748 %!assert (bessely (-alpha,x), yx, 100*eps)
749 %!assert (besseli (-alpha,x), ix, 100*eps)
750 %!assert (besselk (-alpha,x), kx, 100*eps)
751 %!assert (besselh (-alpha,1,x), jx + I*yx, 100*eps)
752 %!assert (besselh (-alpha,2,x), jx - I*yx, 100*eps)
753 %!
754 %!assert (besselj (-alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
755 %!assert (bessely (-alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
756 %!assert (besseli (-alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
757 %!assert (besselk (-alpha,x,1), kx*exp(x), 100*eps)
758 %!assert (besselh (-alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
759 %!assert (besselh (-alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
760 %!
761 %! x *= -1;
762 %! yx = -1.193199310178553861283790424 + 0.3421822624810469647226182835*I;
763 %! kx = 0.9410016167388185767085460540 - 0.6974915263814386815610060884*I;
764 %!
765 %!assert (besselj (alpha,x), jx, 100*eps)
766 %!assert (bessely (alpha,x), yx, 100*eps)
767 %!assert (besseli (alpha,x), ix, 100*eps)
768 %!assert (besselk (alpha,x), kx, 100*eps)
769 %!assert (besselh (alpha,1,x), jx + I*yx, 100*eps)
770 %!assert (besselh (alpha,2,x), jx - I*yx, 100*eps)
771 %!
772 %!assert (besselj (alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
773 %!assert (bessely (alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
774 %!assert (besseli (alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
775 %!assert (besselk (alpha,x,1), kx*exp(x), 100*eps)
776 %!assert (besselh (alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
777 %!assert (besselh (alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
778 %!
779 %! ## Bessel functions, odd order, positive and negative x
780 %! alpha = 3; x = 2.5;
781 %! jx = 0.2166003910391135247666890035;
782 %! yx = -0.7560554967536709968379029772;
783 %! ix = 0.4743704087780355895548240179;
784 %! kx = 0.2682271463934492027663765197;
785 %!
786 %!assert (besselj (alpha,x), jx, 100*eps)
787 %!assert (bessely (alpha,x), yx, 100*eps)
788 %!assert (besseli (alpha,x), ix, 100*eps)
789 %!assert (besselk (alpha,x), kx, 100*eps)
790 %!assert (besselh (alpha,1,x), jx + I*yx, 100*eps)
791 %!assert (besselh (alpha,2,x), jx - I*yx, 100*eps)
792 %!
793 %!assert (besselj (alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
794 %!assert (bessely (alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
795 %!assert (besseli (alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
796 %!assert (besselk (alpha,x,1), kx*exp(x), 100*eps)
797 %!assert (besselh (alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
798 %!assert (besselh (alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
799 %!
800 %!assert (besselj (-alpha,x), -jx, 100*eps)
801 %!assert (bessely (-alpha,x), -yx, 100*eps)
802 %!assert (besseli (-alpha,x), ix, 100*eps)
803 %!assert (besselk (-alpha,x), kx, 100*eps)
804 %!assert (besselh (-alpha,1,x), -(jx + I*yx), 100*eps)
805 %!assert (besselh (-alpha,2,x), -(jx - I*yx), 100*eps)
806 %!
807 %!assert (besselj (-alpha,x,1), -jx*exp(-abs(imag(x))), 100*eps)
808 %!assert (bessely (-alpha,x,1), -yx*exp(-abs(imag(x))), 100*eps)
809 %!assert (besseli (-alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
810 %!assert (besselk (-alpha,x,1), kx*exp(x), 100*eps)
811 %!assert (besselh (-alpha,1,x,1), -(jx + I*yx)*exp(-I*x), 100*eps)
812 %!assert (besselh (-alpha,2,x,1), -(jx - I*yx)*exp(I*x), 100*eps)
813 %!
814 %! x *= -1;
815 %! jx = -jx;
816 %! yx = 0.7560554967536709968379029772 - 0.4332007820782270495333780070*I;
817 %! ix = -ix;
818 %! kx = -0.2682271463934492027663765197 - 1.490278591297463775542004240*I;
819 %!
820 %!assert (besselj (alpha,x), jx, 100*eps)
821 %!assert (bessely (alpha,x), yx, 100*eps)
822 %!assert (besseli (alpha,x), ix, 100*eps)
823 %!assert (besselk (alpha,x), kx, 100*eps)
824 %!assert (besselh (alpha,1,x), jx + I*yx, 100*eps)
825 %!assert (besselh (alpha,2,x), jx - I*yx, 100*eps)
826 %!
827 %!assert (besselj (alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
828 %!assert (bessely (alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
829 %!assert (besseli (alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
830 %!assert (besselk (alpha,x,1), kx*exp(x), 100*eps)
831 %!assert (besselh (alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
832 %!assert (besselh (alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
833 %!
834 %! ## Bessel functions, fractional order, positive and negative x
835 %!
836 %! alpha = 3.5; x = 2.75;
837 %! jx = 0.1691636439842384154644784389;
838 %! yx = -0.8301381935499356070267953387;
839 %! ix = 0.3930540878794826310979363668;
840 %! kx = 0.2844099013460621170288192503;
841 %!
842 %!assert (besselj (alpha,x), jx, 100*eps)
843 %!assert (bessely (alpha,x), yx, 100*eps)
844 %!assert (besseli (alpha,x), ix, 100*eps)
845 %!assert (besselk (alpha,x), kx, 100*eps)
846 %!assert (besselh (alpha,1,x), jx + I*yx, 100*eps)
847 %!assert (besselh (alpha,2,x), jx - I*yx, 100*eps)
848 %!
849 %!assert (besselj (alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
850 %!assert (bessely (alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
851 %!assert (besseli (alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
852 %!assert (besselk (alpha,x,1), kx*exp(x), 100*eps)
853 %!assert (besselh (alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
854 %!assert (besselh (alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
855 %!
856 %! nix = 0.2119931212254662995364461998;
857 %!
858 %!assert (besselj (-alpha,x), yx, 100*eps)
859 %!assert (bessely (-alpha,x), -jx, 100*eps)
860 %!assert (besseli (-alpha,x), nix, 100*eps)
861 %!assert (besselk (-alpha,x), kx, 100*eps)
862 %!assert (besselh (-alpha,1,x), -I*(jx + I*yx), 100*eps)
863 %!assert (besselh (-alpha,2,x), I*(jx - I*yx), 100*eps)
864 %!
865 %!assert (besselj (-alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
866 %!assert (bessely (-alpha,x,1), -jx*exp(-abs(imag(x))), 100*eps)
867 %!assert (besseli (-alpha,x,1), nix*exp(-abs(real(x))), 100*eps)
868 %!assert (besselk (-alpha,x,1), kx*exp(x), 100*eps)
869 %!assert (besselh (-alpha,1,x,1), -I*(jx + I*yx)*exp(-I*x), 100*eps)
870 %!assert (besselh (-alpha,2,x,1), I*(jx - I*yx)*exp(I*x), 100*eps)
871 %!
872 %! x *= -1;
873 %! jx *= -I;
874 %! yx = -0.8301381935499356070267953387*I;
875 %! ix *= -I;
876 %! kx = -0.9504059335995575096509874508*I;
877 %!
878 %!assert (besselj (alpha,x), jx, 100*eps)
879 %!assert (bessely (alpha,x), yx, 100*eps)
880 %!assert (besseli (alpha,x), ix, 100*eps)
881 %!assert (besselk (alpha,x), kx, 100*eps)
882 %!assert (besselh (alpha,1,x), jx + I*yx, 100*eps)
883 %!assert (besselh (alpha,2,x), jx - I*yx, 100*eps)
884 %!
885 %!assert (besselj (alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
886 %!assert (bessely (alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
887 %!assert (besseli (alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
888 %!assert (besselk (alpha,x,1), kx*exp(x), 100*eps)
889 %!assert (besselh (alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
890 %!assert (besselh (alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
891 %!
892 %! ## Bessel functions, even order, complex x
893 %!
894 %! alpha = 2; x = 1.25 + 3.625 * I;
895 %! jx = -1.299533366810794494030065917 + 4.370833116012278943267479589*I;
896 %! yx = -4.370357232383223896393056727 - 1.283083391453582032688834041*I;
897 %! ix = -0.6717801680341515541002273932 - 0.2314623443930774099910228553*I;
898 %! kx = -0.01108009888623253515463783379 + 0.2245218229358191588208084197*I;
899 %!
900 %!assert (besselj (alpha,x), jx, 100*eps)
901 %!assert (bessely (alpha,x), yx, 100*eps)
902 %!assert (besseli (alpha,x), ix, 100*eps)
903 %!assert (besselk (alpha,x), kx, 100*eps)
904 %!assert (besselh (alpha,1,x), jx + I*yx, 100*eps)
905 %!assert (besselh (alpha,2,x), jx - I*yx, 100*eps)
906 %!
907 %!assert (besselj (alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
908 %!assert (bessely (alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
909 %!assert (besseli (alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
910 %!assert (besselk (alpha,x,1), kx*exp(x), 100*eps)
911 %!assert (besselh (alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
912 %!assert (besselh (alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
913 %!
914 %!assert (besselj (-alpha,x), jx, 100*eps)
915 %!assert (bessely (-alpha,x), yx, 100*eps)
916 %!assert (besseli (-alpha,x), ix, 100*eps)
917 %!assert (besselk (-alpha,x), kx, 100*eps)
918 %!assert (besselh (-alpha,1,x), jx + I*yx, 100*eps)
919 %!assert (besselh (-alpha,2,x), jx - I*yx, 100*eps)
920 %!
921 %!assert (besselj (-alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
922 %!assert (bessely (-alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
923 %!assert (besseli (-alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
924 %!assert (besselk (-alpha,x,1), kx*exp(x), 100*eps)
925 %!assert (besselh (-alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
926 %!assert (besselh (-alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
927 %!
928 %! ## Bessel functions, odd order, complex x
929 %!
930 %! alpha = 3; x = 2.5 + 1.875 * I;
931 %! jx = 0.1330721523048277493333458596 + 0.5386295217249660078754395597*I;
932 %! yx = -0.6485072392105829901122401551 + 0.2608129289785456797046996987*I;
933 %! ix = -0.6182064685486998097516365709 + 0.4677561094683470065767989920*I;
934 %! kx = -0.1568585587733540007867882337 - 0.05185853709490846050505141321*I;
935 %!
936 %!assert (besselj (alpha,x), jx, 100*eps)
937 %!assert (bessely (alpha,x), yx, 100*eps)
938 %!assert (besseli (alpha,x), ix, 100*eps)
939 %!assert (besselk (alpha,x), kx, 100*eps)
940 %!assert (besselh (alpha,1,x), jx + I*yx, 100*eps)
941 %!assert (besselh (alpha,2,x), jx - I*yx, 100*eps)
942 %!
943 %!assert (besselj (alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
944 %!assert (bessely (alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
945 %!assert (besseli (alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
946 %!assert (besselk (alpha,x,1), kx*exp(x), 100*eps)
947 %!assert (besselh (alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
948 %!assert (besselh (alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
949 %!
950 %!assert (besselj (-alpha,x), -jx, 100*eps)
951 %!assert (bessely (-alpha,x), -yx, 100*eps)
952 %!assert (besseli (-alpha,x), ix, 100*eps)
953 %!assert (besselk (-alpha,x), kx, 100*eps)
954 %!assert (besselh (-alpha,1,x), -(jx + I*yx), 100*eps)
955 %!assert (besselh (-alpha,2,x), -(jx - I*yx), 100*eps)
956 %!
957 %!assert (besselj (-alpha,x,1), -jx*exp(-abs(imag(x))), 100*eps)
958 %!assert (bessely (-alpha,x,1), -yx*exp(-abs(imag(x))), 100*eps)
959 %!assert (besseli (-alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
960 %!assert (besselk (-alpha,x,1), kx*exp(x), 100*eps)
961 %!assert (besselh (-alpha,1,x,1), -(jx + I*yx)*exp(-I*x), 100*eps)
962 %!assert (besselh (-alpha,2,x,1), -(jx - I*yx)*exp(I*x), 100*eps)
963 %!
964 %! ## Bessel functions, fractional order, complex x
965 %!
966 %! alpha = 3.5; x = 1.75 + 4.125 * I;
967 %! jx = -3.018566131370455929707009100 - 0.7585648436793900607704057611*I;
968 %! yx = 0.7772278839106298215614791107 - 3.018518722313849782683792010*I;
969 %! ix = 0.2100873577220057189038160913 - 0.6551765604618246531254970926*I;
970 %! kx = 0.1757147290513239935341488069 + 0.08772348296883849205562558311*I;
971 %!
972 %!assert (besselj (alpha,x), jx, 100*eps)
973 %!assert (bessely (alpha,x), yx, 100*eps)
974 %!assert (besseli (alpha,x), ix, 100*eps)
975 %!assert (besselk (alpha,x), kx, 100*eps)
976 %!assert (besselh (alpha,1,x), jx + I*yx, 100*eps)
977 %!assert (besselh (alpha,2,x), jx - I*yx, 100*eps)
978 %!
979 %!assert (besselj (alpha,x,1), jx*exp(-abs(imag(x))), 100*eps)
980 %!assert (bessely (alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
981 %!assert (besseli (alpha,x,1), ix*exp(-abs(real(x))), 100*eps)
982 %!assert (besselk (alpha,x,1), kx*exp(x), 100*eps)
983 %!assert (besselh (alpha,1,x,1), (jx + I*yx)*exp(-I*x), 100*eps)
984 %!assert (besselh (alpha,2,x,1), (jx - I*yx)*exp(I*x), 100*eps)
985 %!
986 %! nix = 0.09822388691172060573913739253 - 0.7110230642207380127317227407*I;
987 %!
988 %!assert (besselj (-alpha,x), yx, 100*eps)
989 %!assert (bessely (-alpha,x), -jx, 100*eps)
990 %!assert (besseli (-alpha,x), nix, 100*eps)
991 %!assert (besselk (-alpha,x), kx, 100*eps)
992 %!assert (besselh (-alpha,1,x), -I*(jx + I*yx), 100*eps)
993 %!assert (besselh (-alpha,2,x), I*(jx - I*yx), 100*eps)
994 %!
995 %!assert (besselj (-alpha,x,1), yx*exp(-abs(imag(x))), 100*eps)
996 %!assert (bessely (-alpha,x,1), -jx*exp(-abs(imag(x))), 100*eps)
997 %!assert (besseli (-alpha,x,1), nix*exp(-abs(real(x))), 100*eps)
998 %!assert (besselk (-alpha,x,1), kx*exp(x), 100*eps)
999 %!assert (besselh (-alpha,1,x,1), -I*(jx + I*yx)*exp(-I*x), 100*eps)
1000 %!assert (besselh (-alpha,2,x,1), I*(jx - I*yx)*exp(I*x), 100*eps)
1001
1002 Tests contributed by Robert T. Short.
1003 Tests are based on the properties and tables in A&S:
1004 Abramowitz and Stegun, "Handbook of Mathematical Functions",
1005 1972.
1006
1007 For regular Bessel functions, there are 3 tests. These compare octave
1008 results against Tables 9.1, 9.2, and 9.4 in A&S. Tables 9.1 and 9.2
1009 are good to only a few decimal places, so any failures should be
1010 considered a broken implementation. Table 9.4 is an extended table
1011 for larger orders and arguments. There are some differences between
1012 Octave and Table 9.4, mostly in the last decimal place but in a very
1013 few instances the errors are in the last two places. The comparison
1014 tolerance has been changed to reflect this.
1015
1016 Similarly for modified Bessel functions, there are 3 tests. These
1017 compare octave results against Tables 9.8, 9.9, and 9.11 in A&S.
1018 Tables 9.8 and 9.9 are good to only a few decimal places, so any
1019 failures should be considered a broken implementation. Table 9.11 is
1020 an extended table for larger orders and arguments. There are some
1021 differences between octave and Table 9.11, mostly in the last decimal
1022 place but in a very few instances the errors are in the last two
1023 places. The comparison tolerance has been changed to reflect this.
1024
1025 For spherical Bessel functions, there are also three tests, comparing
1026 octave results to Tables 10.1, 10.2, and 10.4 in A&S. Very similar
1027 comments may be made here as in the previous lines. At this time,
1028 modified spherical Bessel function tests are not included.
1029
1030 % Table 9.1 - J and Y for integer orders 0, 1, 2.
1031 % Compare against excerpts of Table 9.1, Abramowitz and Stegun.
1032 %!test
1033 %! n = 0:2;
1034 %! z = (0:2.5:17.5)';
1035 %!
1036 %! Jt = [[ 1.000000000000000, 0.0000000000, 0.0000000000];
1037 %! [-0.048383776468198, 0.4970941025, 0.4460590584];
1038 %! [-0.177596771314338, -0.3275791376, 0.0465651163];
1039 %! [ 0.266339657880378, 0.1352484276, -0.2302734105];
1040 %! [-0.245935764451348, 0.0434727462, 0.2546303137];
1041 %! [ 0.146884054700421, -0.1654838046, -0.1733614634];
1042 %! [-0.014224472826781, 0.2051040386, 0.0415716780];
1043 %! [-0.103110398228686, -0.1634199694, 0.0844338303]];
1044 %!
1045 %! Yt = [[-Inf, -Inf, -Inf ];
1046 %! [ 0.4980703596, 0.1459181380, -0.38133585 ];
1047 %! [-0.3085176252, 0.1478631434, 0.36766288 ];
1048 %! [ 0.1173132861, -0.2591285105, -0.18641422 ];
1049 %! [ 0.0556711673, 0.2490154242, -0.00586808 ];
1050 %! [-0.1712143068, -0.1538382565, 0.14660019 ];
1051 %! [ 0.2054642960, 0.0210736280, -0.20265448 ];
1052 %! [-0.1604111925, 0.0985727987, 0.17167666 ]];
1053 %!
1054 %! J = besselj (n,z);
1055 %! Y = bessely (n,z);
1056 %! assert (Jt(:,1), J(:,1), 0.5e-10);
1057 %! assert (Yt(:,1), Y(:,1), 0.5e-10);
1058 %! assert (Jt(:,2:3), J(:,2:3), 0.5e-10);
1059
1060 Table 9.2 - J and Y for integer orders 3-9.
1061
1062 %!test
1063 %! n = (3:9);
1064 %! z = (0:2:20).';
1065 %!
1066 %! Jt = [[ 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00];
1067 %! [ 1.2894e-01, 3.3996e-02, 7.0396e-03, 1.2024e-03, 1.7494e-04, 2.2180e-05, 2.4923e-06];
1068 %! [ 4.3017e-01, 2.8113e-01, 1.3209e-01, 4.9088e-02, 1.5176e-02, 4.0287e-03, 9.3860e-04];
1069 %! [ 1.1477e-01, 3.5764e-01, 3.6209e-01, 2.4584e-01, 1.2959e-01, 5.6532e-02, 2.1165e-02];
1070 %! [-2.9113e-01,-1.0536e-01, 1.8577e-01, 3.3758e-01, 3.2059e-01, 2.2345e-01, 1.2632e-01];
1071 %! [ 5.8379e-02,-2.1960e-01,-2.3406e-01,-1.4459e-02, 2.1671e-01, 3.1785e-01, 2.9186e-01];
1072 %! [ 1.9514e-01, 1.8250e-01,-7.3471e-02,-2.4372e-01,-1.7025e-01, 4.5095e-02, 2.3038e-01];
1073 %! [-1.7681e-01, 7.6244e-02, 2.2038e-01, 8.1168e-02,-1.5080e-01,-2.3197e-01,-1.1431e-01];
1074 %! [-4.3847e-02,-2.0264e-01,-5.7473e-02, 1.6672e-01, 1.8251e-01,-7.0211e-03,-1.8953e-01];
1075 %! [ 1.8632e-01, 6.9640e-02,-1.5537e-01,-1.5596e-01, 5.1399e-02, 1.9593e-01, 1.2276e-01];
1076 %! [-9.8901e-02, 1.3067e-01, 1.5117e-01,-5.5086e-02,-1.8422e-01,-7.3869e-02, 1.2513e-01]];
1077 %!
1078 %! Yt = [[ -Inf, -Inf, -Inf, -Inf, -Inf, -Inf, -Inf];
1079 %! [-1.1278e+00,-2.7659e+00,-9.9360e+00,-4.6914e+01,-2.7155e+02,-1.8539e+03,-1.4560e+04];
1080 %! [-1.8202e-01,-4.8894e-01,-7.9585e-01,-1.5007e+00,-3.7062e+00,-1.1471e+01,-4.2178e+01];
1081 %! [ 3.2825e-01, 9.8391e-02,-1.9706e-01,-4.2683e-01,-6.5659e-01,-1.1052e+00,-2.2907e+00];
1082 %! [ 2.6542e-02, 2.8294e-01, 2.5640e-01, 3.7558e-02,-2.0006e-01,-3.8767e-01,-5.7528e-01];
1083 %! [-2.5136e-01,-1.4495e-01, 1.3540e-01, 2.8035e-01, 2.0102e-01, 1.0755e-03,-1.9930e-01];
1084 %! [ 1.2901e-01,-1.5122e-01,-2.2982e-01,-4.0297e-02, 1.8952e-01, 2.6140e-01, 1.5902e-01];
1085 %! [ 1.2350e-01, 2.0393e-01,-6.9717e-03,-2.0891e-01,-1.7209e-01, 3.6816e-02, 2.1417e-01];
1086 %! [-1.9637e-01,-7.3222e-05, 1.9633e-01, 1.2278e-01,-1.0425e-01,-2.1399e-01,-1.0975e-01];
1087 %! [ 3.3724e-02,-1.7722e-01,-1.1249e-01, 1.1472e-01, 1.8897e-01, 3.2253e-02,-1.6030e-01];
1088 %! [ 1.4967e-01, 1.2409e-01,-1.0004e-01,-1.7411e-01,-4.4312e-03, 1.7101e-01, 1.4124e-01]];
1089 %!
1090 %! n = (3:9);
1091 %! z = (0:2:20).';
1092 %! J = besselj (n,z);
1093 %! Y = bessely (n,z);
1094 %!
1095 %! assert (J(1,:), zeros (1, columns (J)));
1096 %! assert (J(2:end,:), Jt(2:end,:), -5e-5);
1097 %! assert (Yt(1,:), Y(1,:));
1098 %! assert (Y(2:end,:), Yt(2:end,:), -5e-5);
1099
1100 Table 9.4 - J and Y for various integer orders and arguments.
1101
1102 %!test
1103 %! Jt = [[ 7.651976866e-01, 2.238907791e-01, -1.775967713e-01, -2.459357645e-01, 5.581232767e-02, 1.998585030e-02];
1104 %! [ 2.497577302e-04, 7.039629756e-03, 2.611405461e-01, -2.340615282e-01, -8.140024770e-02, -7.419573696e-02];
1105 %! [ 2.630615124e-10, 2.515386283e-07, 1.467802647e-03, 2.074861066e-01, -1.138478491e-01, -5.473217694e-02];
1106 %! [ 2.297531532e-17, 7.183016356e-13, 4.796743278e-07, 4.507973144e-03, -1.082255990e-01, 1.519812122e-02];
1107 %! [ 3.873503009e-25, 3.918972805e-19, 2.770330052e-11, 1.151336925e-05, -1.167043528e-01, 6.221745850e-02];
1108 %! [ 3.482869794e-42, 3.650256266e-33, 2.671177278e-21, 1.551096078e-12, 4.843425725e-02, 8.146012958e-02];
1109 %! [ 1.107915851e-60, 1.196077458e-48, 8.702241617e-33, 6.030895312e-21, -1.381762812e-01, 7.270175482e-02];
1110 %! [ 2.906004948e-80, 3.224095839e-65, 2.294247616e-45, 1.784513608e-30, 1.214090219e-01, -3.869833973e-02];
1111 %! [ 8.431828790e-189, 1.060953112e-158, 6.267789396e-119, 6.597316064e-89, 1.115927368e-21, 9.636667330e-02]];
1112 %!
1113 %! Yt = [[ 8.825696420e-02, 5.103756726e-01, -3.085176252e-01, 5.567116730e-02, -9.806499547e-02, -7.724431337e-02]
1114 %! [-2.604058666e+02, -9.935989128e+00, -4.536948225e-01, 1.354030477e-01, -7.854841391e-02, -2.948019628e-02]
1115 %! [-1.216180143e+08, -1.291845422e+05, -2.512911010e+01, -3.598141522e-01, 5.723897182e-03, 5.833157424e-02]
1116 %! [-9.256973276e+14, -2.981023646e+10, -4.694049564e+04, -6.364745877e+00, 4.041280205e-02, 7.879068695e-02]
1117 %! [-4.113970315e+22, -4.081651389e+16, -5.933965297e+08, -1.597483848e+03, 1.644263395e-02, 5.124797308e-02]
1118 %! [-3.048128783e+39, -2.913223848e+30, -4.028568418e+18, -7.256142316e+09, -1.164572349e-01, 6.138839212e-03]
1119 %! [-7.184874797e+57, -6.661541235e+45, -9.216816571e+29, -1.362803297e+18, -4.530801120e-02, 4.074685217e-02]
1120 %! [-2.191142813e+77, -1.976150576e+62, -2.788837017e+42, -3.641066502e+27, -2.103165546e-01, 7.650526394e-02]
1121 %! [-3.775287810e+185, -3.000826049e+155, -5.084863915e+115, -4.849148271e+85, -3.293800188e+18, -1.669214114e-01]];
1122 %!
1123 %! n = [(0:5:20).';30;40;50;100];
1124 %! z = [1,2,5,10,50,100];
1125 %! J = besselj (n.', z.').';
1126 %! Y = bessely (n.', z.').';
1127 %! assert (J, Jt, -1e-9);
1128 %! assert (Y, Yt, -1e-9);
1129
1130 Table 9.8 - I and K for integer orders 0, 1, 2.
1131
1132 %!test
1133 %! n = 0:2;
1134 %! z1 = [0.1;2.5;5.0];
1135 %! z2 = [7.5;10.0;15.0;20.0];
1136 %! rtbl = [[ 0.9071009258 0.0452984468 0.1251041992 2.6823261023 10.890182683 1.995039646 ];
1137 %! [ 0.2700464416 0.2065846495 0.2042345837 0.7595486903 0.9001744239 0.759126289 ];
1138 %! [ 0.1835408126 0.1639722669 0.7002245988 0.5478075643 0.6002738588 0.132723593 ];
1139 %! [ 0.1483158301 0.1380412115 0.111504840 0.4505236991 0.4796689336 0.57843541 ];
1140 %! [ 0.1278333372 0.1212626814 0.103580801 0.3916319344 0.4107665704 0.47378525 ];
1141 %! [ 0.1038995314 0.1003741751 0.090516308 0.3210023535 0.3315348950 0.36520701 ];
1142 %! [ 0.0897803119 0.0875062222 0.081029690 0.2785448768 0.2854254970 0.30708743 ]];
1143 %!
1144 %! tbl = [besseli(n,z1,1), besselk(n,z1,1)];
1145 %! tbl(:,3) = tbl(:,3) .* (exp (z1) .* z1.^(-2));
1146 %! tbl(:,6) = tbl(:,6) .* (exp (-z1) .* z1.^(2));
1147 %! tbl = [tbl;[besseli(n,z2,1),besselk(n,z2,1)]];
1148 %!
1149 %! assert (tbl, rtbl, -2e-8);
1150
1151 Table 9.9 - I and K for orders 3-9.
1152
1153 %!test
1154 %! It = [[ 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00];
1155 %! [ 2.8791e-02 6.8654e-03 1.3298e-03 2.1656e-04 3.0402e-05 3.7487e-06 4.1199e-07];
1156 %! [ 6.1124e-02 2.5940e-02 9.2443e-03 2.8291e-03 7.5698e-04 1.7968e-04 3.8284e-05];
1157 %! [ 7.4736e-02 4.1238e-02 1.9752e-02 8.3181e-03 3.1156e-03 1.0484e-03 3.1978e-04];
1158 %! [ 7.9194e-02 5.0500e-02 2.8694e-02 1.4633e-02 6.7449e-03 2.8292e-03 1.0866e-03];
1159 %! [ 7.9830e-02 5.5683e-02 3.5284e-02 2.0398e-02 1.0806e-02 5.2694e-03 2.3753e-03];
1160 %! [ 7.8848e-02 5.8425e-02 3.9898e-02 2.5176e-02 1.4722e-02 8.0010e-03 4.0537e-03];
1161 %! [ 7.7183e-02 5.9723e-02 4.3056e-02 2.8969e-02 1.8225e-02 1.0744e-02 5.9469e-03];
1162 %! [ 7.5256e-02 6.0155e-02 4.5179e-02 3.1918e-02 2.1240e-02 1.3333e-02 7.9071e-03];
1163 %! [ 7.3263e-02 6.0059e-02 4.6571e-02 3.4186e-02 2.3780e-02 1.5691e-02 9.8324e-03];
1164 %! [ 7.1300e-02 5.9640e-02 4.7444e-02 3.5917e-02 2.5894e-02 1.7792e-02 1.1661e-02]];
1165 %!
1166 %! Kt = [[ Inf Inf Inf Inf Inf Inf Inf];
1167 %! [ 4.7836e+00 1.6226e+01 6.9687e+01 3.6466e+02 2.2576e+03 1.6168e+04 1.3160e+05];
1168 %! [ 1.6317e+00 3.3976e+00 8.4268e+00 2.4465e+01 8.1821e+01 3.1084e+02 1.3252e+03];
1169 %! [ 9.9723e-01 1.6798e+00 3.2370e+00 7.0748e+00 1.7387e+01 4.7644e+01 1.4444e+02];
1170 %! [ 7.3935e-01 1.1069e+00 1.8463e+00 3.4148e+00 6.9684e+00 1.5610e+01 3.8188e+01];
1171 %! [ 6.0028e-01 8.3395e-01 1.2674e+00 2.1014e+00 3.7891e+00 7.4062e+00 1.5639e+01];
1172 %! [ 5.1294e-01 6.7680e-01 9.6415e-01 1.4803e+00 2.4444e+00 4.3321e+00 8.2205e+00];
1173 %! [ 4.5266e-01 5.7519e-01 7.8133e-01 1.1333e+00 1.7527e+00 2.8860e+00 5.0510e+00];
1174 %! [ 4.0829e-01 5.0414e-01 6.6036e-01 9.1686e-01 1.3480e+00 2.0964e+00 3.4444e+00];
1175 %! [ 3.7411e-01 4.5162e-01 5.7483e-01 7.7097e-01 1.0888e+00 1.6178e+00 2.5269e+00];
1176 %! [ 3.4684e-01 4.1114e-01 5.1130e-01 6.6679e-01 9.1137e-01 1.3048e+00 1.9552e+00]];
1177 %!
1178 %! n = (3:9);
1179 %! z = (0:2:20).';
1180 %! I = besseli (n,z,1);
1181 %! K = besselk (n,z,1);
1182 %!
1183 %! assert (abs (I(1,:)), zeros (1, columns (I)));
1184 %! assert (I(2:end,:), It(2:end,:), -5e-5);
1185 %! assert (Kt(1,:), K(1,:));
1186 %! assert (K(2:end,:), Kt(2:end,:), -5e-5);
1187
1188 Table 9.11 - I and K for various integer orders and arguments.
1189
1190 %!test
1191 %! It = [[ 1.266065878e+00 2.279585302e+00 2.723987182e+01 2.815716628e+03 2.93255378e+20 1.07375171e+42 ];
1192 %! [ 2.714631560e-04 9.825679323e-03 2.157974547e+00 7.771882864e+02 2.27854831e+20 9.47009387e+41 ];
1193 %! [ 2.752948040e-10 3.016963879e-07 4.580044419e-03 2.189170616e+01 1.07159716e+20 6.49897552e+41 ];
1194 %! [ 2.370463051e-17 8.139432531e-13 1.047977675e-06 1.043714907e-01 3.07376455e+19 3.47368638e+41 ];
1195 %! [ 3.966835986e-25 4.310560576e-19 5.024239358e-11 1.250799736e-04 5.44200840e+18 1.44834613e+41 ];
1196 %! [ 3.539500588e-42 3.893519664e-33 3.997844971e-21 7.787569783e-12 4.27499365e+16 1.20615487e+40 ];
1197 %! [ 1.121509741e-60 1.255869192e-48 1.180426980e-32 2.042123274e-20 6.00717897e+13 3.84170550e+38 ];
1198 %! [ 2.934635309e-80 3.353042830e-65 2.931469647e-45 4.756894561e-30 1.76508024e+10 4.82195809e+36 ];
1199 %! [ 8.473674008e-189 1.082171475e-158 7.093551489e-119 1.082344202e-88 2.72788795e-16 4.64153494e+21 ]];
1200 %!
1201 %! Kt = [[ 4.210244382e-01 1.138938727e-01 3.691098334e-03 1.778006232e-05 3.41016774e-23 4.65662823e-45 ];
1202 %! [ 3.609605896e+02 9.431049101e+00 3.270627371e-02 5.754184999e-05 4.36718224e-23 5.27325611e-45 ];
1203 %! [ 1.807132899e+08 1.624824040e+05 9.758562829e+00 1.614255300e-03 9.15098819e-23 7.65542797e-45 ];
1204 %! [ 1.403066801e+15 4.059213332e+10 3.016976630e+04 2.656563849e-01 3.11621117e-22 1.42348325e-44 ];
1205 %! [ 6.294369360e+22 5.770856853e+16 4.827000521e+08 1.787442782e+02 1.70614838e-21 3.38520541e-44 ];
1206 %! [ 4.706145527e+39 4.271125755e+30 4.112132063e+18 2.030247813e+09 2.00581681e-19 3.97060205e-43 ];
1207 %! [ 1.114220651e+58 9.940839886e+45 1.050756722e+30 5.938224681e+17 1.29986971e-16 1.20842080e-41 ];
1208 %! [ 3.406896854e+77 2.979981740e+62 3.394322243e+42 2.061373775e+27 4.00601349e-13 9.27452265e-40 ];
1209 %! [ 5.900333184e+185 4.619415978e+155 7.039860193e+115 4.596674084e+85 1.63940352e+13 7.61712963e-25 ]];
1210 %!
1211 %! n = [(0:5:20).';30;40;50;100];
1212 %! z = [1,2,5,10,50,100];
1213 %! I = besseli (n.', z.').';
1214 %! K = besselk (n.', z.').';
1215 %! assert (I, It, -5e-9);
1216 %! assert (K, Kt, -5e-9);
1217
1218 The next section checks that negative integer orders and positive
1219 integer orders are appropriately related.
1220
1221 %!test
1222 %! n = (0:2:20);
1223 %! assert (besselj (n,1), besselj (-n,1), 1e-8);
1224 %! assert (-besselj (n+1,1), besselj (-n-1,1), 1e-8);
1225
1226 besseli (n,z) = besseli (-n,z);
1227
1228 %!test
1229 %! n = (0:2:20);
1230 %! assert (besseli (n,1), besseli (-n,1), 1e-8);
1231
1232 Table 10.1 - j and y for integer orders 0, 1, 2.
1233 Compare against excerpts of Table 10.1, Abramowitz and Stegun.
1234
1235 %!test
1236 %! n = (0:2);
1237 %! z = [0.1;(2.5:2.5:10.0).'];
1238 %!
1239 %! jt = [[ 9.9833417e-01 3.33000119e-02 6.6619061e-04 ];
1240 %! [ 2.3938886e-01 4.16212989e-01 2.6006673e-01 ];
1241 %! [-1.9178485e-01 -9.50894081e-02 1.3473121e-01 ];
1242 %! [ 1.2507e-01 -2.9542e-02 -1.3688e-01 ];
1243 %! [ -5.4402e-02 7.8467e-02 7.7942e-02 ]];
1244 %!
1245 %! yt = [[-9.9500417e+00 -1.0049875e+02 -3.0050125e+03 ];
1246 %! [ 3.2045745e-01 -1.1120588e-01 -4.5390450e-01 ];
1247 %! [-5.6732437e-02 1.8043837e-01 1.6499546e-01 ];
1248 %! [ -4.6218e-02 -1.3123e-01 -6.2736e-03 ];
1249 %! [ 8.3907e-02 6.2793e-02 -6.5069e-02 ]];
1250 %!
1251 %! j = sqrt ((pi/2)./z) .* besselj (n+1/2,z);
1252 %! y = sqrt ((pi/2)./z) .* bessely (n+1/2,z);
1253 %! assert (jt, j, -5e-5);
1254 %! assert (yt, y, -5e-5);
1255
1256 Table 10.2 - j and y for orders 3-8.
1257 Compare against excerpts of Table 10.2, Abramowitzh and Stegun.
1258
1259 Important note: In A&S, y_4(0.1) = -1.0507e+7, but Octave returns
1260 y_4(0.1) = -1.0508e+07 (-10507503.75). If I compute the same term using
1261 a series, the difference is in the eighth significant digit so I left
1262 the Octave results in place.
1263
1264 %!test
1265 %! n = (3:8);
1266 %! z = (0:2.5:10).'; z(1) = 0.1;
1267 %!
1268 %! jt = [[ 9.5185e-06 1.0577e-07 9.6163e-10 7.3975e-12 4.9319e-14 2.9012e-16];
1269 %! [ 1.0392e-01 3.0911e-02 7.3576e-03 1.4630e-03 2.5009e-04 3.7516e-05];
1270 %! [ 2.2982e-01 1.8702e-01 1.0681e-01 4.7967e-02 1.7903e-02 5.7414e-03];
1271 %! [-6.1713e-02 7.9285e-02 1.5685e-01 1.5077e-01 1.0448e-01 5.8188e-02];
1272 %! [-3.9496e-02 -1.0559e-01 -5.5535e-02 4.4501e-02 1.1339e-01 1.2558e-01]];
1273 %!
1274 %! yt = [[-1.5015e+05 -1.0508e+07 -9.4553e+08 -1.0400e+11 -1.3519e+13 -2.0277e+15];
1275 %! [-7.9660e-01 -1.7766e+00 -5.5991e+00 -2.2859e+01 -1.1327e+02 -6.5676e+02];
1276 %! [-1.5443e-02 -1.8662e-01 -3.2047e-01 -5.1841e-01 -1.0274e+00 -2.5638e+00];
1277 %! [ 1.2705e-01 1.2485e-01 2.2774e-02 -9.1449e-02 -1.8129e-01 -2.7112e-01];
1278 %! [-9.5327e-02 -1.6599e-03 9.3834e-02 1.0488e-01 4.2506e-02 -4.1117e-02]];
1279 %!
1280 %! j = sqrt ((pi/2)./z) .* besselj (n+1/2,z);
1281 %! y = sqrt ((pi/2)./z) .* bessely (n+1/2,z);
1282 %!
1283 %! assert (jt, j, -5e-5);
1284 %! assert (yt, y, -5e-5);
1285
1286 Table 10.4 - j and y for various integer orders and arguments.
1287
1288 %!test
1289 %! jt = [[ 8.414709848e-01 4.546487134e-01 -1.917848549e-01 -5.440211109e-02 -5.247497074e-03 -5.063656411e-03];
1290 %! [ 9.256115861e-05 2.635169770e-03 1.068111615e-01 -5.553451162e-02 -2.004830056e-02 -9.290148935e-03];
1291 %! [ 7.116552640e-11 6.825300865e-08 4.073442442e-04 6.460515449e-02 -1.503922146e-02 -1.956578597e-04];
1292 %! [ 5.132686115e-18 1.606982166e-13 1.084280182e-07 1.063542715e-03 -1.129084539e-02 7.877261748e-03];
1293 %! [ 7.537795722e-26 7.632641101e-20 5.427726761e-12 2.308371961e-06 -1.578502990e-02 1.010767128e-02];
1294 %! [ 5.566831267e-43 5.836617888e-34 4.282730217e-22 2.512057385e-13 -1.494673454e-03 8.700628514e-03];
1295 %! [ 1.538210374e-61 1.660978779e-49 1.210347583e-33 8.435671634e-22 -2.606336952e-02 1.043410851e-02];
1296 %! [ 3.615274717e-81 4.011575290e-66 2.857479350e-46 2.230696023e-31 1.882910737e-02 5.797140882e-04];
1297 %! [7.444727742e-190 9.367832591e-160 5.535650303e-120 5.832040182e-90 1.019012263e-22 1.088047701e-02]];
1298 %!
1299 %! yt = [[ -5.403023059e-01 2.080734183e-01 -5.673243709e-02 8.390715291e-02 -1.929932057e-02 -8.623188723e-03]
1300 %! [ -9.994403434e+02 -1.859144531e+01 -3.204650467e-01 9.383354168e-02 -6.971131965e-04 3.720678486e-03]
1301 %! [ -6.722150083e+08 -3.554147201e+05 -2.665611441e+01 -1.724536721e-01 1.352468751e-02 1.002577737e-02]
1302 %! [ -6.298007233e+15 -1.012182944e+11 -6.288146513e+04 -3.992071745e+00 1.712319725e-02 6.258641510e-03]
1303 %! [ -3.239592219e+23 -1.605436493e+17 -9.267951403e+08 -1.211210605e+03 1.375953130e-02 5.631729379e-05]
1304 %! [ -2.946428547e+40 -1.407393871e+31 -7.760717570e+18 -6.908318646e+09 -2.241226812e-02 -5.412929349e-03]
1305 %! [ -8.028450851e+58 -3.720929322e+46 -2.055758716e+30 -1.510304919e+18 4.978797221e-05 -7.048420407e-04]
1306 %! [ -2.739192285e+78 -1.235021944e+63 -6.964109188e+42 -4.528227272e+27 -4.190000150e-02 1.074782297e-02]
1307 %! [-6.683079463e+186 -2.655955830e+156 -1.799713983e+116 -8.573226309e+85 -1.125692891e+18 -2.298385049e-02]];
1308 %!
1309 %! n = [(0:5:20).';30;40;50;100];
1310 %! z = [1,2,5,10,50,100];
1311 %! j = sqrt ((pi/2)./z) .* besselj ((n+1/2).', z.').';
1312 %! y = sqrt ((pi/2)./z) .* bessely ((n+1/2).', z.').';
1313 %! assert (j, jt, -1e-9);
1314 %! assert (y, yt, -1e-9);
1315 */
1316