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