GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
besselj.cc
Go to the documentation of this file.
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 
39 {
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 
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 
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 
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 
165  = x_arg.xfloat_complex_column_vector_value ("%s: X must be a scalar or matrix", fn);
166 
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 
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 
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 
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 
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 
255  = x_arg.xcomplex_column_vector_value ("%s: X must be a scalar or matrix", fn);
256 
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 
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 
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 
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 
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 
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 */
octave_value_list do_bessel(enum bessel_type type, const char *fn, const octave_value_list &args, int nargout)
Definition: besselj.cc:84
#define DO_BESSEL(type, alpha, x, scaled, ierr, result)
Definition: besselj.cc:48
bessel_type
Definition: besselj.cc:39
@ BESSEL_Y
Definition: besselj.cc:41
@ BESSEL_K
Definition: besselj.cc:43
@ BESSEL_H2
Definition: besselj.cc:45
@ BESSEL_H1
Definition: besselj.cc:44
@ BESSEL_J
Definition: besselj.cc:40
@ BESSEL_I
Definition: besselj.cc:42
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:95
octave_idx_type numel(int n=0) const
Number of elements that a matrix with this dimensions would have.
Definition: dim-vector.h:401
octave_idx_type length(void) const
Definition: ovl.h:113
bool bool_value(bool warn=false) const
Definition: ov.h:838
FloatComplexNDArray xfloat_complex_array_value(const char *fmt,...) const
bool isnumeric(void) const
Definition: ov.h:703
bool is_scalar_type(void) const
Definition: ov.h:697
Complex xcomplex_value(const char *fmt,...) const
ComplexNDArray xcomplex_array_value(const char *fmt,...) const
ComplexColumnVector xcomplex_column_vector_value(const char *fmt,...) const
FloatComplex xfloat_complex_value(const char *fmt,...) const
bool is_single_type(void) const
Definition: ov.h:651
FloatComplexColumnVector xfloat_complex_column_vector_value(const char *fmt,...) const
double double_value(bool frc_str_conv=false) const
Definition: ov.h:794
bool islogical(void) const
Definition: ov.h:688
OCTINTERP_API void print_usage(void)
Definition: defun.cc:53
#define DEFUN(name, args_name, nargout_name, doc)
Macro to define a builtin function.
Definition: defun.h:56
void error(const char *fmt,...)
Definition: error.cc:968
void scale(Matrix &m, double x, double y, double z)
Definition: graphics.cc:5846
F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE const F77_INT F77_INT & ierr
F77_RET_T const F77_DBLE * x
Complex besselj(double alpha, const Complex &x, bool scaled, octave_idx_type &ierr)
Definition: lo-specfun.cc:814
Complex bessely(double alpha, const Complex &x, bool scaled, octave_idx_type &ierr)
Definition: lo-specfun.cc:815
Complex besseli(double alpha, const Complex &x, bool scaled, octave_idx_type &ierr)
Definition: lo-specfun.cc:816
Complex airy(const Complex &z, bool deriv, bool scaled, octave_idx_type &ierr)
Definition: lo-specfun.cc:131
Complex besselk(double alpha, const Complex &x, bool scaled, octave_idx_type &ierr)
Definition: lo-specfun.cc:817
Complex biry(const Complex &z, bool deriv, bool scaled, octave_idx_type &ierr)
Definition: lo-specfun.cc:1377
std::complex< double > Complex
Definition: oct-cmplx.h:33
std::complex< float > FloatComplex
Definition: oct-cmplx.h:34
octave_value::octave_value(const Array< char > &chm, char type) return retval
Definition: ov.cc:811