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